Settings for Spaceborne 3-D Scattering Tomography of Liquid-Phase Clouds by the CloudCT Mission

We introduce a comprehensive method for space-borne 3-D volumetric scattering-tomography of cloud microphysics, developed for the CloudCT mission. The retrieved microphysical properties are the liquid-water-content (LWC) and effective droplet radius within a cloud. We include a model for a perspective polarization imager and an assumption of 3-D variation of the effective radius. Elements of our work include computed tomography initialization by a parametric horizontally uniform microphysical model. This results in smaller errors than the prior art. The mean absolute errors of the retrieved LWC and effective radius are reduced from 62% and 28% to 40% and 9%, respectively. The parameters of this initialization are determined by a grid search of a cost function. Furthermore, we add viewpoints in the cloudbow region, to better sample the polarized scattering phase function. The suggested advances are evaluated by retrieval of a set of clouds generated by large-eddy simulations.

statistics. The statistics will be used to study cloud trends in 37 changing environmental conditions. Based on these expected 38 studies, improved parameterizations of climate models will be 39 developed. 40 In CloudCT a formation of ten satellites will capture simul-41 taneously multiview images of cloud fields. These will be 42 analyzed by scattering tomography, based on a 3-D radiative-43 transfer (RT) solver that includes multiple scattering. Images 44 will be captured in daytime, exploiting solar radiation as the 45 light source. The satellites will be in a low Earth orbit (LEO) 46 of ≈500 km altitude. In each orbit cycle (≈94 min), several 47 cloud-fields will be imaged. The clouds of major interest are 48 small liquid-phase clouds, typically hundreds of meters wide. 49 Therefore, ground resolution should be finer than ≈50 m at 50 nadir. 51 Classic remote-sensing assumes a plane-parallel geome-52 try [2], [3], [4], [5], [6], [7]. This effectively degenerates, 53 to some extent, 3-D RT to a 1-D form, as if effectively there 54 is neither horizontal heterogeneity nor RT. This assumption is 55 particularly invalid at edges of clouds. This leads to biased 56 retrievals of the cloud microphysics [2], [3], [4], [5], [6], [7]. 57 In particular, this leads to high uncertainties regarding the 58 microphysics and radiative effects of small clouds, regardless 59 of the sensing resolution. 60 To allow better understanding of microphysics and RT in 61 small clouds, a 3-D retrieval approach is essential. 3-D mod-62 eling currently receives growing attention [8], [9], [10], [11]. 63 Levis et al. [12], [13], [14], introduced 3-D scattering com-64 puted tomography (CT), based on the spherical harmonic dis-65 crete ordinate method (SHDOM) for RT [15]. Their method, 66 pySHDOM [16], retrieves cloud properties by fitting multiview 67 light intensity images to a physics-based forward model. 68 This is a generalization of CT to recover scattering by 69 passive sensing, relying only on the Sun as an illumination 70 source.  [7] present 126 a positive bias of 50%, mainly due to 3-D radiative effects. 127 Matar et al. [6] demonstrate dual-wavelength retrievals 128 using multiview polarized measurements, taken by the airborne 129 Observing System Including PolaRization in the Solar Infrared 130 Spectrum (OSIRIS). Here too, analysis used the plane-parallel 131 assumption. They emphasize the significance of a heteroge-132 neous column model, pointing to an error of more than 10% 133 solely due to a homogeneous column assumption. 134 Liquid-water-content (LWC) can be retrieved using cloud 135 radar reflectivity. Ebell et al. [39] estimate errors in retrieval 136 of the LWC to be larger than 60% for shallow nondrizzling 137 clouds. Zhu et al. [40] estimate the uncertainty in radar dual-138 wavelength retrievals of the LWC. They estimate an uncer-139 tainty of between 0.1 (g/cm 3 ) and 0.65 (g/cm 3 ) for shallow 140 clouds, decreasing with cloud thickness. 141 Operational methods do not consider 3-D RT effects, nor 142 do they perform 3-D volumetric retrieval of cloud parameters. 143 Levis et al. [18] derive and define a polarimetric 3-D scatter-144 ing tomography method for retrieval of cloud droplet micro-145 physics, using multiview multispectral image measurements, 146 based on vSHDOM. This method is the basis for our current 147 work. They present demonstrations based on the Airborne 148 Multiangle SpectroPolarimetric Imager (AirMSPI) [18], [41]. 149 Doicu et al. [42] presents another method for tomographic 150 cloud retrievals, also based on SHDOM. Doicu et al. [43] [11], [46], [47]. Forster et al. [26] 159 examine a basic limitation of retrieval. They explore a 160 region of the cloud in which microphysical properties can-161 not be resolved by scattering tomography, termed the veiled 162 core. Their demonstration, however, was without considering 163 polarization.

164
CloudCT offers very different spaceborne retrievals than the 165 operational state of the art. As in [18], CloudCT suggests 3-D 166 retrieval of cloud microphysics using scattering tomography. 167 Operational spaceborne instruments designed for global cov-168 erage, have retrieval spatial resolution that is coarse relative 169 to the clouds of interest. CloudCT will image clouds with 170 spatial resolution of finer than 50 m, and retrieve at volumetric 171 resolution of ≈50 × 50 × 50 m 3 .

172
In this article, we present a comprehensive method for 173 retrieval, highlighting adjustments and improvements which 174 we implement. We make advances needed for realistic space-175 borne perspective optical imaging. We model the imager as 176 having a polarized sensor as in [48]. We introduce a new 177 method for initialization of the retrieval medium, based on 178 a parameterized horizontally uniform model. This medium is 179 updated by an optimization process. Air molecular density is assumed to be known. Thus, 183 we do not try to retrieve air density. By cloud tomogra-184 phy, we retrieve the cloud-droplet microphysics. We assume 185 small warm clouds, i.e., liquid-phase clouds. The droplets 186 have spherical geometry 1 of variable radius values r . One 187 of the unknowns is the 3-D distribution of the LWC. 188 Let the droplet size distribution at a volume element be 189 n(r ) (1/μm) · (1/m 3 ). The droplet size distribution is 190 parameterized by a gamma distribution model as in [49].

191
The LWC is 259 For the retrieval of a set of N c separate clouds, a mean relative 260 error is  The temporal resolution of the retrieval depends on the 271 frame rate of data acquisition, while high synchronicity is 272 achieved. If synchronicity is not practical, temporal accumu-273 lation of data of up to ≈30 s may be applied [51].   1) The topocentric Earth (TE) coordinate system. The 282 origin of the TE system is on the Earth's surface. The TE 283 coordinate system is expressed as East, North, Up (ENU) 284 (Fig. 2). In our convention, the east axis is labeledŶ , 285 the northX, and the zenithẐ. 2) The Earth-centered (EC) coordinate system. The origin 287 of the EC system is at the Earth's center, as shown in 288 Fig. 3. The axes of this system areX earth ,Ŷ earth ,Ẑ earth . 289 AxisX earth is parallel toX,Ŷ earth is parallel toŶ and 290 Z earth is parallel toẐ. With respect to the EC system, the 291 center of the TE system is at Z earth = R earth , Y earth = 0, 292 X earth = 0, as shown in Fig. 3. The relation between EC 293 and TE is further detailed in Appendix I.

294
3) The camera coordinate system. The origin of the system 295 is at the center pixel of the sensor, as shown in Fig. 4. 296 The axes of this system areX cam ,Ŷ cam ,Ẑ cam . Axes 297 X cam andŶ cam are parallel to the pixel rows and 298 columns of the imager sensor, respectively. AxisẐ cam 299 aligns with the optical axis of the camera. The satellite 300 316 Letb 1 andb 2 be the normalized (to unit vector) pro-317 jections ofb 1 andb 2 on the camera plane, respectively.

318
The vectorb 2 aligns withX cam . The vectorb 1 does not 319 align withŶ cam for most of the pixels. There is a slight 320 deviation which can be expressed as

322
As an example, for a camera with a ≈5 • field-of-323 view (FOV), as intended for the CloudCT mission, 324 ≤ 0.143 • . Therefore, we assume the Stokes vector in 325 the camera coordinate system is practically identical to 326 that which is expressed in the pixel coordinate system.

338
Imagers are positioned by a satellite formation. The satel-339 lites' configuration is constant in all our demonstrations. A set 340 of ten imagers in separate satellites is considered. The satellites 341 are assumed to be in a trailing formation (string of pearls), 342 moving northward consecutively, as illustrated in Fig. 6. The 343 altitude of the satellites is R orbit = 500 km. The uniform 344 distance between each pair of neighboring satellites is 100 km, 345 on the orbit arc. The viewing angles are between −46 • and 346 39.3 • relative to the zenith.

347
The imager used in the demonstrations herein is a mono-348 chrome polarization camera, having a Sony IMX250MZR 349 sensor, and a spectral filter. The Sony sensor [48] has four 350 types of wire-grid polarizing filters which are formed on the 351 chip in blocks of four pixels, as illustrated in Fig. 7. Each 352 filter in a block has a polarization angle ψ ∈ [0 • ,45 • ,90 • ,135 • ], 353 relative toX cam in the image-sensor coordinate system (Fig. 4). 354 High spatial resolution is needed for retrievals at high vol-355 umetric resolution. State-of-the-art remote-sensing retrievals 356 are restricted by radiative smoothing [53], due to unresolved 357 horizontal fluxes. These fluxes have been neglected in tra-358 ditional methods due to the assumption of a plane parallel 359 cloud structure. A 3-D RT model obviates this limitation. 360 Therefore, the only limitation we consider is that of the imag-361 ing optics. In the demonstration herein, the imaging payload 362 optics lead to nadir resolution of 20 m. Tomography relies 363 on integrated measurements along multiangular ray paths. 364  positioned to view. Therefore, the chosen imager depends on 400 the solar illumination direction. The method for choosing the 401 most suitable satellite for the task is detailed in Appendix II. 402 In the setup which is described in Section III-B (see Fig. 6), 403 the satellite which scans the cloudbow is satellite number 9. 404 In the simulations presented in this article, a single imager 405 captures ten additional views within the scattering angles of 406 ≈135 • ≤ θ sca ≤ 150 • . On average, we find the cloudbow-407 scan decreases the mean errors LWC and r e by 2% and 408 1%, respectively. We find this is a consistent improvement. 409 However, further sensitivity studies may improve the choice 410 of scanned scattering angles and resolution.

411
Appendix II describes the calculations of the time-scale 412 for the cloudbow-scan. The time-scale may reach ≈30 s. 413 As demonstrated by [51], this is an acceptable time scale 414 regarding evolution of cloud droplet size. There is a linear relation between grayscale levels and photo-418 electrons detected at sensor pixels (See Appendix III). This 419 section describes the conversion of the expected number 420 of measured photo-electrons to Stokes vector components 421 (radiance units) in the meridian coordinate systems (see 422 Section III-A). 3 The pipeline of the conversion is illustrated 423 in Fig. 10.

424
Assume a polarization sensor, as described in Section III-B 425 (see Fig. 7). To allow evaluation of the Stokes vector per sensor 426 pixel, demosaicing is used [54], [55]. This way, each image 427 pixel k is associated with four expected numbers of measured 428 photo-electrons at a pixel, N measured          We denote this first method as H Typical . In addition to the 493 homogeneous model, we describe a horizontally homogeneous 494 parametric cloud model. The parameter values used for each 495 initialization are found by preoptimization of an initial cost 496 function, as we describe. We emphasize that these models are 497 suggested strictly for Stage 2. Once initialized, the optimiza-498 tion in Stage 3 is no longer restricted to either a homogeneous 499 or a horizontally homogeneous assumption.

501
Let Z 0 be the cloud-base height. The value of Z 0 can be 502 evaluated by methods of space-carving [57]. Assume the initial 503 cloud model has vertical monotonic profiles of LWC and r e 504 at altitude Z ≥ Z 0 within a cloud. Following [58], set: (24) 507 Here, α l and α r are parameters of the LWC and r e , respec-508 tively. Thus, = [α l , α r ]. We set LWC min = 0.0001 (g/m 3 ), 509 r e min = 2.5 μm. These settings are tested on large-eddy-510 simulation (LES) cloud data and found suitable for steady-state 511 cumulus convection.   Examples of the cost for a particular cloud, using the 542 monotonous model, are presented in Fig. 12. The cost D DoLP 543 has a better-defined minimum. This example indicates that 544 there is a potential advantage to use the DoLP for setting . 545 The different methods are summarized in Table I. An example 546 of initialization profiles of LWC and r e is plotted in Fig. 13. 547

V. PRECONDITIONING AND ALTERNATING OPTIMIZATION 548
This section relates to Stage 3: high-dimensional optimiza-549 tion of CT (see Section II-B, and Fig. 1). The optimization 550 attempts to solve a problem, which consists of unknowns 551 whose commonly used numerical scale differs by orders 552 of magnitude. In small clouds (up to a domain of 1 × 1 × However, we find that the preconditioning is cloud-sensitive. 562 Levis et al.
[18] assume a 1-D structure of r e retrieval. For r e 563 with 3-D spatial variation, we find preconditioning factors of 564 LWC = 10 and r e = 0.1 to be useful. 565 We also consider a different approach that requires no pre-  Table II.  mask, which is extracted from the 3-D extinction field of the 607 cloud data.

608
For each initialization method, LWC and r e are estimated 609 (see 7,8). A comparison (presented in Fig. 14) demonstrates 610 the superiority of M Stokes and M DoLP for initialization. These 611 methods use the monotonous cloud model (see 23, 24).

612
The difference in errors of M Stokes and M DoLP is very small, 613 with a slight advantage to M DoLP . For some of the clouds 614 sampled, these two methods result in the same parameters 615 for initialization of Stage 3. For the few clouds that yielded 616 significantly different parameters, we find an improvement in 617 favor of M DoLP . For this reason, it appears there is a benefit 618 to use M DoLP for initialization. The advantage of M DoLP is 619 in clouds having (vertical) optical-depth values below 30. 620 In a scattering medium, a lower optical-depth corresponds to 621 less multiple scattering. This is equivalent to a higher weight 622 of single scattering, in which polarization is significant. The 623 scientific goals of the CloudCT mission are focused on small 624 warm clouds. Therefore, it is a reasonable assumption that a 625 significant part of the actual data acquired in the mission has 626 the order of these low optical-depth values. 627 We find that the 1-D model assumed for r e (see 24) 628 proves to be a good assumption. As demonstrated in Fig. 14 Fig. 15 illustrates that a small 654 cloud rendered by an IPA approach is brighter than that  4 This orthographic projection is not used in the inversion. 5 A plane parallel atmosphere assumption also neglects the curved geometry of the atmosphere. For retrieval of clouds in a wide medium, this is significant. However, in the scale of the medium that we examine, this effect is negligible. As Forster et al.
[26] demonstrate, a veiled core of a cloud 663 limits tomography. We will continue to explore the error in 664 clouds of focus in CloudCT.

665
Retrieval accuracy depends also on the signal-to-noise 666 (SNR) ratio. A thorough analysis of the effects of electronic 667 noise, airlight, and stray light on the SNR will be discussed 668 in our future work.

669
In algorithmic aspects, as demonstrated, retrieval highly 670 depends on initialization. In addition, there is uncertainty 671 stemming from the process of space carving. In space-carving, 672 image processing and stereo imaging are used to set a domain 673 of interest. However, the cloud base is not directly viewed by 674 spaceborne imagers. So, there may be a higher uncertainty of 675 the cloud geometry there.

676
An additional biasing effect can be seen in the retrieval 677 of the LWC. This can be seen in the scatter plots presented 678 in Fig. 16. Here, a tendency of underestimation of LWC is 679 apparent. These errors will be further examined. We believe 680 that algorithmic advances will counter this error.

682
We introduce a comprehensive method for spaceborne 3-D 683 volumetric scattering-tomography of cloud microphysics. The 684 method is tailored for the CloudCT space mission, of 3-D 685 scattering tomography of warm clouds. It includes adjustments 686 to the pySHDOM 3-D microphysical scattering tomography. 687 The major adjustments are implementation of a realistic 688 polarized imager model, and a new initialization method. 689 We demonstrate the superiority of an initialization method 690 based on a parameterized horizontally uniform model, under 691 the constraints of the implemented imager model. For small 692 clouds with low values of vertical OD, we find that initializa-693 tion based on the DoLP has an advantage. For the CloudCT 694 mission which focuses on small clouds, this may be significant.

695
In addition, we suggest cloudbow scanning.

696
Future research toward the CloudCT mission will further 697 address the potential of the cloudbow-scan. The final orbit 698 plan will be taken into consideration, along with a statistical 699 analysis of the satellite pointing accuracy, maneuver capabili-700 ties, and payload frame rate. Section III-A). Both are illustrated in Fig. 3.

734
Let ξ sat be the zenith angle of the satellite, in the EC 735 coordinate system (Fig. 3). Let L be the distance (on the orbit 736 arc) between two adjacent satellites (L = 100 km as described 737 in Section III-B). The angular distance between two adjacent 738 satellites is 746 So, in the EC coordinate system In the TE coordinate system, the satellite coordinates are In TE Earth system (Fig. 3), the zenith angle of the 753 satellite is Here, we show that the time-scale for the cloudbow-758 scan may reach ≈30 s. As demonstrated in [51], this is an 759 acceptable time scale regarding cloud development. In the 760 simulations presented in this article, the range of the scanned 761 scattering angles in the cloudbow is 135 • ≤ θ sca ≤ 150 • . 762 We uniformly sample ten angles in this range.

763
Using the same convention as described in Figs. 2 and 3 764 (φ sat = 0 • ), the direction vector from the Sun to the domain is 765 The direction from the origin of the TE coordinate system to 767 the satellite is The angles ξ sat and θ sat (and generally, the azimuth angles φ sat ) 776 can be set to a specific angular range relative tod sun , specif-777 ically within the cloudbow. This requires inverting (40), and 778 extracting θ sat as a function of θ sca . We solve this numerically. 779 First, angle samples ξ sat are converted to angle samples θ sat 780 (using 33, 35 and 36). Then, by numerical search of (40), 781 we seek samples θ sat that provide the required samples of θ sca . 782 Let ξ sat be a zenith angle of the satellite, in the EC 783 coordinate system, sampled at much higher resolution than 784 ξ sat (Fig. 17). The goal of using angles ξ sat is to discretize the 785 position on the orbit path with high resolution. We generate 786 1310 samples of ξ sat ∈ [ξ sat [1], ξ sat [10]]. ten such angles, we extract the two extreme and denote them 799 by ξ 1 and ξ 2 (Fig. 17). These two extreme angles are used to 800 assess the time that is needed to scan the cloudbow by one 801 satellite.

802
Let v sat be the velocity of the satellite (Fig. 17). For orbit 803 radius (R earth + R orbit ), the velocity is where [65] μ g = 3.986004418 · 10 5 (km 3 /s 2 ) is the stan-806 dard gravitational parameter. In a low-Earth orbit, v sat ≈ 807 7.6 (km/s). The distance that a satellite travels (Fig. 17) on 808 the orbital arc between the two angles ξ 1 and ξ 2 is

810
The time T in which the satellite scans the cloudbow 811 angular region is

813
As an example, consider satellite 9 in the setup described 814 in Section III-B (see Fig. 6). Let it capture ten samples cloudbow-scan satellite is nadir observing (as in Fig. 8 824 There, the wavelength λ is between [λ 1 , λ 2 ]. Let us denote 825 spectral radiance at wavelength λ by I λ . It is calculated by 826 pySHDOM and has units of (W/m 2 · sr · nm). The camera 827 yields perspective projection. The perspective projection is 828 implemented in pySHDOM.

829
The camera system efficiency due to optical losses at wave-830 length λ is τ λ . Consider a camera with a lens of diameter D at 831 distance f from the focal plane. Lens distortions are assumed 832 to be known and compensated for. The camera is focused on 833 the object (i.e., clouds). For simplicity assume I λ is uniform 834 within the area of a pixel footprint on the cloud. Geometrically, 835 this is equivalent to uniform irradiance on the detector pixel 836 (area p 2 ). We consider the pixels to be close to the optical 837 axis of the camera. 6

838
To simulate a readout of the imager, we convert I λ within 839 to photo-generated electrons in the sensor, as follows. 7 Light 840 energy is converted to the expected number of photons at 841 wavelength λ by the factor (λ/h c) (photons/Joule). A pixel 842 on a sensor responds to the photons in a spectral band . The 843 pixel response depends on the sensor's quantum efficiency, 844 QE λ (electrons/photons). It is a measure of the probability 845 for a photo-electron to be created per incident photon with 846 wavelength λ. 847 Let t be the exposure time of the imager. The expected 848 number of photo-electrons that are created in a pixel is Ideally, to calculate the integral in (46), RT would be run 856 multiple times to calculate I λ in the spectral band , at a high 857 wavelength resolution. There is a common approximation that 858 simplifies these calculations [13], [61]. Instead of multiple 859 calculations of I λ within , in this approximation, RT is 860 simulated only once. LetĨ λ be the spectral radiance which is 861 a result of this single RT run i.e.,Ĩ λ is an approximation of I λ . 862 The fieldĨ λ is used instead of I λ in (46), and has the same 863 units as I λ . When λ is significant, the approximationĨ λ ≈ I λ 864 is valid if wavelength dependencies within a spectral band are 865 weak. This condition is met when narrow bands are considered 866 (e.g., up to 40 nm [13]). Under these assumptions, the optical 867 quantities can be spectrally averaged over . Then spectrally 868 averaged optical quantities are used as constants within . 869 Next, we describe how to expressĨ λ . Let I be the spectral

884
To simulate raw measurements, we introduce noise to N .

916
To convert the measured photo-electron count N measured to 917 measured radiance I measured , we use (49)