MONTE CARLO METHOD
The Monte Carlo method is a general algorithm that can be applied to any physical problem that can be broken down into a sequence of events, each of which having a certain probability of occurrence. Before it was applied to radiative transfer (RT) problems, the Monte Carlo method had found applications in a diversity of physical and mathematical problems (Hammersley and Handscomb, 1964). It was first applied to calculations of light scattering in the atmosphere by Collins and Wells (1965) and Plass and Kattawar (1968). Following that, it was applied to light scattering in a coupled atmosphere-ocean system by Plass and Kattawar (1969, 1972).
As a stochastic method, the Monte Carlo RT method does not solve the radiative transfer equations directly. Instead, it tracks the transport of a huge number of photons in a scattering and absorbing system. A photon is usually characterized by its location and direction of propagation. The Mueller matrix is also included if the polarized light field is of interest. The concept of a photon in the Monte Carlo method differs from that in particle physics because it just represents an energy packet that is transported in the scattering and absorbing system. In this sense, it is usually replaced with the concept of a photon packet that carries an additional characteristic, i.e., the photon weight, which represents the fractional number of photons left after one or more interactions in the medium. For every photon packet, its initial photon weight is set to unity, and is reduced when the photon packet is absorbed by a particle or refracted/reflected by the air-sea interface. For each scattering event, the probability that the photon packet is scattered into a desired detector direction can be determined by the scattering phase function at the scattering position, and the corresponding detector response is calculated accordingly. This process is called estimation. After initiation, the transport of each photon packet includes the following steps.
- A photon packet propagates along its current direction until an interaction occurs.
- If the photon packet hits the air-sea interface before collision with a particle, it is refracted or reflected by the interface. The new propagation direction (and Mueller matrix if desired) is determined by Fresnel formulas. Go back to step 1, otherwise proceed to step 3.
- If the interaction is absorptive (single-scattering albedo ω0 < 1), the fraction of the photon packet energy remaining after the interaction is accounted for by reducing the photon weight by a factor of ω0.
- For each detector of interest, when it is possible that the scattered light can reach the detector, the detector response due to this scattering event is estimated.
- A propagation direction after the scattering event is determined, and the photon packet continues propagating along this direction.
- The above five steps are repeated until the photon weight is lower than a preset threshold value, or the photon packet escapes the scattering system.
In other words, the Monte Carlo method works by emulating the real physical processes involved in the transport of radiation energy, represented by the energy carried by an ensemble of photon packets in a scattering and absorbing system. The statistical assumption behind the Monte Carlo method is that the ensemble average of such photon packet contributions represents the true scattered light field as seen by the detectors. Under this assumption, one can keep tracking the history of a sufficiently large number of photon packets and collecting the contributions to the desired detectors. The averaged contribution is equivalent to the true light field.
Sampling techniques and random numbers are used in the above process to determine how each photon packet moves. Suppose there is a physical quantity x, whose distribution is given by a probability density function p(x), such that its expected value can be written as
where x0 and x1 are the lower and upper bounds of the possible range of x, respectively. The corresponding cumulative probability is then
which satisfies P(x0) = 0 and P(x1) = 1. To sample a value for x, one starts with a random number ξ from 0 to 1, and takes
where P-1 is the inverse function of the cumulative probability in Eq. (2). With a sufficiently large number of random numbers, the population of xsample should satisfy the probability density function p(x) in Eq. (1). For example, based on the Bouguer-Lambert-Beer law of extinction, the probability that the photon packet has an interaction at an optical depth τ along the propagation direction is p(τ) = e-τ. Following Eqs. (1)-(3), at step 1, the distance l that the photon travels before hitting a particle is determined by
where c is the extinction coefficient and 1 - ξ has the same distribution as ξ. At steps 4 and 5, the scattering angle that determines where the photon goes after the scattering event is also given by the inverse cumulative probability function, which is derived from the scattering phase function.
Because of its statistical nature, the biggest challenge for the Monte Carlo RT model is the statistical variation associated with the results. In some cases, for instance, when the scattering phase function is highly anisotropic or the extinction coefficient is too small, the size of the ensemble (i.e., the number of photon packets) required to reach a reasonable error may be extremely large. To improve the computational efficiency, one can rely on several variance reduction techniques, among which the most commonly used is called biased sampling. It is based on the fact that Eq. (1) can be alternatively written as
Where w(x') = p(x')/p'(x') is a weighting function, or bias, that relates two probability functions. In case p(x') is not promising for direct sampling, one can instead use p'(x') as the sampling probability function, and multiply the photon weight by the weighting function w(x'). The weighted average given by this sampling scheme is equivalent to that from the direct sampling using p(x') with no weighting function, and suffers less statistical variation provided that w(x') ≤ 1 for any x'. As an example, when the extinction coefficient c is too small, the collision distance lsample sampled by Eq. (4) may be too large, the result of which is that a large fraction of photon packets will escape the scattering system before a collision occurs. In this case, the collision position can be instead sampled by applying a bias,
where τ0 is the optical depth between the current position of the photon packet and the boundary of the scattering system. Note that w(τ) is the probability of a collision within τ0. This implies that the biased sampling forces a collision to happen before the photon packet escapes the scattering system. Now, we can sample τ using the biased probability function p'(τ) = e-τ /(1 -e-τ0), which gives the sampled distance
With this biased sampling, the photon packet will always have an interaction within the scattering system, which reduces the statistical variation.
The Monte Carlo RT method is inherently a 3D method. A photon packet has to be initiatized at a position with three coordinates. Its propagation is then tracked in a 3D scattering medium. A straightforward approach to obtain the radiation field in a plane-parallel scenario is to sample the horizontal position of the incident radiation as well. However, this approach has proven inefficient. An alternative and more effective approach is to fix the incident radiation at one position in the horizontal domain and make use of translational invariance, i.e., the resultant radiation field Iz(x,y;θ,ϕ) corresponding to a single point of incident radiation at I0(0,0;θ0,ϕ0) is the same as Iz(0,0;θ,ϕ) corresponding to an incident radiation at I0(-x,-y;θ0,ϕ0). Therefore, the spatial integration of the radiation field ∫-∞+∞Iz(x,y;θ,ϕ)dxdy resulting from an incident radiation at I0(0,0;θ0,ϕ0) is equivalent to the radiation field Iz(0,0;θ,ϕ) resulting from an integration of point incident radiation ∫-∞+∞I0(x,y;θ0,ϕ0)dxdy, which is plane-parallel incident radiation.
Most RT methods are limited to plane-parallel systems since these are virtually the only ones for which analytic solutions are available. On the other hand, the Monte Carlo method is much more versatile. First, it can be used to calculate the RT in scattering systems with a variety of geometries, such as the spherical shell geometry (Adams and Kattawar, 1978; Kattawar and Adams, 1978; Kattawar, 1979), 3D cases, and Lidar backscatter studies (Adams and Kattawar, 1996). Second, since the Monte Carlo method keeps track of photon packets, it is straightforward to separate detector responses into various components and study the component of interest. For example, water-leaving radiance can be conveniently separated from the total reflected light field. This is a useful case in ocean color studies from remote sensing. It can also be used for spectroscopy studies, where it is important to know the level of line formation in planetary atmospheres.
Adams, C. N. and Kattawar, G. W., Radiative Transfer in Spherical Shell Atmospheres, I. Rayleigh Scattering, Icarus, vol. 35, pp. 139-151, 1978.
Adams, J. T. and Kattawar, G. W., Polarimetric Lidar Returns in the Ocean: A Monte Carlo Simulation, Proc. Soc. Photo-Opt. Instrum. Eng., vol. 2963, pp. 54-59, 1996.
Collins, D. G. and Wells, M. B., Monte Carlo Codes for the Study of Light Transport in the Atmosphere, Radiation Research Associates, Inc., Fort Worth, 1965.
Hammersley, J. M. and Handscomb, D. C., Monte Carlo Methods, Wiley, Hoboken, NJ, 1964.
Kattawar, G. W. and Adams, C. N., Radiative Transfer in Spherical Shell Atmospheres. II. Asymmetric Phase Functions, Icarus, vol. 35, pp. 436-449, 1978.
Kattawar, G. W., Radiative Transfer in Spherical Shell Atmospheres, III. Application to Venus, Icarus, vol. 40, pp. 60-66, 1979.
Plass, G. N. and Kattawar, G. W., Monte Carlo Calculations of Light Scattering from Clouds, Appl. Opt., vol. 7, pp. 415-419, 1968.
Plass, G. N. and Kattawar, G. W., Radiative Transfer in an Atmosphere-Ocean System, Appl. Opt., vol. 8, pp. 455-466, 1969.
Plass, G. N. and Kattawar, G. W., Monte Carlo Calculations of Radiative Transfer in the Earth’s Atmosphere-Ocean System: I. Flux in the Atmosphere and Ocean, J. Phys. Oceanog., vol. 2, pp. 139-145, 1972.
- Adams, C. N. and Kattawar, G. W., Radiative Transfer in Spherical Shell Atmospheres, I. Rayleigh Scattering, Icarus, vol. 35, pp. 139-151, 1978.
- Adams, J. T. and Kattawar, G. W., Polarimetric Lidar Returns in the Ocean: A Monte Carlo Simulation, Proc. Soc. Photo-Opt. Instrum. Eng., vol. 2963, pp. 54-59, 1996.
- Collins, D. G. and Wells, M. B., Monte Carlo Codes for the Study of Light Transport in the Atmosphere, Radiation Research Associates, Inc., Fort Worth, 1965.
- Hammersley, J. M. and Handscomb, D. C., Monte Carlo Methods, Wiley, Hoboken, NJ, 1964.
- Kattawar, G. W. and Adams, C. N., Radiative Transfer in Spherical Shell Atmospheres. II. Asymmetric Phase Functions, Icarus, vol. 35, pp. 436-449, 1978.
- Kattawar, G. W., Radiative Transfer in Spherical Shell Atmospheres, III. Application to Venus, Icarus, vol. 40, pp. 60-66, 1979.
- Plass, G. N. and Kattawar, G. W., Monte Carlo Calculations of Light Scattering from Clouds, Appl. Opt., vol. 7, pp. 415-419, 1968.
- Plass, G. N. and Kattawar, G. W., Radiative Transfer in an Atmosphere-Ocean System, Appl. Opt., vol. 8, pp. 455-466, 1969.
- Plass, G. N. and Kattawar, G. W., Monte Carlo Calculations of Radiative Transfer in the Earth’s Atmosphere-Ocean System: I. Flux in the Atmosphere and Ocean, J. Phys. Oceanog., vol. 2, pp. 139-145, 1972.