Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Multi sky-view 3D aerosol distribution recovery

Open Access Open Access

Abstract

Aerosols affect climate, health and aviation. Currently, their retrieval assumes a plane-parallel atmosphere and solely vertical radiative transfer. We propose a principle to estimate the aerosol distribution as it really is: a three dimensional (3D) volume. The principle is a type of tomography. The process involves wide angle integral imaging of the sky on a very large scale. The imaging can use an array of cameras in visible light. We formulate an image formation model based on 3D radiative transfer. Model inversion is done using optimization methods, exploiting a closed-form gradient which we derive for the model-fit cost function. The tomography model is distinct, as the radiation source is unidirectional and uncontrolled, while off-axis scattering dominates the images.

© 2013 Optical Society of America

1. Introduction

Optical lightfield and integral imaging [13] sample the optical radiance distribution in location and direction. These are mainly used in small-scale setups. This paper deals with a far larger scale, to estimate scatterer 3D spatial distribution in the sky. The atmosphere allows light to pass through in multiple locations and directions. Light is affected by this medium. Therefore, in this paper we lay out a principle to recover this 3D distribution using measured and modeled light-fields (Fig. 1). 3D recovery of this medium has direct implications to various scientific communities that either rely on remotely-sensed imagery, study the atmosphere, or overcome the medium to see beyond. These include meteorology, atmospheric sciences, volcanology, and climatology. Aerosol retrieval is important for understanding climate evolution [4, 5] and monitoring air quality. Mapping aerosol density is significant to aviation safety, which needs real time assessment of conditions and visibility around flight paths.

 figure: Fig. 1

Fig. 1 Integral (lightfield) imaging through a volumetric distribution in the atmosphere, using ground-based cameras.

Download Full Size | PDF

In remote sensing [6], imaging through air is associated with atmospheric correction, based on aerosol retrieval. In this discipline, the atmosphere is often assumed to be plane-parallel, using one-dimensional vertical radiative transfer (RT). Consequently state-of-the-art aerosol retrieval is done in distinct large lateral blocks [7] with limited height resolution [5]. We, however, seek 3D recovery. Atmospheric correction is related to dehazing [8]. In this paper, however, the medium itself, at all relevant altitudes is the object of interest.

We rely on sky imaging from multiple directions and locations. Such projection is similar to tomography in other scientific domains [9, 10]. However, the situation here is distinct. Most tomography setups have a controlled and/or multidirectional radiation source [11, 12]. In contrast, our source is the uncontrolled, unidirectional Sun. Moreover, typical tomography relies on a linear model [13]: the pixel value is a linear combination of components along a line of sight (LOS), or a multiplicative combination (linearized by a logarithm). Linear tomography can detect gases [14], which absorb UV or IR radiation. However, aerosol attenuation of visible light is typically dominated by scattering, rather than absorption. Since the radiation source is single and fixed in space and time, we cannot rely on direct illumination for tomographic recovery of attenuation fields [15]. We may only use sunlight scattered into the LOS. The model is nonlinear, yet tractable. We model passive optical tomographic imaging of 3D atmospheric scatterer distributions in cloudless conditions. Then, we solve this tomography problem, to recover the distribution. Recovery is formulated as an optimization that minimizes a cost function. We derive the gradient of this cost function, to enable efficient optimization.

2. Theoretical background

Extinction: Sun rays (SR) irradiate a small volume that includes particles of a certain type. Each particle has an extinction cross section for interacting with the irradiance. An aerosol extinction cross section is denoted σaerosol. The aerosol density is n. Per unit volume, the extinction coefficient due to aerosols is βaerosol = σaerosoln. In addition, the atmosphere contains air molecules. The extinction coefficient due to the molecules is βair. The volume has infinitesimal length dl. Then, the relative portion of extinct SR irradiance is the unitless differential optical depth, = (βaerosol + βair)dl. The optical depth aggregates in extended propagation:

τ=dτ=(βaerosol+βair)dl=(σaerosoln+βair)dl=τair+σaerosolndl,
where τair = ∫ βairdl. Through an attenuating atmosphere, the transmittance exponentially decays with the optical depth:
t=exp(τ).

Scattering: Interaction of a single particle with the irradiance is by absorption and scattering. The weight of scattering (to all directions), relative to the total extinction is given by the unitless single scattering albedo of the particle. For an aerosol, the single scattering albedo is denoted ϖaerosol. The scattering coefficient due to aerosols in the volume is

αaerosol=ϖaerosolβaerosol=ϖaerosolσaerosoln.
The angular distribution of scattering by an aerosol is determined by a phase function Paerosol, which is normalized: its integral over all solid angles is unit. Part of the light scatters towards a camera’s LOS, as illustrated in Fig. 1. The angle between the SR and LOS is the scattering angle Φscatter. The angular scattering coefficient due to aerosols is
α˜aerosol(Φscatter)=αaerosolPaerosol(Φscatter)=ϖaerosolσaerosolnPaerosol(Φscatter).
The phase function is often approximated by a parametric expression. Specifically, the Henyey-Greenstein [16] function is governed by an anisotropy parameter g:
Pgaerosol(Φscatter)38π(1g2)[1+(cosΦscatter)2](2+g2)(1+g22gcosΦscatter)32

Scattering by air molecules follows the Rayleigh model. The phase function is

Pair(Φscatter)316π(1+cos2Φscatter)
and the single scattering albedo in visible light is ϖair=1. Air molecular density falls off approximately exponentially with altitude h, with a characteristic [17] falloff height Hair = 8 km. Thus, the coefficients for extinction and scattering by air molecules can be modeled by [17],
αair(h,λ)=βair(h,λ)1.09×103λ4exp(h/Hair),
where λ is the wavelength, in microns.

Inverse transform sampling: A tool for RT modeling is the Monte-Carlo (MC) method. MC considers scattering and extinction as random phenomena, that are sampled from their probability distributions. Inverse transform sampling [18] is used for sampling random numbers that comply with a specified probability density function. Let u be a random number drawn from a uniform distribution in the interval [0, 1]. The number u can be transformed into a random variable X, whose cumulative distribution function (CDF) is F(X). The transform is X = F−1(u), where F−1 denotes the inverse function of F.

For example, consider a photon propagating in the atmosphere. The photon has high probability of propagating as long as t is high, and the probability diminishes as t → 0. Thus (Eq. (2)) can be viewed as a probability density function, whose CDF is

F(τ)=0τexp(τ)dτ=1exp(τ).
Thus a photon propagates to a random optical depth
τrandom=F1(u)=ln(1u).
If the atmosphere is uniform, then from Eq. (1) the photon reaches a random distance
lrandom=(βaerosol+βair)1τrandom=(βaerosol+βair)1ln(1u),
before interacting with a particle. If the number of photons launched is very high, their average number falls off in consistency with Eqs. (1) and (2). Analogous transforms generate scattering at random angles, according to the phase function.

3. 3D recovery approach and its assessment

A set of ground-based wide-angle cameras looks upwards (Fig. 1), performing integral imaging [19] of the sky. Per viewpoint, the view azimuth and elevation angles (relative to the zenith and north) are encapsulated in vector Θ. Any camera c is at a distinct fixed location, capturing raw image ic(Θ). This paper reports a first step: formulating a principle for estimating 3D atmospheric aerosols distributions, using such data. To theoretically study the feasibility and cross validate it, we perform RT using two different forward models. A forward model takes as input a 3D aerosol distribution, and outputs images as if captured at various viewpoints. One model uses MC (Sec. 6): it is stochastic, slow but naturally expresses multiple-scattering. Hence, MC rigorously simulates realistic scenes of arbitrary complexity. The second model uses the single-scattering approximation (Sec. 4). It is less accurate than MC, but solves RT in an analytic, closed form. This form enables simple inversion of the model, to estimate the underlying aerosol distribution (Sec. 7).

4. Single-scattering forward model

As illustrated in Figs. 1 and 2, the sky is discretized into a grid of Nvoxels rectangular cuboid volume elements (voxels), indexed by k, m or q. In a single-scattering model, any SR changes direction only once on the way to a camera. This model is valid in atmospheres that are not very dense (inside fog and clouds this model does not apply). Here, RT has three steps (Fig. 1): ( i) Attenuation of a SR propagating from the top of the atmosphere (TOA) to voxel k. ( ii) Light scattering at voxel k, towards a camera. ( iii) Light attenuation on the LOS from voxel k to the camera. Steps i, iii involve the optical depth along the light path to and from voxel k (analyzed in Sec. 4.1). Step ii involves the scattering coefficient at voxel k (Sec. 4.2).

 figure: Fig. 2

Fig. 2 Simulated 3D aerosol distributions: [a] Haze blobs and [b] Haze Front. Here, an aerosol density unit is 106 particles/m3. The atmopsheric domain has area of 50 × 50km2, extending from the ground up to altitude of 10km. The domain is divided into a rectilinear grid having Nvoxels voxels.

Download Full Size | PDF

4.1. Optical depth

Denote a SR line segment between the TOA and the center of voxel k by [SR, k] (See Fig. 1). This SR intersects voxel q. The intersection line segment has geometric length lSR(qz, ΦSR, k), where Δz is a voxel’s vertical geometric thickness, and ΦSR is the Sun angle. Define a Nvoxels × Nvoxels matrix DSun→voxel, whose element (k, q) represents lSR(q):

DSunvoxel(k,q)={lSR(q)ifq[SR,k]0otherwise
Matrix DSun→voxel is sparse. Analogously, for camera c, denote by [LOSc, k] a LOS bounded between the camera and the center of voxel k (see Fig. 1). The LOS zenith angle is ΦcLOS. Suppose this LOS intersects voxel q. The geometric length of this intersection line segment is lLOSc(mz, ΦcLOS, k). As in Eq. (11), define a Nvoxels × Nvoxels sparse matrix whose element (k, m) is
Dcvoxelcam(k,m)={lLOSc(m)ifm[LOSc,k]0otherwise.

We focus on situations where, in addition to air molecules, there is a single type of aerosol over a scene. Hence, the three-element vector [σaerosol, ϖaerosol, g] is uniform across the scene, but the density distribution n(k) is variable. As a numerical approximation, assume that within any voxel, the molecular parameters {βair(k), αair(k)} and the aerosol density n(k) are constants, e.g., corresponding to the value at each voxel center. Following Eq. (1), the optical depths between the TOA and voxel k, and from voxel k to camera c are

τSR(k)=q[SR,k]lSR(q)[βair(q)+σaerosoln(q)]τLOSc(k)=m[SR,k]lLOSc(m)[βair(m)+σaerosoln(m)].
Let n, τSR and βair be column stack vector representations of n(k), τSR(k) and βair(k), respectively. Then, we can write Eq. (13) using matrix notation:
τSR=DSunvoxelβair+σaerosolDSunvoxelnτLOSc=Dcvoxelcamβair+σaerosolDcvoxelcamn.

The SR and LOS paths joining at each voxel form a single path from the TOA to camera c. Hence, in matrix notation, Eq. (1) yields the total optical depths corresponding to LOSs (of camera c) and SRs that cross all voxels:

τc=τcair+σaerosolDcn,
where Dc=DSunvoxel+Dcvoxelcam and τcair=Dcβair.

4.2. Scattering

Per camera c and voxel k, the lines [SR, k] and [LOSc, k] intersect at a fixed angle Φcscatter(k). Column-stacking all angles yields the vector representation of all scattering angles in the domain, Φcscatter per camera c. Using Eqs. (4) and (6), the angular scattering coefficients across the domain are expressed in vector form by

α˜cair=βairPair(Φcscatter),α˜caerosol=ϖaerosolσaerosolnPgaerosol(Φcscatter).
Here the operator ⊙ denotes the Hadamard (element-wise) product.

4.3. Image capture

Compounding the attenuation of irradiance along both [SR, k] and [LOSc, k], and scattering by voxel k towards the camera (Eqs. (2), (4) and (15)) the voxel contributes radiant power

pc(k)=LTOA[α˜cair(k)+α˜caerosol(k)]exp[τc(k)],
where LTOA is the radiance at the TOA. A column-stack vector of all voxel contributions is
pc=LTOA(α˜cair+α˜caerosol)exp(τc).

A camera sensor comprises Npix pixels. Each pixel collects light from a narrow cone in the atmosphere. The cone contains or intersects some voxels, while being oblivious to all the rest. Overall, light power captured at the pixel is a weighted sum of the power pc(k) radiating from all voxels. This sum is formulated by a matrix operation Πc over pc:

ic=Πcpc,
where ic is the image, column-stacked to a vector Npix long. Combining Eqs. (15), (16), (18) and (19), the captured image is thus
ic=LTOAΠc{[α˜cair+ϖaerosolσaerosolPgaerosol(Φcscatter)n]exp[(τcair+σaerosolDcn)]}.

5. Examples

Figure 3 shows examples of photographs rendered using the single-scattering model. To create such photographs we simulated several scenarios, whose details are as follows.

 figure: Fig. 3

Fig. 3 Different viewpoints of a sky that includes Haze Blobs. Images are rendered using the single-scattering model. A yellow dot marks the Sun location. Dashed white circles mark zenith angles.

Download Full Size | PDF

Geometry: The sun is at zenith angle ΦSR = 45°. Its LTOA is obtained from [20, 21], with respective red-green-blue ratios of LRTOA:LGTOA:LBTOA=255:236:224.

The atmospheric domain is described in Fig. 2. We tested domain division to various grids, from 10 × 10 × 20 voxels, up to 50 × 50 × 100. The corresponding voxels had width, length and thickness of 5km × 5km × 0.5km and 1km × 1km × 0.1km. In all cases, the vertical resolution was higher than horizontal, via use of voxels that are vertically thin but laterally wide. This expresses the natural higher variability of air and aerosols in the vertical dimension. Each ground-based camera has a hemispherical field of view, and low resolution, 128px × 128px, as cloudless (yet hazy) sky images are rather smooth.

Aerosol: We used particle-type 6 from the aerosol list in [22]. This is a spherical non-absorbing sea-salt/organic particle, whose anisotropy parameter per color channel is [gR, gG, gB] = [0.763, 0.775, 0.786]. Its corresponding extinction cross sections are [σRaerosol,σGaerosol,σBaerosol]=[16.5,16.2,15.9]μm2. At all channels, ϖaerosol = 1. We simulated spatial distributions using a product of two functions. To express the general trend of reduced density with altitude h(k) of voxel k, we follow [17] and set the first function as

f1[h(k)]=nsealevelexp[h(k)/Haerosol],
where Haerosol = 2 km. To express a clustered nature of aerosol distributions (as clouds), we define blobs in the form of 3D ellipsoids. There may be multiple ellipsoids suspended. Then
f2(k)={1ifkanyblobcluster0otherwise.
The true aerosol distribution is then ntrue = 𝒮 {f1f2}, where 𝒮 is 3D spatial smoothing, obtained by narrow 3D Gaussian filtering. It expresses non-sharp boundaries of typical distributions. The Haze Blobs scene in Fig. 2(a) uses two ellipsoids: one is 32km wide, 2.8km thick and centered at altitude 2.5km; the other is 24km wide, 2.1km thick and centered at altitude 5km. Horizontal widths are much larger than the vertical thickness, in consistency with atmospheric scales. Figure 3 shows photographs of the Haze Blobs scene, as simulated by Eq. (20). For clarity of display, the brightness of the displayed pictures in this paper underwent the same standard contrast enhancement (stretching), including gamma-correction. The recovery algorithm, of course, uses the raw (not brightness enhanced) images.

The Haze Front scene in Fig. 2(b) uses an ellipsoid degenerated to an elliptic cylinder that partly enters the analyzed volume. It is 32km long (only part of which enters the volume), 4km thick and centered at altitude 5km. The front stretches across the domain width.

6. Multiple-scattering forward model by Monte Carlo

The MC simulations rely on the same grid as Sec. 4. We use backward MC: from a detector pixel, photons are ‘launched’ via the camera lens; then the photons are traced back through the atmosphere. Photons that happen to be back-traced to the Sun are counted, per pixel. In this model, RT takes the following steps, per pixel (Fig. 4a):

  1. Launch a photon-packet from camera c to direction Θ. This is the initial ray, denoted 0. The packet has intensity I0.

    Per iteration i:

  2. Find the distance on ray i to which the photon-packet propagates. To do this, Eq. (9) yields τrandom. Then, Eq. (10) generalizes to a non-uniform atmosphere: using Eq. (1), numerically seek lrandom s.t.
    0lrandom(σaerosoln+βair)dl=τrandom
    Distance lrandom along i yields 3D position Xi. If Xi is outside the atmospheric grid, the packet is terminated. If, in addition, i || SR, the packet is counted as contributing to the image pixel.
  3. The type of particle (air molecule or aerosol) that the photon-packet interacts with at point Xi is determined randomly. The relative probabilities of this random step are set by the respective extinction coefficients (βair vs. βaerosol) at the voxel containing Xi.
  4. If the particle is an aerosol, the photon-packet intensity is attenuated to Ii+1 = ϖaerosolIi. If Ii+1 is lower than a threshold, the packet is stochastically terminated, following [23].
  5. The photon-packet is scattered to a new direction, determined by inverse transform sampling, according to the phase function of the particle (Eqs. 5 or 6). This yields a new ray, denoted i+1, emanating from Xi. Local estimation [24] derives the amount of light back-traced to the sun at each scattering. The next iteration of propagation (denoted ii above) proceeds.

 figure: Fig. 4

Fig. 4 [a] Multiple-scattering forward model by MC. [b] Photograph of the same sky and view point as in Fig. 3(a), rendered using MC.

Download Full Size | PDF

The quality of MC increases with the number of photons launched. Figure 4(b) shows a photograph of the Haze Blobs scene, derived by the MC process, using 10000 photons per pixel, from the same viewpoint as Fig. 3(a). Similarly, Fig. 5(a) corresponds to the same viewpoint as Fig. 3(d). Figure 6 shows the contribution by successive scattering orders. They are derived from the MC image. Photons contributing to any pixel are accumulated in steps ii and v above. The contributing photons that stem from just a single scattering event yield Fig. 6(a). Similarly, contributing photons that had experienced exactly two or three scattering events yield respectively Fig. 6(b) and Fig. 6(c). First-order scatter yields most of the energy in the MC image. Radiance by higher order scattering is mostly evident at the horizon. A horizontal LOS passes through a long and dense part of the atmosphere, while a vertical LOS cuts short through the dense, lower part. Hence, toward the horizon the probability of high order scattering increases.

 figure: Fig. 5

Fig. 5 Images simulated by MC, from the same viewpoint, but different aerosol characteristics. (a) Haze blobs, characterized in Sec. 5. (b) Partly absorbing aerosol. (c) Aerosol having an isotropic phase function. (d) High aerosol density. Details of the scenarios used in creating (b) and (c) are given in Sec. 8.

Download Full Size | PDF

 figure: Fig. 6

Fig. 6 Separation of the image shown in Fig. 4(b) to contributions by successive orders of scattering: (a) first, (b) second, (c) third and (d) forth scattering order.

Download Full Size | PDF

Figure 7 plots the image intensity profile along a straight line from left to right, via the image center (zenith). The plot compares photographs simulated using the single-scattering and the MC forward models. Figure 7(a) is extracted from Fig. 3(d) and Fig. 5(a) of the Haze Blobs scene. This plot demonstrates consistency in this scene, for most pixels. Figure 7(b) plots the same profiles, where the aerosol density was increased tenfold (see Sec. 8 and the corresponding MC photograph, Fig. 5(d)). In Fig. 7(b) the plots differ more than in Fig. 7(a), due to high order scattering.

 figure: Fig. 7

Fig. 7 Cross-sections along the X-axis of photographs rendered by the single-scattering (dotted) and MC (solid line) forward models. (a) cross sections of the Haze Blobs scene (described in Sec. 5). The difference between the models is much more pronounced (b) when the scatterers are ten times denser (see Sec. 8).

Download Full Size | PDF

7. Inverse problem

When the type of aerosol above an area is known [22], estimation of the density distribution n is needed. The data are Nviews measured photographs {icmeasured}c=1Nviews. Recovery is phrased as optimization of a cost function that fits the image-formation model to the data. Let 𝒞 be the set of all distributions complying with some constraints. Particularly, n is non-negative, and its spatial support is bounded between the ground and the TOA. The optimization problem is

n^=argminn𝒞E(n)
where
E(n)=c=1Nviews[icmeasuredic(n)]22+ηΨ(n).
Here masks the area around the sun: in real-world images it is indeed blocked. We mask the horizon, by using a camera whose field of view is a little narrower than hemispherical. This reduces the influence of LOSs that are more strongly affected by high-order scattering (Fig. 6(b)). Consequently, the fit of the imaging model focuses on the bulk image region, where the single scattering model is more valid. In Eq. (25), Ψ(n) is a regularization term that expresses some prior knowledge on the distribution, while η is the regularization weight.

Since aerosol distributions are usually fuzzy, useful regularization is by a smoothness term, which penalizes for energy in the second order spatial derivatives (3D Laplacian), 2n22. Regularization is not required when data is sufficient and reliable. Voxels at high altitude have good coverage by many cameras at multiple directions (Fig. 1), while low altitude voxels are mainly observed by local cameras. Hence, regularization can be weakened with altitude. We accomplish this weakening using a weight w(k) = exp[−h(k)/Hsmooth]. Overall we use

Ψ(n)=Wn22
where is a matrix representation of the 3D Laplacian operator, and matrix W is diagonal, whose elements are w(k).

We solve Eq. (24) using standard optimization tools. The gradient of E with respect to n is

nE=2c=1Nviews[Jic(n)][icmeasuredic(n)]+2ηWWn.
Here the matrix Jic(n) is the Jacobian of the vector ic with respect to n. Element (Θ, k) of the Jacobian differentiates the intensity in pixel Θ (in viewpoint c) with variation of the aerosol density at voxel k, i.e., ∂ic(Θ)/∂n(k).

We now detail the derivation of the Jacobian Jic(n). First, we provide some results relating to differentiation. Let a(n), u(n) be vector functions: each outputs a vector of length r. Let C be a r × Nvoxels matrix, where Nvoxels is the length of n. Then,

C(au)n=[an𝔻{u}+un𝔻{a}]C.
Here 𝔻{v} is a conversion of a general vector v into a diagonal matrix, whose main diagonal elements correspond to the elements of v. Now, let a(n) = exp(−Cn), where the exponent is element-wise (not raising an operator to some power). Then
an=C𝔻{exp(Cn)}.
Using Eq. (29),
exp[(τcair+σaerosolDcn)]n=σaerosolDc𝔻{exp[(τcair+σaerosolDcn)]},
where τcair is independent of n. Using Eq. (28),
[ϖaerosolσaerosolPgaerosol(Φcscatter)n]n=ϖaerosolσaerosol𝔻{Pgaerosol(Φcscatter)}.
Based on Eqs. (28), (30) and (31), we derive the Jacobian of Eq. (20) in close-form,
Jic(n)=LTOAσaerosol(AB)𝔻{exp[(τcair+σaerosolDcn)]}Πc,
where
A=ϖaerosol𝔻{Pgaerosol(Φcscatter)},B=Dc𝔻{[α˜cair+ϖaerosolσaerosolPgaerosol(Φcscatter)n]}

8. Recovery simulations

To demonstrate recovery, we performed simulations, and tested effects of the following variations: {1} Multiple or only single scatter in the forward model; {2} Spatial distribution; {3} Density of viewpoints; {4} Aerosol characteristics; {5} Density scale.

In all cases, images were created as described in Secs. 5 and 6. These images served as input to the reconstruction algorithm. The optimization used an L-BFGS-B solver [25] on a computer cluster. Each computer core was dedicated to rendering a modeled image. The optimization was initialized by nnn = 0. Satisfactory convergence occurred in several hundred iterations. Depending on the resolution it took between minutes to a couple of hours to complete, in total. The total estimation error can be quantified by the aerosol mass that is over and under-estimated in all voxels, relative to the total aerosol mass in the scene. Using the 1 norm, the total mass relative difference is δmass = (||||1 − ||ntrue||1)/||ntrue||1. To also sense local errors, we use ε = ||ntrue||1/||ntrue||1. These results are listed in Table 1.

Tables Icon

Table 1. Relative errors in various simulations. Here 𝒪 denotes order of magnitude.

Forward models and distributions: Figure 8 shows estimation results, corresponding to the original distributions of Fig. 2. Here, the atmosphere domain was defined over a 20 × 20 × 40 grid. Cameras were placed ∼ 7km apart on a 6 × 6 grid. In Figs. 8(a) and 8(b), the single-scattering model was used both in rendering of the raw images, and in the estimation algorithm. In Figs. 8(c) and 8(d), the MC model (multiple scattering) rendered the raw images, over which our recovery (Sec. 7) was tested. Overall, the distributions, the density order of magnitude and values are consistent. Errors in Figs. 8(a) and 8(b) stem from random noise, which had been induced in the raw images. Errors are larger when MC is used to render the images, since high-order scattering is not accounted for in the algorithm of Sec. 7.

 figure: Fig. 8

Fig. 8 Recoveries of scenes shown in Fig. 2. Color represents aerosol density, in units of 106 particles/m3. In (a) and (b) images simulated by single-scattering are used as input for the recovery. In (c) and (d), MC simulated images are used as input for the recovery.

Download Full Size | PDF

Density of viewpoints: Keeping the domain and scene (Haze blobs) fixed, we repeated the reconstruction many times, each using a different subset of the above mentioned 36 simulated cameras. Each test randomly selected viewpoints, for various fixed values of Nviews < 36. All rendered images used MC. Figure 9 plots the performance. Diminishing returns are obtained beyond Nviews ≈ 20, which corresponds to ≈ 11km separation between neighboring viewpoints. These tests point to fundamental questions about spatial-angular sampling: how spatially dense should cameras be, and how dense should a camera’s field of view be sampled (in pixels)? This likely depends on the finesse of atmospheric structure sought, and the level of optical diffusion achieved by various aerosols. The novel tomography domain introduced here thus raises new theoretical questions that will require extensive further research.

 figure: Fig. 9

Fig. 9 The relative error ε measure decreases with the number of cameras. Bars represent the standard deviation of ε, over our random tests.

Download Full Size | PDF

Aerosol characteristics: We used the Haze blobs scene, but varied the aerosol characteristics. Compared to the aerosol described in Sec. 5, we tested two variations. In one variation, the phase function was isotropic, i.e., [gR, gG, gB] = [0, 0, 0]. In the other variation, the aerosol is partly absorbing: its single scattering albedo per color channel is [ϖRaerosol,ϖGaerosol,ϖBaerosol]=[0.69,0.74,0.78]. All rendered images used MC (See Figs. 5(a)–5(c)). Recovery results using Nviews = 36 are shown in Figs. 10(a) and 10(b) and in Table 1. The results indicate that if the aerosol’s phase function has smaller anisotropy, recovery of the aerosol’s density is easier. We hypothesize that tomography is better served by scattering isotropy, since then significant radiance from any voxel reaches cameras in a wide range of directions. This may make back-projection (recovery from cameras to voxels) more accurate. On the other hand, if the phase function is too strongly peaked around one direction, most cameras would not receive much scattered radiance from a voxel, undermining tomography. A very broad research is needed to verify and quantify such dependencies.

 figure: Fig. 10

Fig. 10 Recovered distributions when the aerosol is either party absorbing (a), has isotropic phase function (b) or has high density (c).

Download Full Size | PDF

Density scale: We used the aerosol Haze blobs distribution and characteristics as described in Sec. 5, but increased the aerosol density tenfold everywhere. This significantly increases the effects of high order scattering, as shown in Fig. 7(b) and Fig. 5(d). Consequently, ε increases (Table 1), when using Nviews = 36. Still, Fig. 10(c) shows that the spatial distribution outlines are recovered similarly to the other tests.

9. Discussion

We describe a way to sense the 3D atmosphere. It can help remote sensing depart from common one-dimensional aerosol distribution models, in favor of detailed 3D RT analysis. The paper describes both a novel data acquisition system concept for this recovery task (ground-based camera network), and a dedicated algorithm for reconstructing aerosol distributions. Using more prior knowledge on the nature of distributions will improve the recovery.

The estimation algorithm uses a single scattering model. This approximation is valid when the optical depth is small. However, as MC simulations show, it is important to generalize the recovery algorithm so that high-order scattering is included in it. Such generalization would be computationally complex: in our tests, the MC forward model was slower by 2–3 orders of magnitude than a single-scattering analytic model. Thus, inverting a 3D multiple-scattering model requires research into numerical schemes for increasing the solution efficiency of both forward and inverse problems. The remote sensing community has devised several approaches to better approximate multiple-scattering problems, mainly in one-dimensional models (e.g., successive orders of scattering and the discrete ordinates method). We may tap into some of these approaches. Deploying an experimental system will be needed: camera specifications [26] should be set to optimally recover a wide range of aerosol distributions.

Acknowledgments

We are grateful to David J. Diner for fruitful discussions and his support at all levels. We thank the reviewers for useful comments. Yoav Schechner is a Landau Fellow - supported by the Taub Foundation. This work is supported by the US-Israel Binational Science Foundation (BSF grant 2012202) and conducted in the Ollendorff Minerva Center. Minerva is funded through the BMBF. Part of this work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA and with support by NASA’s Earth Science Technology Office Advanced Information Systems Technology Program.

References and links

1. A. Stern and B. Javidi, “Three-dimensional image sensing, visualization, and processing using integral imaging,” Proceedings of the IEEE 94.3 (2006).

2. J. Kim, D. Lanman, Y. Mukaigawa, and R. Raskar, “Descattering transmission via angular filtering,” in Proc. ECCV’ 10, (Springer-Verlag, Berlin, Heidelberg, 2010), pp. 86–99.

3. M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz, “Light field microscopy,” ACM Transactions on Graphics (TOG) 25, 924–934 (2006). [CrossRef]  

4. U. Dayan, B. Ziv, T. Shoob, and Y. Enzel, “Suspended dust over southeastern Mediterranean and its relation to atmospheric circulations,” International Journal of Climatology 924, 915–924 (2008). [CrossRef]  

5. O. V. Kalashnikova, M. J. Garay, A. B. Davis, D. J. Diner, and J. V. Martonchik, “Sensitivity of multi-angle photo-polarimetry to vertical layering and mixing of absorbing aerosols: Quantifying measurement uncertainties,” Journal of Quantitative Spectroscopy and Radiative Transfer 112, 2149–2163 (2011). [CrossRef]  

6. M. I. Mishchenko and I. V. Geogdzhayev, “Satellite remote sensing reveals regional tropospheric aerosol trends,” Opt. Express 15, 7423–38 (2007). [CrossRef]   [PubMed]  

7. J. Martonchik, D. D. J. Diner, J. M. David, R. Kahn, T. P. Ackerman, M. M. Verstraete, B. Pinty, and H. R. Gordon, “Techniques for the Retrieval of Aerosol Properties over Land and Ocean Using Multi-angle Imaging,” IEEE Trans. on Geoscience and Remote Sensing 36, 1212–1227 (1998). [CrossRef]  

8. E. Namer, S. Shwartz, and Y. Y. Schechner, “Skyless polarimetric calibration and visibility enhancement,” Opt. Express 17, 472–93 (2009). [CrossRef]   [PubMed]  

9. A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9, 272–86 (2001). [CrossRef]   [PubMed]  

10. I. Ihrke, K. N. Kutulakos, H. P. Lensch, M. Magnor, and W. Heidrich, “State of the art in transparent and specular object reconstruction,” in EUROGRAPHICS 2008 STAR–STATE OF THE ART REPORT , (2008).

11. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” Sig.Proc. Magazine 18, 57–75 (2001). [CrossRef]  

12. H. Messer, A. Zinevich, and P. Alpert, “Environmental sensor networks using existing wireless communication systems for rainfall and wind velocity measurements,” IEEE Instrumentation & Measurement Magazine 15, 32–38 (2012). [CrossRef]  

13. J. Gregson, M. Krimerman, M. B. Hullin, and W. Heidrich, “Stochastic tomography and its applications in 3D imaging of mixing fluids,” ACM Trans. Graph. 31, 52:1–52 (2012). [CrossRef]  

14. B. R. Cosofret, D. Konno, A. Faghfouri, H. S. Kindle, C. M. Gittins, M. L. Finson, T. E. Janov, M. J. Levreault, R. K. Miyashiro, and W. J. Marinelli, “Imaging sensor constellation for tomographic chemical cloud mapping,” Appl. Opt. 48, 1837–52 (2009). [CrossRef]   [PubMed]  

15. J. A. Aviles, “The Development and Validation of a First Generation X-Ray Scatter Computed Tomography Algorithm for the Reconstruction of Electron Density Breast Images Using Monte Carlo Simulation,” Ph.D. thesis (2011).

16. W. M. Cornette and J. G. Shanks, “Physically reasonable analytic expression for the single-scattering phase function,” Appl. Opt. 31, 3152–3160 (1992). [CrossRef]   [PubMed]  

17. L. Levi, Applied Optics (John Wiley & Sons, Inc., 1980).

18. L. Devroye, “Sample-based non-uniform random variate generation,” in Proceedings of the 18th conference on Winter simulation , (ACM, 1986), pp. 260–265. [CrossRef]  

19. S.-H. Hong, J.-S. Jang, and B. Javidi, “Three-dimensional volumetric object reconstruction using computational integral imaging,” Opt. Express 12, 483–491 (2004). [CrossRef]   [PubMed]  

20. M. Charity, “Blackbody color datafile,” (2001), http://www.vendian.org/mncharity/dir3/blackbody/UnstableURLs/bbr\_color.html.

21. Wikipedia, “Sunlight — Wikipedia, The Free Encyclopedia,” (2012), http://en.wikipedia.org/w/index.php?title=Sunlight\&oldid=502554571.

22. J. V. Martonchik, R. A. Kahn, and D. J. Diner, “Retrieval of aerosol properties over land using MISR observations,” in Satellite Aerosol Remote Sensing over Land,, A. A. Kokhanovsky and G. Leeuw, eds. (Springer BerlinHeidelberg, 2009), pp. 267–293. [CrossRef]  

23. H. Iwabuchi, “Efficient Monte Carlo Methods for Radiative Transfer Modeling,” Journal of the Atmospheric Sciences 63, 2324–2339 (2006). [CrossRef]  

24. A. Marshak and A. Davis, 3D Radiative Transfer in Cloudy Atmospheres, Physics of Earth and Space Environments (Springer, 2005). [CrossRef]  

25. C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization,” ACM Trans. Math. Softw. 23, 550–560 (1997). [CrossRef]  

26. N. J. Pust, A. R. Dahlberg, M. J. Thomas, and J. a. Shaw, “Comparison of full-sky polarization and radiance observations to radiative transfer simulations which employ AERONET products,” Opt. Express 19, 18602–18613 (2011). [CrossRef]   [PubMed]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)

Fig. 1
Fig. 1 Integral (lightfield) imaging through a volumetric distribution in the atmosphere, using ground-based cameras.
Fig. 2
Fig. 2 Simulated 3D aerosol distributions: [a] Haze blobs and [b] Haze Front. Here, an aerosol density unit is 106 particles/m3. The atmopsheric domain has area of 50 × 50km2, extending from the ground up to altitude of 10km. The domain is divided into a rectilinear grid having Nvoxels voxels.
Fig. 3
Fig. 3 Different viewpoints of a sky that includes Haze Blobs. Images are rendered using the single-scattering model. A yellow dot marks the Sun location. Dashed white circles mark zenith angles.
Fig. 4
Fig. 4 [a] Multiple-scattering forward model by MC. [b] Photograph of the same sky and view point as in Fig. 3(a), rendered using MC.
Fig. 5
Fig. 5 Images simulated by MC, from the same viewpoint, but different aerosol characteristics. (a) Haze blobs, characterized in Sec. 5. (b) Partly absorbing aerosol. (c) Aerosol having an isotropic phase function. (d) High aerosol density. Details of the scenarios used in creating (b) and (c) are given in Sec. 8.
Fig. 6
Fig. 6 Separation of the image shown in Fig. 4(b) to contributions by successive orders of scattering: (a) first, (b) second, (c) third and (d) forth scattering order.
Fig. 7
Fig. 7 Cross-sections along the X-axis of photographs rendered by the single-scattering (dotted) and MC (solid line) forward models. (a) cross sections of the Haze Blobs scene (described in Sec. 5). The difference between the models is much more pronounced (b) when the scatterers are ten times denser (see Sec. 8).
Fig. 8
Fig. 8 Recoveries of scenes shown in Fig. 2. Color represents aerosol density, in units of 106 particles/m3. In (a) and (b) images simulated by single-scattering are used as input for the recovery. In (c) and (d), MC simulated images are used as input for the recovery.
Fig. 9
Fig. 9 The relative error ε measure decreases with the number of cameras. Bars represent the standard deviation of ε, over our random tests.
Fig. 10
Fig. 10 Recovered distributions when the aerosol is either party absorbing (a), has isotropic phase function (b) or has high density (c).

Tables (1)

Tables Icon

Table 1 Relative errors in various simulations. Here 𝒪 denotes order of magnitude.

Equations (33)

Equations on this page are rendered with MathJax. Learn more.

τ = d τ = ( β aerosol + β air ) d l = ( σ aerosol n + β air ) d l = τ air + σ aerosol n d l ,
t = exp ( τ ) .
α aerosol = ϖ aerosol β aerosol = ϖ aerosol σ aerosol n .
α ˜ aerosol ( Φ scatter ) = α aerosol P aerosol ( Φ scatter ) = ϖ aerosol σ aerosol n P aerosol ( Φ scatter ) .
P g aerosol ( Φ scatter ) 3 8 π ( 1 g 2 ) [ 1 + ( cos Φ scatter ) 2 ] ( 2 + g 2 ) ( 1 + g 2 2 g cos Φ scatter ) 3 2
P air ( Φ scatter ) 3 16 π ( 1 + cos 2 Φ scatter )
α air ( h , λ ) = β air ( h , λ ) 1.09 × 10 3 λ 4 exp ( h / H air ) ,
F ( τ ) = 0 τ exp ( τ ) d τ = 1 exp ( τ ) .
τ random = F 1 ( u ) = ln ( 1 u ) .
l random = ( β aerosol + β air ) 1 τ random = ( β aerosol + β air ) 1 ln ( 1 u ) ,
D Sun voxel ( k , q ) = { l SR ( q ) if q [ SR , k ] 0 otherwise
D c voxel cam ( k , m ) = { l LOS c ( m ) if m [ LOS c , k ] 0 otherwise .
τ SR ( k ) = q [ SR , k ] l SR ( q ) [ β air ( q ) + σ aerosol n ( q ) ] τ LOS c ( k ) = m [ SR , k ] l LOS c ( m ) [ β air ( m ) + σ aerosol n ( m ) ] .
τ SR = D Sun voxel β air + σ aerosol D Sun voxel n τ LOS c = D c voxel cam β air + σ aerosol D c voxel cam n .
τ c = τ c air + σ aerosol D c n ,
α ˜ c air = β air P air ( Φ c scatter ) , α ˜ c aerosol = ϖ aerosol σ aerosol n P g aerosol ( Φ c scatter ) .
p c ( k ) = L TOA [ α ˜ c air ( k ) + α ˜ c aerosol ( k ) ] exp [ τ c ( k ) ] ,
p c = L TOA ( α ˜ c air + α ˜ c aerosol ) exp ( τ c ) .
i c = Π c p c ,
i c = L TOA Π c { [ α ˜ c air + ϖ aerosol σ aerosol P g aerosol ( Φ c scatter ) n ] exp [ ( τ c air + σ aerosol D c n ) ] } .
f 1 [ h ( k ) ] = n sealevel exp [ h ( k ) / H aerosol ] ,
f 2 ( k ) = { 1 if k any blob cluster 0 otherwise .
0 l random ( σ aerosol n + β air ) d l = τ random
n ^ = arg min n 𝒞 E ( n )
E ( n ) = c = 1 N views [ i c measured i c ( n ) ] 2 2 + η Ψ ( n ) .
Ψ ( n ) = W n 2 2
n E = 2 c = 1 N views [ J i c ( n ) ] [ i c measured i c ( n ) ] + 2 η W W n .
C ( a u ) n = [ a n 𝔻 { u } + u n 𝔻 { a } ] C .
a n = C 𝔻 { exp ( Cn ) } .
exp [ ( τ c air + σ aerosol D c n ) ] n = σ aerosol D c 𝔻 { exp [ ( τ c air + σ aerosol D c n ) ] } ,
[ ϖ aerosol σ aerosol P g aerosol ( Φ c scatter ) n ] n = ϖ aerosol σ aerosol 𝔻 { P g aerosol ( Φ c scatter ) } .
J i c ( n ) = L TOA σ aerosol ( A B ) 𝔻 { exp [ ( τ c air + σ aerosol D c n ) ] } Π c ,
A = ϖ aerosol 𝔻 { P g aerosol ( Φ c scatter ) } , B = D c 𝔻 { [ α ˜ c air + ϖ aerosol σ aerosol P g aerosol ( Φ c scatter ) n ] }
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.