Patent Grant
7593835
This invention relates to a radiative transport band model for application to high spectral resolution simulations.
Calculation of in-band radiances via the line-by-line (LBL) approach is computationally intensive. Not only does it require the calculation of spectral molecular absorption coefficients for the series of pressure and temperature pairs encountered in the specified atmosphere, but it also requires the radiation transport equation (RTE) be solved for each spectral point. Because of the fine structure of molecular absorption, especially at higher altitudes (lower pressures), a very finely spaced spectral grid is required. Spectral step size depends upon scenario and accuracy requirements, but values of between 0.01 and 0.00001 cm−1 is common.
The MODTRAN™ solution to this problem is to use a narrow statistical band model. In this approach, the RTE is solved for a relatively narrow spectral band all at once. A number of physical and statistical approximations are made. The quality of the MODTRAN™ results depend largely upon the validity of these approximations for any specific scenario. Over the past 30 years of MODTRAN™ development, much effort has been focused on fine tuning the model to improve accuracy for the terrestrial atmosphere. An earlier version of MODTRAN is described in U.S. Pat. No. 5,884,226.
The integral RTE defines the line-of-sight (LOS) monochromatic diffuse radiant intensity Iv(Ω0;0) at spectral frequency v along a refracted path in the direction specified by solid angle Ω0:
This equation states that the observed spectral radiant intensity at the sensor (s=0) equals the spectral radiant intensity directed towards the sensor a distance s away along the refracted path attenuated by the foreground path transmittance e−τ
The extinction coefficients consists of absorption bvabs(s) and scattering bvsct(s) components, which are computed as sums over atmospheric molecular and particulate species contributions bv,i(s).
The individual species absorption and scattering coefficients are themselves the product of a cross-section κv,i(s) [units of area] and a density ρi(s) [number per unit volume].
The atmospheric source term Jv(Ωs;s) consists of three distinct components: local thermal emission, single scattered directly transmitted solar (and/or lunar) irradiance and diffuse radiation scattered into the LOS
The first term, the thermal emission, is the product of the Planck blackbody function Bv(Ts) at local temperature Ts and the local emissivity, which is represented as a ratio of absorption to extinction coefficients. This expression for thermal emission assumes local thermodynamic equilibrium (LTE) conditions—an assumption, valid in the lower terrestrial atmosphere, that molecular collisions occur in sufficient frequency to control upper state population densities.
The latter two terms in Eq. (4) are the path scattered radiance. The single scatter solar radiation is the product of the top-of-atmosphere (TOA) solar irradiance Fvsun, the scattering point to sun transmittance e−τ
The phase functions of Eq. (4) are all unit normalized, with dimensions of reciprocal steradians (1/sr). It is somewhat more common to define phase functions as dimensionless with a 4π sr normalization, but MODTRAN™ uses the unit normalization to eliminate all constant factors from the source function equation.
Substituting the source term of Eq. (4) into the integral RTE, the spectral radiant intensity of Eq. (1) can be rewritten as absorption optical depth τvabs(s) and scattering optical depth τvsct(s) integrals:
For band model theory, the key atmospheric attenuation parameter is the transmittance itself, not the optical depth. For that reason, it is convenient to again rewrite the RTE in terms of transmittance, tv(s)=exp [−τv(s)], integrals:
Within MODTRAN™, the boundary radiant intensity term Iv(Ωs;s) for the full LOS path is typically zero unless the path terminates at the spherical Earth surface. There is an option to include grey body emission with no reflected source at the path final altitude. For a down-looking sensor whose LOS intersects the ground, the surface radiance includes both emissive and reflective components:
In this equation, μ′ is the cosine of the zenith angle of incoming radiation measured from the ground, fv(Ωs,Ω′;s) is the surface bi-direction reflectance distribution function (BRDF) in units of reciprocal steradians at ground location s, εv(Ωs;s) is the dimensionless directional emissivity in the LOS path direction at s, and the surface reflected downward flux is computed as the upper hemisphere integral (2π sr) over incoming diffuse radiation. The BRDF defines the reflected radiance resulting from surface illumination as a function of both incoming and outgoing angles, and it is normalized to the surface albedo αv(s):
The directional emissivity, the complement of the directional reflectivity ρv(Ωs), dictates the angular distribution of the surface thermal emission, and is computed directly from the BRDF by the following relationship:
Finally, the Earth viewing radiance is obtained by substituting the surface radiance term, Eq. (8), into the RTE, Eq. (7):
In writing this final set of equations, the explicit dependence of the solid angle Ω on position s has been dropped for notational simplicity and the Dirac delta function δ(Ω−Ωsun) has been introduced to move the surface direct solar illumination term inside the angular integration. There is also an implicit understanding that the boundary term Iv(Ω;s) is zero unless the path actually terminates at the ground, s=sg. Eq. (11) serves as the starting point for the derivation of the MODTRAN™ band model.
The MODTRAN™ Band Model
The MODTRAN™ band model approach is concerned with calculation of narrow spectral bin radiant intensities <I0(Ω)>:
The width Δv of the band model spectral bin should be selected sufficiently narrow to insure no more than a modest variation in all spectral data other than molecular absorption. This definition enables all the spectrally smoothly varying properties to be modeled with an effective band value, computed either as an interval average or central value. A smaller band width typically produces finer accuracy. The reformulated band model incorporated into MODTRAN™, described in the parent application which is incorporated herein by reference, provides 4 band model resolution options, 0.1, 1.0, 5.0 and 15.0 cm−1. With an over-line used to represent an effective band value, the band model complement to the LBL LOS radiance expression, Eq. (11), has the following form:
In these equations, the LOS integrands are expressed as functions of the refractive path length s′. Furthermore, the attenuation resulting from absorption is partitioned into its molecular and particulate contributions. The labeling of the band-averaged absorption as arising solely from particulate sources is somewhat misleading. Detailed line compilation data is unavailable for some molecular species, especially large molecules such as the chiorofluorocarbons (CFC's). For these species, absorption is modeled via relatively coarse spectral resolution cross-section data. The resulting transmittances are band averaged values and included in the
Eq. (13) is the fundamental narrow band model radiance equation. It explicitly partitions the spectral data for which effective band values are used from the highly spectrally varying molecular absorption contributions. The effective coarse spectral resolution band data falls into three categories.
The one term for which spectral averaging is somewhat suspect is the TOA solar irradiance. The sun is essentially a ˜6000K blackbody; however, molecular absorption within the solar atmosphere generates Fraunhofer structure. Two facts help support the spectral averaging of the solar data approximation: (1) The spectral correlation between the solar spectrum and the Earth atmospheric absorption is well represented as randomly correlated for most spectral bins; and (2) The current state-of-the-art in measurement and calculation of TOA solar irradiances suggests that spectral structure assignments finer than ˜1 cm−1 are unreliable.
Currently, the multiple scattering (MS) algorithms within MODTRAN™ only support specification of surface reflectance via an angularly independent spectral albedo, fv(Ω,Ω′;s)=αv(s)/π. Given this limitation, the surface reflected downward flux is approximated as a product of the directional reflectivity and the hemisphere integrated downward flux, Fv↓(s):
This expression for the reflected downward flux is quite reasonable if the distribution of downward diffuse radiance is nearly homogeneous or the surface reflectance nearly Lambertian. On the other hand, for clear sky solar illumination of a specular surface such as still water, Eq. (14) is a poor approximation; however, in this case, the observed radiance will be dominated by the direct solar illumination and the error in total observed radiance may still be small. Merging the flux approximation into the narrow band model radiance expression, Eq. (13), the surface reflected diffuse radiance term factors as follows:
When MODTRAN™ utilizes its correlated-k (CK) option, an algorithm that recasts the molecular band model into weighted sum of monochromatic calculations, the path radiant intensity computation follows the prescription of Eq. (13) with the caveat of Eq. (15). However, further approximations are necessary if the band model algorithm is used without CK. The issue here is that the MODTRAN™ atmospherically scattered diffuse radiance and the downward surface flux terms are computed from the MS algorithms within MODTRAN™. These MS algorithms were developed to solve the coupled monochromatic radiation transport problem, and they rely heavily upon Beer's Law. Beer's Law states that the spectral transmittance of a multi-segmented path equals the product of the individual segment transmittances:
tv(0→s′→s)=tv(0→s′) tv(s′→s). (16)
Alternatively stated, segment extinction optical depths are additive. This relationship does hold for CK transmittances but not for band model in-band transmittances:
<t0→s′→s>≠<t0→s><ts′→s>. (17)
When the CK option is not invoked, MODTRAN™ proceeds with multiple scatter source term and level flux calculations computed from a Beer's Law model even though the layer optical depths are poorly defined. The inaccuracy of this approach is limited somewhat by combining the MS source terms with foreground transmittances computed for the inhomogeneous path. Thus, the LOS radiance terms containing diffuse radiation are calculated as follows:
Together, Eqs. (13), (15) and (18) define the MODTRAN™ path radiance formulation. However, implementation of the approach still requires a description of three key model components: the MODTRAN™ statistical band model molecular transmittance calculation, the multiple scattering algorithms, and the surface water model. Descriptions of these components follow.
This invention comprises an improved band model method for modeling atmospheric propagation at arbitrarily fine spectral resolution, sometimes termed herein as “Improvements to MODTRAN4™”. Featured in the invention is a method for predicting atmospheric radiation transport parameters including transmittance and radiance, comprising defining one or more conditions under which the atmospheric transmittance and radiance will be predicted, providing atomic and molecular transition data for a given spectral range and atmospheric conditions, dividing the spectral region into a number of spectral bins that determine the spectral resolution, each bin having a width of less than 1.0 cm−1, calculating atomic and molecular species line center absorption from at least the equivalent width of the atomic and molecular transitions centered within each spectral bin, calculating line tail absorption within each spectral bin from atomic and molecular transitions not centered within the bin, and determining atomic and molecular species spectral transmittances for each spectral bin, the spectral transmittance having a value which is a function of at least the calculated line center absorptions and the calculated line tail absorptions.
The spectral bins may have a width of about 0.1 cm−1. The step of calculating line center absorption may include calculating, from an exact expansion, the bin Voigt equivalent width of atomic and molecular transitions whose centers lie within each spectral bin. The exact expansion may be an exact modified Bessel functions expansion. The calculating line tail absorption step may include subtracting line-tail absorption as calculated from the column strength, the Lorentz half-width, the Doppler half-width, and the line tail spectral displacement. The calculating line center absorption step may include determining the Voigt line-shape function computed at specific frequencies.
The line tail calculation step may include calculating line tail absorption within each bin from atomic and molecular transitions centered outside of the bin using Padé approximant spectral fits to Voigt absorption coefficient curves. The line tail absorption calculation step may include determining a database of temperature and pressure dependent Padé approximant spectral fits to Voigt absorption coefficient curves. There may be five Padé parameters. The Padé parameters may be determined from summed line tail spectral absorption coefficients. One Padé parameter may be determined at the center of the bin, and one at each edge of the bin. One Padé parameter may be the derivative of the absorption coefficient with respect to the normalized spectral variable at the line center. One Padé parameter may be the integral of the spectral absorption coefficient over the spectral band. The Padé parameters database may be generated for a plurality of temperatures. The Padé parameters database may be determined for a plurality of pressures. The line center absorptions may be calculated from atomic and molecular transitions centered no more than half a spectral bin width from the bin, and the line tail absorptions are calculated from atomic and molecular transitions not centered within a half spectral bin from the bin.
The conditions may comprise atmospheric species for which transition data is provided. Band model data may be provided for at least some of the species. Absorption cross section data may be provided for at least some of the species. The conditions may comprise an arbitrary quantity of species for which transition data is provided. The conditions may comprise a surface water model. The surface water model may introduce a spectral attenuation layer. The surface water model may predict surface leaving radiances due to surface water.
The defining step may comprise allowing the user to define the radiation refractive path through the atmosphere. The refractive path may comprise a circular arc. The conditions may comprise the aerosol spectral optical properties. The properties can preferably be varied as a function of altitude. The properties may comprise spectral extinction, spectral scattering and spectral phase function. The method may further comprise determining the medium embedded diffuse transmittance.
Featured in another embodiment of the invention is a method for predicting atmospheric radiation transport parameters including transmittance and radiance, comprising defining one or more conditions under which the atmospheric transmittance and radiance will be predicted, wherein the conditions comprise an arbitrary number of atmospheric species for which transition data is provided, a surface water model, a user-defined radiation refractive path through the atmosphere, and the aerosol spectral optical properties, providing atomic and molecular transition data for a given spectral range and atmospheric conditions, dividing the spectral region into a number of spectral bins that determine the spectral resolution, each bin having a width of less than 1.0 cm−1, calculating atomic and molecular species line center absorption from at least the equivalent width of the atomic and molecular transitions centered within each spectral bin, calculating line tail absorption within each spectral bin from atomic and molecular transitions not centered within the bin, determining atomic and molecular species spectral transmittances for each spectral bin, the spectral transmittance having a value which is a function of at least the calculated line center absorptions and the calculated line tail absorptions, and determining the medium embedded diffuse transmittance.
A primary use of the invention is spectral exploitation. More specifically, this invention greatly enhances the ability to perform the radiation transport necessary to accurately atmospherically correct airborne and satellite imagery of the Earth. Thus, a user can accurately translate the at-sensor observed radiances (i.e. an image captured by an airborne or space-based image sensor) into at-ground spectral reflectances and emissivities. When accurately determined, the spectral exploitation scientists can use the retrieved surface reflectances to classify and identify materials on the Earth. Applications are many including tracking of crop health, monitoring of the littoral zone, mapping of mineral fields, and detecting military targets.
Other objects, features and advantages will occur to those skilled in the art from the following description of the preferred embodiments and the accompanying drawings in which:
IMPROVEMENTS TO MODTRAN4™ Molecular Transmittance
SERTRAN, the Spectral Enhanced Resolution band model IMPROVEMENT TO MODTRAN4™, is statistical. SERTRAN is described in detail in the parent patent application incorporated herein by reference. MODTRAN™ does not account for the exact location of molecular transitions and the distribution of their line strengths within each spectral bin. Instead, these distributions are modeled statistically. In this spirit, the formulation of the molecular transmittance calculation begins with the stipulation that the spectral bin molecular absorption from distinct species is randomly correlated, so that the combined species molecular transmittance <tmol
This approximation works well for a number of reasons. As illustrated in
Within MODTRAN™, the individual species spectral bin transmittance is itself partitioned into three distinct components. These components are:
The starting point for the calculation of the MODTRAN™ in-band molecular line center transmittance is the Plass expression. Plass proved that the homogeneous path transmittance resulting from n identical molecular lines randomly distributed within each spectral interval of width Δv for a contiguous collection of spectral bins is given by the following power law expression:
where Wsl is the finite-bin equivalent width (EW) of a single Voigt line,
The Voigt spectral line-shape fv(γc,γd) [cm] at line center displacement frequency v, the convolution of Lorentz and Doppler line-shape functions, characterizes the spectral dependence of molecular line absorption. In these equations, γc is a collision (Lorentz) half-width at half maximum [cm−1], γd is a Doppler half-width at 1/e of the maximum [cm−1], S is a line strength [cm−2/Atm] and u is a path column density [Atm cm].
The Plass expression, Eq. (21), defines the band model curve-of-growth (COG) for the randomly distributed collection of identical molecular lines, i.e., the rate at which the molecular absorption increases with column density. The actual distribution of molecular lines in a specific spectral interval will exhibit a COG unique to the spectral band. The objective of band model theory in general is to find an optimal fit to each COG. In MODTRAN™, the Plass expression is used as the COG functional form. The fitting parameters are the average line strength of each line S, the effective number of molecular lines in the spectral interval n, the spectral bin averaged Lorentz half-width γc and the Doppler half-width γd.
Band Model Parameter Determination
The traditional band model equations for the temperature dependent line strength, S(T), and effective line number, n(T), parameters are derived by fitting the EW COG in the weak-line (small column density) and Lorentz strong-line (large column density) limits. The Voigt EW varies linearly with column density in the weak-line limit, and it increases with the square root of the column density in the strong line limit. The accuracy of the interpolation between these limits largely dictates the accuracy of the band model. The expressions for the band model parameters resulting from the weak and strong line limit constraints are:
(Note that the actual temperature dependent parameters stored by MODTRAN™ are (S/d)≡nS/Δv and (1/d)≡n/Δv.) Typically, these sums would be over all the transitions centered within the spectral bin for a given molecular species. For MODTRAN™, the spectral bin sums are over spectral sub-intervals of width 0.01 cm−1, and each Sj(T) is itself the sum of temperature-dependent line strengths for lines centered in sub-interval j. This coagulation of molecular transitions within narrow spectral bins has no effect on the calculation the total line strength, n(T) S(T), but it reduces the effective number of lines, n(T). MODTRAN™ temperature dependent band model parameters are generated for 6 atmospheric temperatures, T=180, 205, 230, 255, 280, and 305K. The Lorentz half-width band model parameter, γc0, is computed as a line strength weighted average of molecular transition half-widths at standard temperature (T0=273.15K) and pressure (P0=1 Atm).
The Doppler half-width at 1/e of maximum is defined using the spectral bin center frequency,
γd=√{square root over (2kT/m)}
where k is the Boltzmann constant (1.3807E-16 erg/K), m is the molecular mass (g/molecule), and c is the speed of light in a vacuum (2.9979246E+10 cm/sec).
Sub-intervals are used in the definition of the band model parameters, Eq. (23), for two reasons: (1) the atmospheric absorption from nearly coincident lines (such as those arising from degenerate transitions) is accurately modeled as a single line of combined strength, and (2) the Plass expression was derived for a random distribution of line centers, but near degeneracy occurs at a higher frequency than a random distribution would predict. On a more practical level, it has been shown by extensive comparisons to LBL calculations that the classical expression for the effective number of lines n is too large; too much absorption is predicted. The sub-interval summing decreases the effective line number without changing the total line strength nS, and thereby lowers the band model absorption.
To increase accuracy in the IMPROVEMENTS TO MODTRAN4™ band model absorption for the coarser resolution band models (Δv=1.0, 5.0 and 15.0 cm−1), the effective line number band model parameter is tuned using the 0.1 cm−1 band model. The two-step approach is relatively straightforward. First, the column amount U is determined for which the 0.1 cm−1 band model predicts a homogeneous path transmittance equal to 1/e when integrated (summed) over a given coarse spectral bin:
Here Δv is the coarse band model bin width and the loop is over the N 0.1 cm−1 spectral sub-bins within the coarse spectral bin. Once the spectrally dependent 1/e transmittance column density has been determined, the coarse resolution effective line number n is perturbed to give the same in-band transmittance:
<tk
These calculations are performed at each of the MODTRAN™ temperatures and for a specific pressure depending upon the molecule. For H2O, which is concentrated near the ground, the reference pressure is set to 1 Atm. The absorption of O3 is most important at higher altitudes, so its effective line number values are tuned for a pressure of 0.1 Atm. For all other species, the homogeneous path is defined with a reference pressure of 0.5 Atm.
Calculation of Equivalent Width
The precise spectral domain of the finite bin EW integral in Eq. (22) was left ambiguous. In traditional band model theory, the finite bin EW is computed for a line located in the center of the spectral bin; thus, the spectral domain of the integral extends from −Δv/2 to +Δv/2. However, MODTRAN™ seeks to determine the EW of lines whose centers are randomly distributed within the spectral bin. For a single pure Lorentzian line in the weak line limit, the exact expression for the offset from bin center D(rc) that produces a finite bin integrated absorption equal to the random distribution value can be derived as a function of the ratio of the Lorentz half-width to the bin width, rc=γc/Δv:
The function D(rc) is monotonic for positive rc, decreasing from D(0)=½Δv(the edge of the spectral bin) to D(rc→∞)=Δv/√12. The terrestrial surface Lorentz half-widths are typically ˜0.07 cm−1. For the MODTRAN™ band model widths of 0.1, 1.0, 5.0 and 15.0 cm−1, the 0.07 cm−1 Lorentz transition weak-line displacement values are D(0.07/0.1)=0.028 cm−1, D(0.07/1.0)=0.34 cm−1, D(0.07/5.0)=1.96 cm−1, and D(0.07/15.0)=6.2 cm−1, respectively.
The weak-line Lorentz transition offsets are only a guide for determining the optimal Voigt line-shape transition values; the actual displacements used in MODTRAN™ are D[0.1 cm−1]=0.025 cm−1, D[1.0 cm−1]=0.30 cm−1, D[5.0 cm−1]=1.97 cm−1, and D[15.0 cm−1]=6.3 cm−1 with the EW integral having the form:
The second line above is derived by noting that fvis symmetric in spectral displacement frequency v(fv=f−v), and the Voigt line tail EW WΔ is defined in the last line. The methodology employed to compute the Voigt line tail EW is described in the parent SERTRAN patent application, incorporated by reference herein.
Spherical Refractive Geometry
From its inception, MODTRAN™ included a spherical refractive geometry package. The modeling of refraction was introduced to model near-horizon sun and/or low-altitude near-horizontal path scenarios. As is well known, the sun is already below the horizon when sunset begins. The sun only appears as if it is above the horizon because of refraction.
With the locally spherical Earth approximation used in MODTRAN™, Snell's Law dictates that the product of three terms, (1) the index of refraction, 1+N
[1+N
Here, μs is the cosine of the path zenith at path length s. The MODTRAN™ refractivity, N
For a path specified by an initial altitude H0 and direction μ0<0, MODTRAN™ follows the line-of-sight (LOS) trajectory through the atmosphere continually updating the path zenith angle according to Eq. (1). With normal refraction, the path bends downward from the straight-line path because of the decrease in refractivity with increasing altitude. A limb path is defined as a downward LOS which passes through a minimum tangent altitude, Htan, before beginning an upward trajectory. At the tangent altitude, the path zenith angle is exactly 90°:
[1+N
The MODTRAN™ refractive geometry package was specifically designed to model this scenario.
When large temperature and/or H2O gradients are present within the atmosphere, non-standard refraction can occur. For optical frequencies, this scenario can arise within the marine boundary layer. If refractivity gradient, dN
To properly model the subrefraction and superrefraction scenarios, the invention (IMPROVEMENTS TO MODTRAN4™) provides an option to bypass the built-in refractive geometry package and explicitly enter refractive path data. Specification of a user-defined path requires that some stringent rules be followed. Four columns of data are required. As illustrated in Table 1 (This example is kindly been provided by Drs. Denis Dion and Vincent Ross, Defence Research and Development Canada, Valcartier, 2005), these columns contain the cumulative (spherical) Earth surface distance in meters, the altitude in meters, the path zenith angle in degrees, and the incremental path segment length in meters. In Table 1, all lines beginning with an exclamation are ignored and the word “NEXT” is used as the delimiter. Every atmospheric profile altitude level between the beginning of the path (H1) and the end of the path (H2) must be explicitly included in the user-defined path file. If the LOS passes through a minimum tangent point, the altitude of the tangent point and the altitude of the minimum of H1 and H2 must be included in the path definition. Similarly, if the LOS passes through a maximum tangent point, the altitude of the tangent point and the altitude of the maximum of H1 and H2 must be included in the path definition. In Table 1, a minimum tangent altitude occurs at 0.16075426526368 m, and the final altitude (H2=25.821244410585 m) is listed as the third entry in the downward leg. If a MODTRAN™ atmosphere includes a non-zero ground altitude or a cloud, it is best to run MODTRAN™ first without a user-defined path to determine the model relayering. Note that it is not necessary that “mirrored” path segments within a common atmospheric layer have identical path lengths or angles.
The full suite of refractive path input requirements do make the use of the user-specified refractive geometry option rather complex, but this is currently necessary for modeling strong index of refraction gradients such as those observed in the marine boundary layer.
Inhomogeneous Paths
This line center transmittance section began with a description of the Plass expression for the transmittance of a homogeneous path and the use of this expression by MODTRAN™ to estimate line center transmittances as a function of both column density and the temperature and pressure dependent band model parameters (pressure for the Lorentz half-width). However, the molecular transmittances of the band model RTE, Eq. (13), are defined for an inhomogeneous path. More specifically, the molecular spectral cross-sections κv,iabs(s) and molecular densities ni(s) are functions of path length s along the line-of-sight, and atmospheric conditions change as the path is traversed. In MODTRAN™, this disconnect is rectified by modeling the Earth atmosphere as stratified, with molecular and particulate densities, total pressure, and ambient temperature defined as a function of altitude (there is an option in MODTRAN™ to alternatively define temperature and density profiles as a function of pressure with altitudes determined from the hydrostatic equation and a user defined ground altitude). The layering is assumed to be relatively finely gridded, so that a homogeneous path can be defined within individual layers.
Implementation of the IMPROVEMENTS TO MODTRAN4™ user-specified refractive geometry option required that a completely new approach for calculating MODTRAN™ column densities be incorporated into the model. In MODTRAN4™, column density calculations are performed coincident with the refractive geometry calculation; as the refractive path is incremented in small steps, ds, and the cosine of the path zenith is updated, column amounts are summed by modeling the short incremental segments as a straight lines. With the user-specified refractive geometry option, the full path length within a layer is provided, but not the detailed mapping of path length as a function of altitude within the layer. Like earlier MODTRAN™ versions, IMPROVEMENTS TO MODTRAN4™ assumes an exponential variation of densities and pressure vertically between altitude levels and a linear variation of temperature vertically between altitude levels. However, IMPROVEMENTS TO MODTRAN4™ also models the refracted path between altitude levels as a circular arc, typically with an extremely large radius of curvature—much larger than the Earth radius for normal refraction:
In these path integrals, s is the path length variable; sσ−1 and sσ are path lengths to the beginning and end of segment σ, H(s) is the altitude at path length s; lσ is index of the vertical layer containing segment σ, Hl
If Beer's Law were applicable, the LOS molecular line center transmittance could be computed by simply stringing together individual segment transmittances for each species. However, as detailed above, the in-band absorption does not retain a linear dependence on column density as that density increases. The Curtis-Godson (CG) approximation provides a solution to the band model inhomogeneous path problem. The approach involves approximating the inhomogeneous path by a quasi-equivalent homogeneous path. First, a total path in-band weak-line optical depth (τcen)
This value is used to normalized a CG path-averaged effective line number,
The CG single-line optical depth parameter
Finally, the CG Lorentz and Doppler line-width parameters are defined as the optical depth weighted sum of the n line widths in each band normalized by the product of the weak-line optical depth and the CG effective line number:
In MODTRAN™, the CG data is inserted into the Plass expression to compute line center transmittance for the multi-segmented paths.
The CG approximation works well when the bulk of the absorption occurs within a sub-section of the LOS over which the line shape function does not vary drastically. For most ambient atmospheric species, this condition is satisfied independent of the particular LOS. An important exception for the terrestrial atmosphere is ozone. The absolute O3 concentration has a relatively constant profile from the Earth's surface up into the stratosphere. However, the character of Voigt line shape changes from Lorentzian to Doppler over this altitude range. Thus, MODTRAN™ models the O3 in-band absorption with a higher-order CG correction.
Applying CG path averaging for each of the solar illumination paths arising in the spherical Earth single scatter solar calculation was deemed computationally too expensive when these calculations were introduced into MODTRAN™. As illustrated in
These single illumination path layers are then coupled to the LOS using CG path averaging to calculate the L-shaped path transmittances.
Line Tail Transmittance
The MODTRAN™ spectral bin line tail transmittance is the attenuation arising from molecular transitions centered outside of the bin but no more than 25 cm−1 from the bin center. Unlike the line center spectral attenuation, the line tail absorption varies smoothly across each spectral bin with at most a single minimum. Thus, the entire line tail spectral absorption cross-section curve can be accurately fit to temperature and pressure dependent Padé approximants. This fitting procedure is described in detail within the parent patent application.
For the IMPROVEMENTS TO MODTRAN4™ invention, the Padé parameter database is generated at pressures of 1.0 and 0.1 Atm and at the 6 evenly spaced MODTRAN™ temperatures between 180 and 305K. The Padé parameter interpolations are linear in temperature and linear in pressure-squared. The P2 dependence arises because of the γc2 term in the denominator of the Voigt line-shape, Eq. (22).
An IMPROVEMENTS TO MODTRAN4™ parameter defines the number of monochromatic line tail absorption cross-section calculations per spectral bin so the line tails can be computed at arbitrarily fine spectral resolution. Nominally, the parameter is set to 5. The monochromatic line tail optical depths are added for multiple segment paths since Beer's Law applies, and these summed optical depths are modeled as varying linearly between spectral points for computing the in-band line tail transmittances.
Line Continuum Transmittance
The primary atmospheric molecules for which lines tails beyond 25 cm−1 from line center sum to make a significant contribution to atmospheric absorption are H2O and CO2. MODTRAN™ adopts the continuum algorithms employed in standard LBL models for both of these molecules.
The temperature dependent H2O continuum spectral absorption cross-section κvH
The spectral foreign (CvH
The form of the CO2 continuum spectral absorption cross-section is simpler than that of H2O. CO2 has a zero dipole moment, so partitioning into foreign and self broadened components is unnecessary:
The CO2 data is tabulated on a 10 cm−1 grid from 0 to 10,000 cm−1. Since the CO2 continuum cross-section values are themselves not functions of mole fraction, they can and are directly incorporated into the CO2 line tail band model parameters.
Auxiliary Species
All versions of MODTRAN™ compute the absorption and emission resulting from a standard set of atmospheric molecules. For MODTRAN4™, this standard set includes 26 molecules. However, the invention (IMPROVEMENTS TO MODTRAN4™) provides a new auxiliary species option which couples an arbitrary number of additional molecular species into MODTRAN's radiation transport. These molecules can be modeled either as Beer's Law absorption cross-section species or via the IMPROVEMENTS TO MODTRAN4™ band model. The band model provides greater accuracy if optical depths are substantial, while cross-section data is generally more readily available. The IMPROVEMENTS TO MODTRAN4™ database includes band model parameter data for 24 HITRAN species (OH, HF, HCl, HBr, HI, ClO, OCS, H2CO, HOCl, N2, HCN, CH3Cl, H2O2, C2H2, C2H6, PH3, COF2, H2S, HCOOH, HO2, O, NO+, HOBr and C2H4) not among the MODTRAN™ set of standard molecules (actually, N2 is one of the MODTRAN™ standard atmospheric molecules, but it is modeled in MODTRAN™ as a cross-section species, not as a band model molecule).
In the IMPROVEMENTS TO MODTRAN4™ invention, auxiliary species molecular densities can be defined via user-specified atmospheric vertical profiles or with ‘default’ model profiles; currently, profiles are provided for 16 of the new band model species (OH through PH3 in the above list). The molecular species temperature-dependent optical data inputs must conform to MODTRAN™ definitions. Thus, if molecular band model data are user-supplied, this data must be defined in accordance with the MODTRAN™ band model formalism [Eqs. (23), (24) and (41)]; similarly, molecular absorption cross-section data must be defined at the 6 MODTRAN™ temperatures (180, 205, 230, 255, 280 and 305K) and in units of cm2/molecule.
The mechanism for introducing auxiliary species is similar to that for MODTRAN™ standard atmospheric molecular species when user-defined profiles or radiosonde data are being entered. For each atmospheric altitude profile level, the density of each absorption species is specified. Unlike the standard molecules, the number of auxiliary species is read from the profile table header. This number is a TOTAL number of auxiliary species in the user's master list. The input to the IMPROVEMENTS TO MODTRAN4™ invention accepts a negative sign “−” prefix for any species name to indicate the species is to be excluded from the current calculation. This input format was instituted to avoid having to completely reconfigure an input stream each time a species is added or removed. A star (*) at the end of species name indicates band model rather than absorption cross-section data is to be used. The optical data for each species is stored in species-dependent auxiliary input files (SF6.BM for a SF6 absorption cross-section data file, SF6.cBM for a SF6 line center band model data file, and SF6.tBM for a SF6 line tail band model data file).
To illustrate the auxiliary species option, IMPROVEMENTS TO MODTRAN4™ Top-Of-Atmosphere (TOA) nadir view calculations were run with a 1 km thick surface layer of HCl embedded into the Mid-Latitude Summer (MLS) model atmosphere. The auxiliary species 0.1 cm−1 band model data files, p1_HCl.cBM and p1_HCl.tBM, were read in. Calculations could have been performed with the default HCl profile, but the current objective was to predict the change in spectral signatures with increasing HCl column density. HCl densities of 0.1, 1.0 and 10.0 ppm were input. The transmittance spectra for the HCl component are illustrated in
The transmittance spectra for the HCl fundamental R2 and P3 vibrational/rotational lines with the HCl surface layer embedded in the MLS atmosphere are shown in
The thermal emission spectra are also plotted in
Spectrally Varying Aerosol Profiles
MODTRAN™ has always included an option to define up to 4 aerosol types, each specified in a distinct altitude region. Nominally, a boundary layer aerosol density profile is defined for altitude levels up to 2 km, a tropospheric aerosol density profile is defined for altitude levels from 3 to 10 km, a stratospheric aerosol density profile is defined for altitude levels from 11 to 30 km, and a high altitude aerosol density profile is defined above the stratopause. Aerosol spectral optical properties within each of these regions do not vary with altitude; the aerosol concentration profile provide the only variation with altitude. A new Spectral Aerosol Profiles (SAP) option has been developed for the IMPROVEMENTS TO MODTRAN4™ invention. The SAP option allows aerosol spectral extinction coefficients, spectral scattering coefficients, and spectral scattering phase functions all to be entered as a functions of altitude. This enables the IMPROVEMENTS TO MODTRAN4™ invention to more realistically capture the true effect of aerosols on atmospheric radiation transport. As part of this development, the Legendre expansion representation of the SAP aerosol phase function is directly coupled into the DISORT multiple scattering (below). In all the other MODTRAN™ aerosol options, default or user-specified, the scattering phase function is always approximated by a Henyey-Greenstein parameterization for use in multiple scattering computations. This is a critical upgrade for accurate inclusion of aerosols in the scattered radiance calculations.
Whenever the SAP option is exercised within IMPROVEMENTS TO MODTRAN4™, spectral aerosol profile data is read in from an external file, a sample of which is included in Table 2(This example is kindly been provided by Drs. Denis Dion and Vincent Ross, Defence Research and Development Canada, Valcartier, 2005). The first line of the SAP file contains 3 integers: the number of spectral grid points, the number of scattering phase function Legendre moments, and the number of scattering phase function angular points. The second line contains the scattering phase function angular grid in degrees. The next group of lines contains the aerosol data for the first MODTRAN altitude level. There is a separate line for each wavelength. The individual lines contain the altitude level, the spectral wavelength, extinction and scattering coefficients, and phase function data. Legendre moments are listed first followed by phase function values at each of the angular grid points. Subsequent groups of lines contain the data for additional altitude levels.
In the RTE, the absorption and scattering coefficients, Eq. (3), and the effective phase function, Eq. (5), are all express as sums of individual species contribution. When the SAP options is exercised in IMPROVEMENTS TO MODTRAN4™, the SAP aerosol is treated as one of these species. The extinction coefficient, scattering coefficient, and required phase function profile data are spectrally interpolated (linearly in wavelength) at the central bin frequency within the band model spectral loop. The extinction and scattering coefficients are modeled as varying exponentially with altitude, while the scattering phase functions and the Legendre coefficient are linearly interpolated in altitude. The result is a first principles treatment of the aerosol contributions.
Multiple Scattering via the DISORT Algorithm
MODTRAN™ computes LOS spectral radiant intensities by integrating the attenuated spectral source function along the path defined by its spherical refractive geometry algorithm and summing the boundary term if termination occurs at the hard Earth. This is the essence of Eq. (1). However, since the multiple scattering component of the source function [Eq. (4)] and the downward flux component [Eq. (14)] of the boundary term both depend on LOS spectral radiant intensities themselves, there is an internal coupling that complicates the path integration procedure. MODTRAN™ resolves the self-consistency problem by using a plane-parallel multiple scattering algorithm to compute the LOS spectral radiant intensity dependent contributions. Underlying this procedure is the assumption that multiple scattering in the terrestrial atmosphere is predominately a lower atmosphere effect and therefore well characterized with the locally flat Earth approximation. This ansatz is least accurate for low sun conditions. Thus, for solar elevation angles below ˜10°, MODTRAN™ solar radiances carry a relatively high uncertainty.
The MODTRAN™ model includes two multiple scattering algorithms. An analytic two-stream algorithm dubbed the Isaacs approach approximates the multiple scattering source terms with minimal computations. These two-stream results are best suited for calculation of spatially integrated radiances (fluxes) and thermal scatter contributions. However, the accuracy of the 2-stream solar scatter radiant intensities is relatively poor, typically ±20%, but factor of 2 errors are not uncommon in the shortwave end of visible spectral region (blue). For higher fidelity results, MODTRAN™ employs the DISORT discrete ordinate algorithm to compute multiple scatter radiances. In this section, the DISORT algorithm and its IMPROVEMENTS TO MODTRAN4™ interface are described.
Chandrasekhar's Solution to the Plane-Parallel Atmosphere RTE
DISORT explicitly solves the integral-differential RTE for a horizontally-homogeneous plane-parallel atmosphere with no diffuse radiation impingent on the top of atmosphere (TOA) and with an emitting and reflecting ground surface:
The top equation in this set is the path optical depth derivative of the integral RTE, Eq. (1), for a plane-parallel atmosphere. Since refraction and spherical Earth effects are excluded, the path solid angle Ω[≡(μ,φ)] is constant for a given LOS. The reference to position has also been dropped as an argument of the BRDF because the Earth surface is modeled as horizontally homogeneous. For notational simplicity, the spectral frequency subscript has also been dropped from Eq. (46) and all subsequent equations in this section, even though all quantities other than the temperatures and angles are spectrally dependent. Finally, the path optical depth is represented as the quotient of the path nadir optical depth measured from the TOA τ↓ and the cosine of the path zenith angle with the ground value, τg↓, equal to the total atmosphere nadir extinction optical depth. Any remaining references to position s have been replaced by references to the nadir optical depth τ↓ because of the lack of horizontal variability.
In order to facilitate solving the plane parallel atmosphere radiation transport problem, two additional model assumptions are made up front: (1) the scattering phase function is modeled as a function of the angle Θ between the incident and outgoing beams, p(Ω,Ω′;τ↓)=p(Θ;τ↓); and (2) the surface BRDF is modeled as a function of the incoming and outgoing zenith angles and of the relative azimuth angle, but not of the incoming and outgoing azimuth angles independently, i.e., f(Ω,Ω′)=f(μ,−μ′,φ−φ′). These assumptions limit DISORT's ability to precisely model physical phenomena such as preferential orientation of ice particles in cirrus clouds and plowed fields on the Earth's surface. However, under most terrestrial atmospheric scenarios, it is quite reasonable to model suspended particles as randomly oriented and the Earth surface as locally isotropic.
DISORT follows a 4-step approach, originally laid out by Chandrasekhar, to solve the integral-differential RTE. In developing the DISORT model, Stamnes and his coworkers solved the previous unsolved and challenging problem of creating a numerically stable realization of the Chandrasekhar approach.
Step 1: Factoring the Azimuth Dependence
The azimuth dependence of the RTE is factored from the integral-differential RTE by expanding the radiance and bidirectional reflectivity in Fourier cosine series, and expanding the scattering phase function in Legendre polynomials, Pl (μ), and unit normalized associated Legendre polynomials, Λlm(μ):
In these equations, the Im(μ;τ↓) are the azimuth moment radiant intensities, the fm(μ;−μ′) are the azimuth moment BRDFs, the gτ
Note that in the absence of solar radiation, only the m=0 problem needs to be solved. This is also true if there is azimuthal symmetry resulting from either a zenith sun or a vertical/nadir LOS.
Step 2: Introducing Discrete Ordinates
In the second step of the Chandrasekhar procedure, the independent integral-differential equations are each converted into a system of 2N coupled ordinary differential equations (ODE's) by applying the discrete ordinate approximation to the cosine zenith integrals.
DISORT employs a double Gaussian quadrature, with paired cosine angles, μj=−μ−j, symmetric in each hemisphere about μ=½ and having a common Gaussian weight, wj=w−j.
Step 3: Simplifying Optical Depth Dependence via Atmospheric Layering
In the third step of the Chandrasekhar procedure, the optical depth dependence of the atmospheric single scattering albedo, the scattering phase function and the Planck emission is parameterized so that an analytic matrix solution to the system of ODE's can be derived. This is accomplished by partitioning the atmosphere into l=1 to L layers with l=0 to L level boundaries. Altitudes z0=TOA to zL=ground, nadir optical depths τ0↓=0 to τL↓=τg↓, and temperatures T0 to TL are specified at each altitude level. Scattering properties are modeled as constants within each layer:
while Planck emission is characterized as varying linearly with nadir optical depth through each layer:
A general solution can be derived for the m-dependent RTE within each atmospheric layer given this parameterization of optical depth dependence:
In Eq. (52), the Y0l(μi) and Y1l(μi) arise as a particular solution to the thermal source problem and are obtained by solving a system 4N linear equations. Similarly, the Zlm(μi) are a particular solution to the solar source problem, obtained by solving a system of 2N linear equations for each value of m. The kjlm and Gjlm are the eigenvalues and eigenvectors, respectively, that arise in solving the hormogeneous plane-parallel RTE problem [Fsun=0; Bl=0], and the Cjlm are the constants of integration that must be determined from boundary and layer continuity conditions.
Step 4: Applying Boundary and Continuity Equations
For each of the azimuth moments m, the L atmospheric layers each contain 2N constants of integration. These values are determined from the (2N)×L linear boundary and continuity equations:
Solving for these constants is the most computationally time intensive step in DISORT.
Once the constants of integration have been calculated, the specific requested spectral radiant intensities and horizontal fluxes can be computed. The interpolated spectral radiant intensity in direction (±μ, φ) and at an arbitrary nadir optical depth τ↓ within the plane-parallel atmosphere is computed by inserting the integral RTE into the radiant intensity expression of Eq. (47):
The azimuth moments of the source function Jm(±μ;τ↓) are defined by Eq. (49). Substitution of Eqs. (50), (51) and (52) for the nadir optical depth dependent terms yields simple exponential integrals that are evaluated analytically. For the upward and downward flux calculations, only m=0 terms contribute. The fluxes are computed from the Gaussian quadrature integral representation:
In the last line of Eq. (55), the relationship between the total (direct+diffuse) downward flux at the ground and the diffuse downward flux of Eq. (14) is explicitly presented.
The MODTRAN™ Interface to DISORT
In principle, MODTRAN™ could be used solely as the optical property and geometry specifier for DISORT, defining the DISORT inputs as a function of wavelength, and DISORT used as the RT engine, the generator of the RT outputs. However, in this mode, the MODTRAN™ ability to incorporate spherical refractive effects would be lost. As described above, DISORT is instead only used to compute the multiple scattering and surface downward flux contributions to Eq. (13). In this sub-section, the precise interface is laid out. The sub-section begins with a description of DISORT convergence issues, the delta-M method and its impact on the MODTRAN™ interface. The IMPROVEMENTS TO MODTRAN4™ method of defining DISORT inputs is described next, and that is followed by a description of special uses of DISORT via the IMPROVEMENTS TO MODTRAN4™.
DISORT Convergence, the Delta-M Method and Radiant Intensities
Atmospheric aerosol and cloud particulates typically possess highly forward peaked scattering phase functions. Accurate representation of these functions using a Legendre expansion, Eq. (47), requires hundreds of terms. Computation times grow rapidly with the number of discrete ordinate streams, which DISORT couples to the number of azimuth terms.
To alleviate the convergence problem, the delta-M method was introduced. In this method, the scattering phase function is represented as a sum of a forward peaked Dirac-delta function δ(1-cos Θ) and a Legendre expansion:
Here, the delta-M method expansion coefficients g′τ
With the delta-M method, accurate spectral radiant intensities and fluxes can generally be computed with as few as 4, 8 or 16 streams. Substituting of the delta-M expansion for the phase function into the RTE transforms it into an identical equation, but with p, τ↓, and ω replaced by p′, τ′↓, and ω′, where the later two quantities are defined by:
DISORT automatically implements these transformations.
The one obvious drawback of the delta-M method is that the transformed radiant intensities are somewhat inaccurate for very near forward scattering scenarios, such as occurs in the study of the solar aureole, because the finite width of the forward peak is not characterized. To improve near forward scattering accuracy, DISORT provides an option to invoke a series of radiant intensity corrections developed by Nakajima and Tanaka. In these corrections, the plane-parallel atmosphere single scatter solar radiant intensity is explicitly subtracted from the diffuse solar radiant intensity and then added back in using a high order Legendre expansion for the scattering phase function. Additional, less dramatic Nakajima and Tanaka corrections improve higher order solar scattering intensity terms.
For DISORT runs from within MODTRAN™, it is important that the Nakajima and Tanaka corrections not be invoked. MODTRAN™ itself subtracts the plane-parallel single scatter solar contribution from the DISORT at-source solar radiant intensity for each layer l, ΔIl(μ,φ):
In this equation, ΔIlmss(±μ,φ) is the at-source multiple (no single) scatter solar radiant intensity from layer l; ωl and p′l(Θ) are the single scattering albedo and delta-M scattering phase function within layer l; Δτl↓ is the nadir optical depth through layer l; tl→sun is the solar path transmittance from the bottom of layer l; and tl−1→sun is the solar path transmittance from the top of layer l. MODTRAN™ subsequently adds in the spherical refractive geometry single scatter solar radiant intensity. MODTRAN™ computes the species scattering phase functions either from an analytic form (for Rayleigh and Henyey-Greenstein representations) or from a scattering angle dependent table. For user-defined scattering phase functions, both the Legendre coefficients and the tabulated values are required inputs to insure that the Legendre expansion does not have to be relied upon for the single scatter calculation.
MODTRAN incorporates the DISORT multiple scattering layer contributions into its LOS radiant intensity by first defining a DISORT thermal and multiple scattering solar layer source terms, ΔJlthm and ΔJlmss:
These source terms are then multiplied by refractive path transmittance differences and added to the single scatter solar contributions:
This procedure enables the band model molecular transmittance COG to be incorporated into the multiple scatter calculation.
DISORT I/O for MODTRAN™
The complete list of inputs to the plane-parallel multiple scattering algorithms in MODTRAN™ is quite modest. As noted above, the atmosphere is stratified with each layer l characterized by a nadir optical depth, Δτl↓, and layer averaged values for the single scattering albedo, ωl, and the scattering phase function Legendre moments, gll. In addition, layer boundaries are characterized by altitude level temperatures, Tl. Other than this profile data, the only inputs are spectral range (required for the Planck function evaluation), TOA solar irradiance, solar angle and observer zenith and relative solar azimuth angles.
In principle, solving for layer source functions, ΔJl, is simplest in the optical thin limit. However, this limit is occasionally problematic for DISORT. To simplify the source code, the developers of DISORT made a conscience effort to avoid branching for special cases. As a result, inclusion of layers with a very small nadir optical depth, Δτl↓, can produce numerical instabilities. This problem is avoided by merging MODTRAN™ layers to insure that Δτl↓ exceeds a 0.0001 threshold. This actually provides the added benefit of speeding up the DISORT calculations, because computation time is proportional to the number of computation layers. The DISORT output of vertical flux profiles are redefined onto the finer gridded MODTRAN™ levels via optical depth weighted interpolations. The coarse layer source function values are assigned without modification to the MODTRAN™ sub-layers with one caveat: the thermal source for thin layers (nadir optical depth of less than 0.001) is set to the layer blackbody emission value (the conservative scattering approximation).
The calculation of the layer source functions is numerically more stable if the solar and thermal DISORT scatter calculations are performed independently, especially for azimuth dependent scenarios. One reason for this is that when DISORT solar scatter calculations are performed and the nadir absorption optical depth is large, DISORT only calculates the solar scattering contributions for the upper atmosphere; atmospheric layers for which the nadir absorption optical depth exceeds 10 are excluded. However, the dark layers are not excluded from the solar calculations if thermal scatter contributions are computed simultaneously. Combining of solar and thermal calculations results in numerical instabilities. Thus, whenever MODTRAN performs azimuth dependent scattering DISORT calculations, the solar and thermal contributions are computed from separate calls to DISORT.
Additionally, special cases can arise that are problematic for DISORT and either require special handling by the MODTRAN™ interface or complete avoidance. Three of these cases involve a solar geometry. If the solar zenith angle is too close to one of the Gaussian double quadrature angles, the IMPROVEMENTS TO MODTRAN4™ package dithers the solar angle to insure that |μi−μsun|>μsun/10,000. If the view zenith angle is near coincident with the solar zenith angle, the view angle is perturbed as follows:
where μ0 is the original view cosine zenith and μ0′ is the view cosine zenith used in the DISORT algorithm. Finally, if the sun is below the horizon (μsun≦0), then solar multiple scattering contributions are set to zero. Spectrally, if the single scattering albedo is too close to unity for any layer, it is lowered to a value that will insure a numerically stable solution; this can lead to inaccurate scattered radiances for a purely Rayleigh atmosphere.
MODTRAN™ Atmospheric Correction Data File
The observed solar scatter radiance at nadir optical depth τ↓ measured from an Earth viewing Visible/Near-Infrared (VNIR) through Shortwave-Infrared (SWIR) hyperspectral imager along the Ω0 direction can be conveniently partitioned into three components: the scattered only by atmosphere radiance Iatm(Ω0;τ↓) the directly transmitted from ground reflected radiance Idirgnd(Ω0;τ↓), and the diffusely transmitted from ground reflected radiance Idifgnd(Ω0;τ↓):
As before, the spectral index has been omitted for notational simplicity. This partitioning is particularly useful for remote sensing applications because it provides a distribution of ground information content: Iatm contains no information about the ground; Idirgnd is directly proportional to the imaged pixel reflectance; and Idifgnd provides data about surrounding ground pixels.
For a specified atmosphere, the scattered only by atmosphere solar radiance is computed in DISORT by setting the upper illumination at the lower boundary, Im(+μ;τg↓), to zero and dropping the thermal contributions from Eq. (46):
In practice, the lower boundary condition is specified by setting a thermal flag to false and setting the surface reflectance to zero.
The directly transmitted from ground reflected radiance, Idirgnd can be expressed as a function of ground-dependent only and atmosphere-dependent only terms if the imaged pixel is modeled as a Lambertian reflector with BRDF, f(+μ,−μ′), equal to α/π:
In this expression, the product μsun Fsun equals the top-of-atmosphere (TOA) solar downward horizontal flux. Multiplying this product by the sun-to-ground total (direct+diffuse) transmittance, ttotalsun→gnd (=tdirectsun→gnd+tdiffusesun→gnd) results in the expression for the surface downward solar flux given a black (i.e., non-reflecting) and cold (i.e., non-emitting) ground. If the ground is actually a Lambertian reflector with a positive spatially-averaged surface albedo, i.e., <α>>0, then the product <α>σgnd equals the fraction of initial downward illumination at the ground that will return to the bottom of the atmosphere after a single ground reflection, where σgnd is the atmospheric spherical albedo from the ground. The atmospheric spherical albedo is defined as the fraction of lower boundary isotropic illumination that will scatter back to the ground. The multiple reflections expression is a geometric sum:
Thus, the left-most quotient on the right-hand side of Eq. (64) equals the total downward solar flux at the ground. A final multiplication of this surface flux by the Lambertian reflectance and the direct transmittance to the observer, tdirectobs→gnd, provides the expression for the directly transmitted ground reflected radiance.
MODTRAN™ itself computes direct path transmittances, tdirectsun→gnd and tdirectobs→gnd. The remaining atmosphere-dependent only terms of Eq. (64) are readily determined using DISORT by appealing to the RT reciprocity principle. As noted by Stamnes, the plane albedo from and diffuse transmissivity through a vertically inhomogeneous planetary atmosphere can be expressed in terms of normalized emergent diffuse intensities from isotropic incident illumination:
where the plane albedo, σgnd(−μ′), is related to the spherical albedo, σgnd, by the integral equation
σgnd=∫2πσgnd(−μ′)μ′dΩ′. (67)
The beauty of these equations is that there is no azimuth dependence, so only the m=0 solution need be computed, and the particular solution does not have to be calculated because the boundary and the external source and the emitted source terms are all zero. Thus, these calculations are computationally modest.
The final term of Eq. (62) is the diffusely transmitted ground reflected radiance:
In this equation, the total downward solar flux, the left-most quotient on the right-hand side, is multiplied by the spatially averaged Lambertian reflectance and a medium embedded diffuse transmittance term, temb
As this equation illustrates, the radiative transport problem that must be solved to calculate temb
The calculation of spherical albedo, the solar path diffuse transmittance and the view angle medium-embedded diffuse transmittance by DISORT has been streamlined as part of the IMPROVEMENTS TO MODTRAN4™. DISORT provides an option to calculate spherical albedo and diffuse transmittance in a dedicated run of the model. The IMPROVEMENTS TO MODTRAN4™ use DISORT to calculate both the line-of-sight (LOS) scattered radiance contributions and the atmospheric correction data. These two calculations required the same sets of Legendre polynomials and atmospheric layer Eigenvalues and Eigenvectors; they differ in their constants of integration because of differing boundary conditions. DISORT in IMPROVEMENTS TO MODTRAN4™ has been modified to generate both the LOS and atmospheric correction data (including the medium-embedded diffuse transmittance) in a single call with redundant computations eliminated. With a vector optimized compiler, the extra computational effort required to generate the atmospheric correction data during a LOS scattered radiance calculation only amounts to a few percent of the overall time.
Surface Moisture
Atmospheric water vapor retrievals from down-looking visible and infrared imaging sensor data are adversely impacted by surface droplets and/or embedded liquid water on otherwise dry backgrounds, as with dew or within chlorophyll structures like the leaves of a forested canopy. The water vapor retrieval algorithms typically assume the background and atmospheric water signatures are uncorrelated; this is untrue if the background contains liquid water. This spectral correlation of the surface signature with atmospheric water also presents a problem for viewing against extended water body backgrounds, as in littoral scenes. The details of the radiation transport (RT) in the background underwater scenario is considerably different from embedded water or water droplet scenario. This section describes approaches recently implemented within IMPROVEMENTS TO MODTRAN4™ for dealing with both situations.
Embedded Water Attenuation Model
The IMPROVEMENTS TO MODTRAN4™ embedded water or moisture attenuation model incorporates the effects of the liquid water extinction directly into its surface spectral reflectance model. For a dry background, MODTRAN™ parameterizes surface reflectance via a spectral grid λi(i=0, . . . ,N) of bidirectional reflectance distribution functions (BRDFs), fi(Ωv,Ωs) and their angularly integrated products. These reflectance products include (a) ρi(Ωv), the surface spectral directional reflectivities in the view solid angle direction, Ωv≡(μv,φv); (b) ρi(Ωs), the surface spectral directional reflectivities in the source solid angle direction, Ωs≡(μs,φs); and (c) αi, the double hemisphere integrated surface spectral albedos:
As indicated above, the solid angles Ωv and Ωs have two components: the cosine of a zenith angle μ and an azimuth angle φ.
The embedded water attenuation model uses a single semi-empirical depth parameter, D, to incorporate liquid water attenuation into the surface spectral reflectance terms:
In these definitions, bi is the liquid water spectral extinction coefficient (1/m) at spectral wavelength λi. The extinction coefficients are computed as the sum of absorption biabs and scattering biscat terms,
In this expression, λ is the spectral wavelength in microns. As evidenced by
Since MODTRAN™ only needs the BRDF and directional reflectivities values at the solar and (downward) viewing angles, a generally coarse spectral grid of values are computed up front, i.e., before the primary MODTRAN spectral loop. For the dry surface, the reflectance quantities are then linearly interpolated in wavelength, as needed:
However, the coarse wavelength grids commonly used for the surface description will not capture the spectral structure of the liquid water absorption,
Incorporating the finer spectral dependence of the water absorption is straightforward for the BRDF term:
Modeling of the directional reflectivity and double hemisphere integrated albedo terms is more difficult because the angular integrations, which are preferentially performed in preprocessing, are spectrally dependent on the water absorption. The following approximation is introduced to compute the spectral directional reflectivities:
and E3(bλD) is the exponential integral of order 3:
Here the mean value theorem is being invoked to pull fi(Ωv,Ωs) out of the integrand at Ωs=<Ωs>i and fi(Ωv, <Ωs>i) is being approximated by a transmittance-weighted average at λ=λi. Thus, even though fi(Ωv, <Ωs>i) is actually defined as a spectrally varying quantity because of its dependence on the extinction coefficient bλ, it is approximated by a spectrally independent quantity. This weighted average insures that ρλ(Ωv; D) indeed equals the precomputed integral ρi(Ωv; D) at the spectral grid point, λ=λi. Since the effective hemispheric transmittance equals twice the exponential integral [note that 2E3(biD)|D=0 equals unity, as it should], a coefficient of 2 is included with each occurrence of the exponential integral in the final line of Eq. (75). With this expression for the directional reflectivity, the only integral that has to be computed on-the-fly spectrally is the exponential integral of order 3, and well established algorithms are readily available for this evaluation.
The double integral mean value theorem and the corresponding double path transmittance weighted average are used to derive an analogous expression for the spectral albedo:
Again, the mean value theorem allows that BRDF to be pulled out of the integrand, but fi is evaluated at a pair of solid angles, (Ωs,Ωv)=<(Ωs,Ωv)>i for the surface albedo calculation.
In the preprocessing step, the ρi(Ω;D) of Eq. (75) and αi(D) of Eq. (78) are computed by numerically integrating the BRDF; these integrals are divided by the spectral grid denominator terms [2 E3(biD) exp(−biD/μ) and 4 E2(biD)2, respectively] and stored for use in the fine spectral resolution calculations. No numeric intergrations over the BRDF are performed during this band model spectral loop—only spectral interpolations and calculations of exponentials and exponential intergrals.
Flat Water Layer Model
The flat water layer model is quite distinct form the embedded layer model in that the liquid water is modeled as a layer that lies above a surface with specified reflectance properties. This introduces three effects in addition to the attenuation embedded water model:
In these equations, ri(<1) is the air-to-water index of refraction ratio. In the first line of Eq. (80), the ri2 and ri−2 terms equal the radiance decrease from the expanded solid angle and increase from the constricted solid angle, respectively, resulting from Snell's Law refraction; these terms cancel because the radiation both enters and exits the water. The t(θ′,θ) terms are the spectral grid Fresnel interface transmittances. Primes (′) are used to indicate underwater angles.
Approximate approaches are used to model the initial above water Fresnel reflection and the multiple under water reflections. Since water bodies are seldom well represented as having a strictly flat surface, the Fresnel reflection is averaged over the full hemisphere of incident angles, θs; this approximation smoothes out the large contributions that occur for the specular reflection geometry. The reflected underwater contributions are modeled by treating the upwelling radiation from below the water-air interface is nearly isotropic. The fraction of photons whose incident angle will exceed the critical angle is (1−ri2). Additional photons are Fresnel reflected, but this contribution is small because the Fresnel transmittance is close to one for incident angles less than the critical angle. Thus, a fraction equal to ri2 of the upwelling flux penetrates the air-water interface, and (1−ri2) is reflected. The next question that must be answered is “what fraction of the photons reflected from the underside of the water-air interface will return to that interface?” This water layer bottom spherical albedo term is approximated by the double path attenuated water layer surface albedo αi(D) of Eq. (71). Two somewhat compensating assumptions are made: (1) that the backscatter from within the water layer is relatively minor and (2) that the reflected downward flux is nearly isotropic even though the total internal reflections only occur for larger incident angles. Given these approximations, the flat water layer model BRDF is expressed as
The first two terms in Eq. (83) are illustrated pictorially by
To calculate of the flat water layer BRDF at an arbitrary wavelength, the underwater surface BRDF and the water layer attenuated surface albedo are spectrally interpolated, Eqs. (73) and (76), respectively:
The angularly averaged Fresnel reflection, <1−tλ(θ′s, θs)>, is solely a function of the air-to-water refractive index ratio, rλ. This ratio falls in a range from 1.09 to 1.56 over the entire far ultra-violet through long-wave infrared spectral region from 0.2 to 32 μm; with this limited range, the angularly averaged Fresnel reflection can be accurately fit to a quadratic in rλ for efficient computation:
<1−tλ(θ′s,θs)>≈−0.3421+0.5033rλ−0.1532rλ2 (86)
The spectral grid directional reflectance in the view direction, ρi(Ωv;D;ri), and surface albedo, αi(D;ri), for the flat water layer model are defined by the single and double hemisphere integrated BRDF, respectively:
where {tilde over (ρ)} is defined as the water leaving radiance contribution to ρ. Since the BRDF is defined for the underwater surface, it is easiest to perform numerical integrations using the underwater source and view angles. A change of variables via Eq. (81) leads to the following expressions:
In these expressions, the angular integral domain itself is spectrally dependent.
For the calculation of directional reflectivity at an arbitrary spectral point, ρλ(Ω′;D;rλ), the spectrally slowly varying BRDF assumption used in the embedded water attenuation model is replaced with the ansatz that the product of the BRDF with the Fresnel interface transmittance is spectrally slowly varying:
As illustrated in
Once again, the mean value theorem is invoked; here it is used to extract the [ft]λ product out of the spectral directional reflectivity angular integration:
Like the full hemisphere water layer transmittance, Eq. (77), the conical water layer transmittance of Eq. (91) is a function of elliptical integrals, but with added complexity resulting from the restricted domain:
The spectral grid components of [ft]λ are approximated by the conical water layer transmittance weighted average evaluated at the spectral grid points:
By first substituting Eqs. (90) and (92) into the expression for the spectral directional reflectivity, Eq. (91); then approximating the [ft]i terms according to Eq. (93); and finally incorporating the spectral grid directional reflectivity expression, Eq. (89), the spectral directional reflectivity in the view (source) direction has the following form:
where the denominator Cλ(ρ) is a transmittance enhanced by multiple reflections
A completely analogous procedure is used to define the flat water layer spectral double hemisphere albedo. The water leaving contribution to the spectral grid value has the following fully expanded form:
Spectral interpolations are performed over a triple product, the BRDF with the Fresnel transmittance in both the view and source directions:
The final expression for the flat water layer spectral albedo contains the double path conical transmittance enhanced by multiple reflections:
During MODTRAN™ pre-processing for the flat water layer model, the BRDF is numerically integrated to compute the spectral grid point directional reflectivities ρi(Ω;D;ri) and albedos αi(D;ri). These values are stored after being divided by transmittances terms Cλ(ρ) and Cλ(α), respectively. No numerical integration of the BRDF is required at the finer spectral resolution of the MODTRAN™ band model. The flat water layer model calculations performed at this finer resolution are quick and straightforward; they include calculation of the liquid water and air refractive indices, calculation of the Fresnel transmittance, calculation of the elliptical integrals and calculation of the spectral interpolations.
Special Cases
It is of interest to examine both the embedded water attenuation and the flat water layer models for the scenario in which the underlying BRDF is Lambertian:
fiLamb(Ωv,Ωs)≡αi/π. (100)
For the embedded water attenuation model, substitution into Eq. (71) gives the following expressions:
Not surprisingly, as the depth D approaches zero, the reflectances each approach the Lambertian limit.
The flat water layer model is more complex because of the Fresnel interface and the multiple underwater reflections. Combining the results of Eq. (101), the definition from Eq. (90), and the BRDF expression of Eq. (82) gives
To obtain approximate analytic expressions for the integrated reflectance terms, the additional simplification is made that the Fresnel transmittance within the angular integrands is unity for zenith angles below the critical angle, sin θ′<ri:
For the Lambertian limit, Eq. (89) becomes
As the depth of the water goes to zero, these reflectance simplify even further:
If one furthers considers the limit of a perfect Lambertian reflector, αi=1, and the Fresnel transmittances are set to one, the Lambertian limit is obtained. The Lambertian limit is also obtained for a black (αi=0) surface. However, Eq. (105) does not reduce to the Lambertian limit for grey bodies. In particular, if the reflectance is 50% (αi=½), then the albedo equals ri2/(1+ri2). For water, this ratio is closer to one-third than to one-half.
Δs
s
k
c
d
Although specific features of the invention are shown in some drawings and not others, this is for convenience only as some feature may be combined with any or all of the other features in accordance with the invention.
This application is a continuation in part of parent application Ser. No. 09/838,801 filed Apr. 20, 2001, now U.S. Pat. No. 7,433,806. Priority is claimed. The disclosure of the parent application is incorporated by reference herein. This application also claims priority of Provisional application Ser. No. 60/668,159 filed on Apr. 4, 2005
This invention was made with government support under contracts F19628-02-C-0078 and FA8718-04-C-0073 awarded by the Department of the Air Force. The government has certain rights in this invention.
| Number | Name | Date | Kind |
|---|---|---|---|
| 5315513 | Abreu et al. | May 1994 | A |
| 5884226 | Anderson et al. | Mar 1999 | A |
| 6584405 | Moncet | Jun 2003 | B1 |
| 7337065 | Adler-Golden et al. | Feb 2008 | B2 |
| 7433806 | Berk et al. | Oct 2008 | B2 |
| Number | Date | Country | |
|---|---|---|---|
| 20070027664 A1 | Feb 2007 | US |
| Number | Date | Country | |
|---|---|---|---|
| 60668159 | Apr 2005 | US |
| Number | Date | Country | |
|---|---|---|---|
| Parent | 09838801 | Apr 2001 | US |
| Child | 11398696 | US |
