Retrieval of Cloud Liquid Water Using Microwave Signals from LEO Satellites: A Feasibility Study through Simulations

A novel approach, using low Earth orbit (LEO) satellite microwave communication links for cloud liquid water measurements, is proposed in this paper. The feasibility of this approach is studied through simulations of the retrieval system including a LEO satellite communicating with a group of ground receivers equipped with signal-to-noise ratio (SNR) estimators, a synthetic cloud attenuation field and a tomographic retrieval algorithm. Rectangular and Gaussian basis functions are considered to define the targeted field. Simulation results suggest that the proposed least-squares based retrieval algorithm produces satisfactory outcomes for both types of basis functions. The root-mean-square error of the retrieved field is around 0.2 dB/km, with the range of the reference field as 0 to 2.35 dB/km. It is also confirmed that the partial retrieval of the cloud field is achievable when a limited number of receivers with restricted locations are available. The retrieval outcomes exhibit properties of high resolution and low error, indicating that the proposed approach has great potential for cloud observations.


Introduction
Understanding the role of clouds is of great importance in many areas such as numerical weather prediction and climate studies, as clouds are not only vital in the hydrological cycle, but also crucial for the Earth's radiation budget. Various instruments have been developed and deployed, such as microwave radiometers, cloud radars and lidars, to measure cloud liquid water content (LWC) on a fine scale [1]. However, due to the high cost and scarce distribution of such instruments, the real-time global seamless high-resolution observation of clouds is still a challenge.
In this paper, we propose a novel cloud estimation approach as an auxiliary means to traditional methods, in which the attenuation of low Earth orbit (LEO) satellite microwave links is used to measure the LWC in clouds. Currently, the number of active LEO satellites in space is well above 1000, and several companies are pursuing the proposal for large constellations of LEO satellites to provide global broadband access, e.g., OneWeb's system with 720 LEO satellites for global coverage and SpaceX's system with more than 4000 LEO satellites [2]. Once this service is available to the public, it can be foreseen that millions of user terminals (ground receivers) will be deployed across the globe. Because of the high demand for bandwidth, Ku-band (12~18 GHz) and Ka-band (26.5~40 GHz) are usually used This model consists of a LEO satellite that passes over a cloud and multiple ground receivers equipped with electronically steerable phased array antennas [14]. The overpass time is determined by the orbital parameters of the satellite and the minimum elevation angle of the ground receivers. The satellite has a direct spot beam to the ground receivers, and the receivers use their antennas' steering capability to track the transit of the satellite. The SNRs at the ground receivers are estimated for the purpose of measuring the path-integrated attenuation.
During a single overpass, the signal link between a satellite and a ground receiver moves along a curved surface, thereby the attenuation caused by water droplets on this surface can be detected by the SNR estimation at the receiver. If we consider the scenario where the satellite passes directly over the receiver, the curved surface can be regarded as flat because the curvature is purely induced by the slight Earth rotation during the overpass and can be ignored. When multiple receivers are placed along the ground track of the satellite, the moving signal links will cover most of the vertical plane [6], and thereby offer the possibility of retrieving the vertical attenuation field. An illustration of the orbital geometry and the targeted attenuation field is shown in Figure 1. Depending on the satellite system, the surface coverage of a single beam could vary from 2800 km 2 to 72,000 km 2 [2]. An example of the −3 dB footprint contour is shown by the black dashed lines in Figure 1. Each satellite can achieve a data-rate above 20 Gbps [2], which means that it can easily serve 10,000 user terminals (assuming 2 Mbps bandwidth each) simultaneously within a beam.
The proposed method can be easily extended to 3D attenuation field retrieval with additional receivers placed in a two-dimensional array on the ground. Previous studies have shown that 3D tomographic attenuation field reconstruction can be achieved through the aggregation of a series of reconstructed vertical planes [6]. This paper will only focus on the retrieval of one vertical cloud attenuation field.  The proposed method can be easily extended to 3D attenuation field retrieval with additional receivers placed in a two-dimensional array on the ground. Previous studies have shown that 3D Atmosphere 2020, 11, 460 4 of 17 tomographic attenuation field reconstruction can be achieved through the aggregation of a series of reconstructed vertical planes [6]. This paper will only focus on the retrieval of one vertical cloud attenuation field.

Signal Attenuation Model
The signal attenuation model used here is the same as those in previous studies [6,15]. The signal power at a ground receiver is obtained by in which the power of the received signal and the power of the transmitted signal at the satellite are denoted by P r and P t , respectively. The power of noise is denoted by P n , which includes both the system noise and the sky noise. The channel gain in Equation (1) is determined by three major factors: ξ I is the path-integrated cloud attenuation, ξ F is the free space power loss (FSPL) and ξ C represents the attenuation due to all other factors such as the gain of antennas and effects of water vapor. Suppose that there are M signal samples available at the receiver over the duration the overpass of the satellite. Each sample is acquired from a different and very short (0.05 s in the simulations) time window. For the k-th sample, the receiver's SNR in decibels is obtained by We define A I (k) as the path-integrated cloud attenuation in decibels as and the free space loss A F (k) is obtained by thus, ρ(k) can be rewritten as Let the noise temperature of the receiving system and the sky temperature be denoted by T sys (k) and T sky (k), in Kelvin, respectively. In this case, noise power P n (k) is obtained by [16] where k B is the Boltzmann constant and B is the bandwidth of the communication system. We define parameter C(k) as and noise figure F n (k) as according to Equation (5), ρ(k) is then obtained by Atmosphere 2020, 11, 460 5 of 17 Parameter C(k) can be regarded as a baseline value for a receiver, which is mainly determined by the noise temperature of the receiving system, the antenna gain and the power of the transmitted signal. In this model, cloud is the primary source for T sky (k), which is generally obtained by [17] T sky (k) = G(k, Ω)T(Ω)dΩ where Ω is the solid angle, G(k, Ω) is the antenna gain pattern and T(Ω) is the noise temperature for a particular direction determined by Ω, which is obtained by where t m is the mean path temperature, in Kelvin, and A I is the total path attenuation for that direction, in dB. We assume that the beam pattern of the antenna is known so the relation between T sky and A I can be established using Equations (10) and (11). Meanwhile, A F (k) is a function of the satellite-to-ground path length L S and the signal frequency f , and is commonly obtained by [18] where the unit of L S (k) is in km and f in GHz. Since A F (k) can be regarded as known, Equation (9) indicates that A I (k) has the following relation with the SNR estimation: whereρ(k) is the estimated SNR at the receivers.

Retrieval Method
With A F (k) andρ(k) known, the retrieval process mainly involves handling the two unknown terms C(k) and F n (k) as well as solving the cloud attenuation field.
Consider a vertical rectangular attenuation field of interest. As the overpass time is very short for the LEO satellite (e.g., 289 s for a satellite at 1000 km of orbital height and a minimum elevation angle of 40 degrees for ground receivers), the attenuation field is regarded as static for the purpose of retrieval. It is described by a function of two variables R(x, z) representing the specific attenuation in dB/km at coordinates (x, z) km, where the x-axis represents the horizontal direction along the ground track of the satellite and the z-axis is the vertical direction starting from sea level. The relation between the path-integrated attenuation and the attenuation field function is obtained by [19] where L k denotes the signal line path for the k-th sample (see Figure 2). It is assumed that R(x, z) can be approximated by a linear combination of N basis functions b n (x, z) N n=1 : where w n is the weight for b n . Equation (14) then becomes Atmosphere 2020, 11, 460 6 of 17 in dB/km at coordinates ( , ) km, where the -axis represents the horizontal direction along the ground track of the satellite and the -axis is the vertical direction starting from sea level. The relation between the path-integrated attenuation and the attenuation field function is obtained by [19] ( ) = ( , ) , where denotes the signal line path for the -th sample (see Figure 2). It is assumed that ( , ) can be approximated by a linear combination of basis functions ( , ) : where is the weight for . Equation (14) then becomes We define and , represents the line integral of the -th basis function along the -th sample path. Once the basis functions have been determined, , can be obtained solely from the geometrical conditions of the satellite and the ground receiver at the time that the signal is sampled. As a result, , can be calculated offline and stored away for later use. We define and q k,n represents the line integral of the n-th basis function along the k-th sample path. Once the basis functions have been determined, q k,n can be obtained solely from the geometrical conditions of the satellite and the ground receiver at the time that the signal is sampled. As a result, q k,n can be calculated offline and stored away for later use.
If the path-integrated attenuation A I can be measured, the retrieval task then becomes a matter of solving the series weights {w n } N n=1 . However, the path-integrated attenuation is only measured indirectly through the SNR estimation. We rewrite Equation (13) as where V(k) is the sum of two known terms: Note that A F (k) can be calculated by Equation (12) but C(k) and F n (k) are unknown. For a similar problem in [6], C(k) is regarded as a constant and solved together with {w n } N n=1 using a least-squares method. However, in reality, parameter C(k) could change during one overpass [15,20], so forcing it to be a constant could introduce considerable errors. An improved approach [15] to eliminate C(k) before solving {w n } N n=1 was introduced as follows. Considering another sample link (k + ∆) which is very close in time to sample link k, Equation (18) becomes The difference between C(k + ∆) and C(k) is ignored because the two links are very close in time, Equations (18) and (20) yield Without prior knowledge of the path-integrated attenuation, noise figure F n (k) is completely unknown. One viable approach is to initially regard F n (k) as zero for the first least-squares computing. According to Equation (8), this means that T sky is zero, which is physically incorrect [21]. However, it has been shown that the retrieved {w n } N n=1 can be used to update F n (k) and in turn correct in an iterative manner [6]. In this paper, we adopt this iterative approach, as well as explore an Atmosphere 2020, 11, 460 7 of 17 alternative way of processing the first retrieval result directly to achieve satisfactory estimations, which is discussed in detail in Section 4. For the first iteration, we regard F n (k) as zero and consider Equation (21) in the matrix form where w is {w n } N n=1 and Q is the (M − ∆) × N differential matrix: We then solve Equation (22) by formulating the following linear least-squares problem where From the second iteration, F n (k) is updated using Equations (8) and (10), , and the least-squares problem (25) is run again to update {w n } N n=1 . Two types of basis functions are investigated in this paper. The first is a two-dimensional rectangular function obtained by where u is the width of the rectangular function and (X n , Z n ) is the center coordinate of the n-th function. Essentially, the rectangular basis functions offer a way to cover the field of interest with a regular grid. The integral q k,n is simplified to the length of the k-th link within the n-th unit square. The second type of basis function is a Gaussian basis function defined by [22] µ n (x, z) = exp where δ x,z,n is the distance from (x, z) to the center of the n-th basis function located at (X n , Z n ) and σ is the width of the Gaussian functions.
Note that no matter which basis function is applied, the differential matrixQ can be calculated offline before using the least-squares.

Synthetic Attenuation Field
The synthetic cloud data used in this paper were produced from an idealized simulation of a mesoscale convective system using the WRF model. More details about the simulation can be found in Ref. [12]. This cumulus type cloud has very large vertical development, the cloud base is approximately 5 km in height, with the top around 10 km. One vertical slice of the cloud field that covers 80 km in width and 0-12.5 km in height was selected as the reference field. This slice contains a very high content of water, which mainly consists of supercooled cloud water. Some areas can produce specific attenuation up to 2.35 dB/km, while there are other areas with relatively thin clouds. Very few rainwater droplets are present in this slice. A 320 × 50 grid is used to define the field with each unit square being 250 m by 250 m, i.e., both the horizontal and vertical resolutions are 250 m.
In the cloud LWC field generated by WRF, each unit square has a uniform cloud water mixing ratio given in g/g, which is defined as where ω w and ω a are the liquid water and dry air density in g/m 3 , respectively. Note that other particles such as ice crystals and snow are not considered as they have minimal effects on the attenuation of microwave signals around 30 GHz [8].
Dry air density is obtained by [23] ω a = 11.612p × 300 in which T c is the cloud liquid water temperature in Kelvin and p is the dry air pressure in kPa. We use the barometric law [24] to calculate p: where z is the altitude above sea level, P 0 is static pressure at sea level (101.325 kPa), g 0 is gravitational acceleration (9.80655 m/s 2 ), M e is molar mass of Earth's air (0.0289644 kg/mol) and R 0 is universal gas constant (8.3144598 J/(mol·K)). For each unit square in the cloud field, temperature T c is calculated according to the altitude of the center of the square, assuming that the temperature lapse rate L b is −0.0065 K/m from the sea level temperature T 0 of 288 K. Note that the constant lapse rate is an approximation for the cloudy atmosphere and will introduce an inaccuracy in the temperature profile. Nevertheless, it only serves as a means to translate the WRF-generated LWC field into the reference attenuation field, and using it does not affect the validation of the retrieval method per se. In order to generate the reference attenuation field, we use the equations provided by the International Telecommunications Union (ITU) [25], where the specific attenuation within a cloud is obtained by where K l is the cloud liquid water specific attenuation coefficient (details can be found on Page 2 in [25]) in (dB/km)/(g/m 3 ), f is the signal frequency in GHz. With the temperature and pressure profile ready for each unit square, ω w can be calculated using Equations (30) and (31) and the specific attenuation R can then be calculated by Equation (33) with the signal frequency f = 30 GHz. The reference attenuation field is shown in Figure 3a.
because is overestimated, as the attenuation field is overestimated in the first iteration. As a result, the retrieved fields will oscillate between iterations, but averaging the attenuation fields of two consecutive iterations will yield satisfactory results for rain field estimation [6]. Figure 3c shows the final retrieved field by averaging the third and fourth iterations. It can be seen that, for both heavy cloud areas where there is large vertical development and thin cloud areas, the retrieved field and the reference field have good agreement.

Simulation of the Estimated SNR
We considered a LEO satellite with a circular Keplerian motion trajectory of an orbital height of 1000 km and an inclination angle of 96 • . The trajectory is on the same vertical plane as the cloud attenuation field of interest. There are 91 receivers evenly placed at the ground level, with any two adjacent receivers being 1 km apart. The first and the last receivers are located at x = −5 km and x = 85 km in Figure 3a, respectively. With the assumption that the minimum elevation angle for the receivers is 40 • , it can be calculated that the overpass time is approximately 289 s. We use 400 signal samples (M = 400) during the satellite overpass time for each receiver. The samples are taken at fixed time intervals, i.e., the overpass time is evenly divided by the 400 samples.
The observation window for each sample is chosen as 0.05 s. According to previous results [6], the error induced by satellite movement in the observation window is significantly lower than the Cramer-Rao Lower Bound (CRLB) of the SNR estimation. Consequently, the path-integrated cloud attenuation A I and the free space loss A F are calculated only for the starting point of the observation window and assumed to remain the same within the window. The estimated SNRρ for each sample of each receiver is simulated in the following steps: • For each receiver, preset the constant C along a sine curve between 214 dB and 216 dB to simulate the differences across different receivers.

•
Calculate the path-integrated attenuation A I (k) and the free space loss A F (k) using Equations (12) and (14), respectively.

•
The sky noise temperature T sky (k) is computed according to the cloud attenuation from nine different directions and a simplified antenna beam pattern (Equations (10) and (11) with t m = 275 K, see [6] for details). Noise figure F n (k) is then calculated using Equation (8) with T sys = 150 K.

•
The true SNR is then calculated using Equation (9).

•
We assume a Binary Phase Shift Keying (BPSK) modulated signal with a bit rate of 10 Mbps and the estimated SNRρ is generated by adding the true SNR to a Gaussian random variable with zero mean and variance being equal to the CRLB.

Retrieval Using Rectangular Basis Functions
In this section, we use N rectangular basis functions obtained by Equation (27) to retrieve the attenuation field. Here, N is assumed to be 320 × 50 and width u for each function to be 250 m. In other words, the retrieved attenuation field is defined by a grid exactly the same as that of the reference field.
Theoretically, parameter ∆ for the differential approach in Equation (21) needs to be small so that the assumption that C(k) stays static can be met. However, simulations show that ∆ = 1 will generally make the differential distances (q k+∆,n − q k,n ) very close to zero and thus make the algorithm sensitive to noise and estimation errors. After some experiments, parameter ∆ is set at 3, which seems to produce the best retrieval outcomes. After calculating the differential matrixQ, a constrained least-squares algorithm using an interior-point method (lsqlin in Matlab, see [26,27]) is employed. The upper and lower bounds for w n are set to 5 and 0, respectively. In addition, we assume that there being no cloud in the bottom 320 × 20 unit squares is prior knowledge. In all of the following simulations, therefore, the weights for the basis functions below 5 km are considered zero and only the weights for the basis functions above 5 km (320 × 30 basis functions in this case) are retrieved. Figure 3b shows the retrieved attenuation field obtained by the least-squares algorithm. Note that the retrieved specific attenuation is greater than the reference value (overestimation) for many of the unit squares. This is because, in the first iteration of the algorithm, by assuming F n in Equation (13) to be zero, the impact of F n is included in the path-integrated attenuation A I and makes it effectively larger. The retrieved attenuation field can be used to update F n and therefore enable the second iteration of the least-squares. The second retrieved field will be underestimated because F n is overestimated, as the attenuation field is overestimated in the first iteration. As a result, the retrieved fields will oscillate between iterations, but averaging the attenuation fields of two consecutive iterations will yield satisfactory results for rain field estimation [6]. Figure 3c shows the final retrieved field by averaging the third and fourth iterations. It can be seen that, for both heavy cloud areas where there is large vertical development and thin cloud areas, the retrieved field and the reference field have good agreement.
Assuming that the ground receiver antenna has high directivity, we can simplify the relation between the sky noise and the path-integrated attenuation by assuming that the sky temperature picked up by the antenna is from one direction only, i.e., T sky (k) = t m 1 − 10 − A I (k) 10 . ( As a result, from Equations (8) and (34), we can conclude that F n increases with A I in a nonlinear manner. The nonlinear relation can be observed in Figure 4, where F n (red line), A I (blue dashed line) and T sky (black line) are plotted over the entire course of the overpass. The left panel is for the receiver located at 40 km and the right panel is for the receiver located at 60 km. The ratio of F n /A I is large for small A I (e.g., F n A I ≈ 1.8 when A I = 0.1 dB) and is small for large A I (e.g., F n A I ≈ 0.6 when A I = 6 dB). Assuming that the ground receiver antenna has high directivity, we can simplify the relation between the sky noise and the path-integrated attenuation by assuming that the sky temperature picked up by the antenna is from one direction only, i.e., As a result, from Equations (8) and (34), we can conclude that increases with in a nonlinear manner. The nonlinear relation can be observed in Figure 4, where (red line), (blue dashed line) and (black line) are plotted over the entire course of the overpass. The left panel is for the receiver located at 40 km and the right panel is for the receiver located at 60 km. The ratio of / is large for small (e.g., 1.8 when = 0.1 dB) and is small for large (e.g., 0.6 when = 6 dB). For the iterative approach that generates Figure 3c, the least-squares computing needs to be run four times, and each time it deals with 400 × 91 = 36,400 equations. The relation between and in the above can be utilized to reduce the number of the least-squares computing to one time only, thereby reduce the computing complexity of the retrieval process. We propose the following piecewise linear (PL) function to process the retrieved after the first least-squares computing: where ( ) = 2.8 for < 0.5 1.9 for 0.5 ≤ < 1 1.8 for 1 ≤ < 1.5 1.6 for 1.5 ≤ < 2 1.5 for ≥ 2 .
The value of ( ) for each interval is chosen empirically from multiple simulations to achieve the best retrieval outcome. The retrieved attenuation field by applying on Equation (15) is shown in Figure 3d, which is very similar to the retrieved field by averaging the third and fourth iterations.
In reality, the basis functions can only be an approximation of the actual field. Different types of basis functions will introduce different levels of errors into the retrieval model. To investigate the effects of larger rectangle basis functions, we consider the retrieved attenuation field represented by For the iterative approach that generates Figure 3c, the least-squares computing needs to be run four times, and each time it deals with 400 × 91 = 36, 400 equations. The relation between A I and F n in the above can be utilized to reduce the number of the least-squares computing to one time only, thereby reduce the computing complexity of the retrieval process. We propose the following piecewise linear (PL) function to process the retrieved {ŵ n } N n=1 after the first least-squares computing: where The value of B(x) for each interval is chosen empirically from multiple simulations to achieve the best retrieval outcome. The retrieved attenuation field by applying w n N n=1 on Equation (15) is shown in Figure 3d, which is very similar to the retrieved field by averaging the third and fourth iterations.
In reality, the basis functions can only be an approximation of the actual field. Different types of basis functions will introduce different levels of errors into the retrieval model. To investigate the effects of larger rectangle basis functions, we consider the retrieved attenuation field represented by 160 × 25 rectangular basis functions for the second set of simulations. The width of each function (u) is 500 m, i.e., each unit square is now four times as big as before in area. We use the same PL function to process the retrievedŵ n . The final retrieved attenuation field is shown in Figure 3e. It can be seen that similar results are achieved but with poor spatial resolutions.
As a global index of the performance for the retrieval, the root-mean-square error (RMSE) between the retrieved attenuation field and the reference field has been considered, which is defined by where W denotes the total number of comparison points, R i and R i are the specific attenuations at a particular comparison point for the reference field and the retrieved field, respectively. For the following results, the comparison points are assigned to the centers of the top 320 × 30-unit squares in the reference field. The RMSE for the retrieved attenuation field after the first iteration (Figure 3b) is 0.323 dB/km. Averaging the third and four iterations (results shown in Figure 3c) reduces the RMSE to 0.101 dB/km and the PL function (results shown in Figure 3d) also brings the RMSE down to 0.107 dB/km. The RMSE for the retrieved attenuation field with enlarged rectangular basis functions (Figure 3e) is 0.215 dB/km.

Retrieval Using Gaussian Basis Functions
In this set of simulations, we replace the 160 × 25 rectangular basis functions in Section 4.2 with the same number of Gaussian basis functions, with the center location of each function unchanged. The width of the Gaussian functions (σ) is set to 1.5 km.
To compute the integral in Equation (17), let us divide the sample link k into a large number (S) of segments of the same length l. Each segment is so small in relation to the width of the Gaussian function that µ(x, z) can be regarded constant for anywhere on a segment. Consequently, Equation (17) becomes where x k,i , z k,i is the coordinates of the center of the i-th segment of sample link k.
After the differential matrixQ is calculated, the same least-squares algorithm is employed to retrieve {ŵ n } N n=1 . However, intuitively, a different PL function B(x) needs to be adopted because of the different nature of the Gaussian basis function, compared with the rectangular basis function. Multiple simulations suggest that a simple PL function such as the following would generate satisfactory results: A contour map of the final retrieved field after the PL function is applied is shown in Figure 5a. It has good agreement with the reference field. The RMSE of this field is 0.207 dB/km, which is smaller than the rectangular case in Figure 3e.
Examining the change in specific attenuation along a vertical line shows in detail the performance of different retrieval methods. For instance, Figure 5b is for a vertical line located at x = 21 km (shown by a black dashed line in Figure 5a) where the cloud is relatively thin. The green line shows how the specific attenuation changes with height in the reference field. The blue and red lines are for the retrieved fields using the PL function with respectively small and large rectangular basis functions. For the retrieved field with Gaussian basis functions (Figure 5a), the specific attenuation is plotted in black. The same plots are done for heavier clouds (vertical line located at x = 45 km, marked by a black dashed line in Figure 5a), shown in Figure 5c. For all three types of basis functions, the retrieved lines are able to closely follow the reference line and the level of error is the lowest for small rectangle basis functions. When the Gaussian basis functions are employed, as they are continuous, continuity is guaranteed in the retrieved field. The peak of the Gaussian retrieved line is generally lower than that of the reference line. This is partly due to the larger value of B(x) chosen to achieve a lower overall RMSE and partly due to the number of basis functions of the retrieved field being less than the number of grids in the reference field.
Atmosphere 2020, 11,460 13 of 17 continuous, continuity is guaranteed in the retrieved field. The peak of the Gaussian retrieved line is generally lower than that of the reference line. This is partly due to the larger value of ( ) chosen to achieve a lower overall RMSE and partly due to the number of basis functions of the retrieved field being less than the number of grids in the reference field. In addition to the RMSE (Equation (37)), the Pearson correlation coefficient (PCC) is also employed to compare the retrieval results across different methods. PCC is defined by where is the mean of and is the mean of . The scatter plots of the retrieved specific attenuations against the reference values at the comparison points can be found in Figure 6. Panel (a) is for the retrieved field by the averaging method (the result in Figure 3c), and the PCC is 0.986. Panel (b) is for the retrieved field by the PL function (the result in Figure 3d), with the PCC being 0.973. There are some horizontal gaps in the plot, indicating that discontinuity is present in the value range of the retrieved field ( ). This is because the retrieved after the first least-squares computing was reduced unevenly by the PL function. Panel (c) shows the scatter plot for the retrieved field in Figure 3e, and the PCC is 0.883. The scatter plot for the retrieved field with Gaussian basis functions (Figure 5a) is in panel (d), and the PCC is 0.892. In addition to the RMSE (Equation (37)), the Pearson correlation coefficient (PCC) is also employed to compare the retrieval results across different methods. PCC is defined by where R is the mean of R i and R is the mean of R i . The scatter plots of the retrieved specific attenuations against the reference values at the comparison points can be found in Figure 6. Panel (a) is for the retrieved field by the averaging method (the result in Figure 3c), and the PCC is 0.986. Panel (b) is for the retrieved field by the PL function (the result in Figure 3d), with the PCC being 0.973. There are some horizontal gaps in the plot, indicating that discontinuity is present in the value range of the retrieved field (R i ). This is because the retrieved {ŵ n } N n=1 after the first least-squares computing was reduced unevenly by the PL function. Panel (c) shows the scatter plot for the retrieved field in Figure 3e, and the PCC is 0.883. The scatter plot for the retrieved field with Gaussian basis functions (Figure 5a) is in panel (d), and the PCC is 0.892.

Effects of Larger Spacing among Receivers
Simulations were also carried out for more sparsely distributed receivers. For instance, Figure 7a shows the retrieval result when there are 37 receivers and the distance between any adjacent two is 2.5 km. The same rectangular basis functions and the same PL function for Figure 3d are employed.
The number of equations will not be sufficient for the least-squares algorithm to function correctly if the number of signal samples stays at 400. Hence, M is increased to 800 for each receiver. Note that this only uses 0.05 × 800/289 = 13.84% of the satellite signals. The RMSE and the PCC are 0.135 dB/km and 0.956, respectively. While it can be seen that the retrieval outcome is not as good as that in Figure 3d with 91 receivers, key features of the cloud can still be identified.

Effects of Larger Spacing among Receivers
Simulations were also carried out for more sparsely distributed receivers. For instance, Figure  7a shows the retrieval result when there are 37 receivers and the distance between any adjacent two is 2.5 km. The same rectangular basis functions and the same PL function for Figure 3d are employed. The number of equations will not be sufficient for the least-squares algorithm to function correctly if the number of signal samples stays at 400. Hence, is increased to 800 for each receiver. Note that this only uses 0.05 × 800/289 = 13.84% of the satellite signals. The RMSE and the PCC are 0.135 dB/km and 0.956, respectively. While it can be seen that the retrieval outcome is not as good as that in Figure 3d with 91 receivers, key features of the cloud can still be identified.
When the number of receivers is reduced to 19, which means that two adjacent ones are 5 km apart, the number of signal samples for each receiver is increased to 1600 in order to obtain a sufficient number of equations. Note that this only uses 0.05 × 1,600/289 = 27.68% of the satellite signals. The retrieval result (Figure 7b) further degrades in spite of the increase of samples. The RMSE and the PCC are 0.210 dB/km and 0.891, respectively. For both simulations with 37 and 19 receivers, further increasing the number of samples does not improve the performance of the retrieval. This indicates that the density of the receivers plays a more important role than the number of samples in the accuracy of the retrieval results.

Partial Retrieval
The variations in cloud horizontal scale are often quite significant [28]. Some cloud covers stretch over hundreds of kilometers and it may be difficult to measure the entire field at the same time. It would be of interest to see how the retrieval method would perform if only part of the cloud is reachable by the signal links. Thus, simulations were carried out for scenarios where ground receivers are restricted to smaller areas. For example, 11 receivers placed from = 35 km to = 45 km When the number of receivers is reduced to 19, which means that two adjacent ones are 5 km apart, the number of signal samples for each receiver is increased to 1600 in order to obtain a sufficient number of equations. Note that this only uses 0.05 × 1600/289 = 27.68% of the satellite signals. The retrieval result (Figure 7b) further degrades in spite of the increase of samples. The RMSE and the PCC are 0.210 dB/km and 0.891, respectively. For both simulations with 37 and 19 receivers, further increasing the number of samples does not improve the performance of the retrieval. This indicates that the density of the receivers plays a more important role than the number of samples in the accuracy of the retrieval results.

Partial Retrieval
The variations in cloud horizontal scale are often quite significant [28]. Some cloud covers stretch over hundreds of kilometers and it may be difficult to measure the entire field at the same time. It would be of interest to see how the retrieval method would perform if only part of the cloud is reachable by the signal links. Thus, simulations were carried out for scenarios where ground receivers are restricted to smaller areas. For example, 11 receivers placed from x = 35 km to x = 45 km (marked by the green bar at the bottom of Figure 8a) are used for retrieving the attenuation field. We firstly identify the targeted area (marked by the red dashed line in Figure 8a) and define it using 120 × 30 (u = 250 m) rectangular basis functions. The final retrieved field after the PL function is shown in Figure 8b. It can be seen that the middle section directly above the receivers has been retrieved successfully while for other areas the retrieved field is erroneous. The limited number of receivers with restricted locations cannot provide enough information to retrieve the left and right sections of the targeted area. Secondly, we redefine the target area (marked by the purple dashed line in Figure 8a) with 30 × 15 Gaussian basis functions. The final retrieved field after the PL function is shown in Figure 8c. The retrieved field in the middle section also agrees with its counterpart in the reference field, but the left area outside the 35 km mark and the right area outside the 45 km mark have a high level of errors. The retrieved attenuation field with 37 receivers (distance between two adjacent receivers being 2.5 km) using small basis functions and the same PL function. Panel (b): The retrieved attenuation field with 19 receivers (distance between two adjacent receivers being 5 km).

Partial Retrieval
The variations in cloud horizontal scale are often quite significant [28]. Some cloud covers stretch over hundreds of kilometers and it may be difficult to measure the entire field at the same time. It would be of interest to see how the retrieval method would perform if only part of the cloud is reachable by the signal links. Thus, simulations were carried out for scenarios where ground receivers are restricted to smaller areas. For example, 11 receivers placed from = 35 km to = 45 km (marked by the green bar at the bottom of Figure 8a) are used for retrieving the attenuation field. We firstly identify the targeted area (marked by the red dashed line in Figure 8a) and define it using 120 × 30 ( = 250 m) rectangular basis functions. The final retrieved field after the PL function is shown in Figure 8b. It can be seen that the middle section directly above the receivers has been retrieved successfully while for other areas the retrieved field is erroneous. The limited number of receivers with restricted locations cannot provide enough information to retrieve the left and right sections of the targeted area. Secondly, we redefine the target area (marked by the purple dashed line in Figure  8a) with 30 × 15 Gaussian basis functions. The final retrieved field after the PL function is shown in Figure 8c. The retrieved field in the middle section also agrees with its counterpart in the reference field, but the left area outside the 35 km mark and the right area outside the 45 km mark have a high level of errors. The results above indicate that partial retrieval is possible when a limited number of receivers or signal links are available. Because most of the signal links have high elevation angles, the best retrieval performance is found in the area directly above the receivers. It is also implied that tomographic retrieval can be carried out in a sectioned manner. For a targeted area within the attenuation field, only the signal links relevant to it are needed for its retrieval.

Conclusions
In this paper, a novel approach of using LEO satellite microwave signal links for cloud measurements is proposed. The feasibility of the approach is investigated using simulations, in which a synthetic cloud attenuation field is generated and the SNR estimation at the ground receivers is simulated. We use two types of basis functions to define the targeted field and employ a constrained least-squares algorithm to retrieve it. We also propose a piecewise linear (PL) function to quickly correct the overestimation issue due to the cloud LWC induced sky noise in the initial least-squares computing. Furthermore, partial retrieval can be achieved when only a limited number of receivers with restricted locations are available. The fact that the retrieval is satisfactory for both thin and heavy cloud areas suggests that the proposed approach has great potential. Its high resolution and low level of error imply that it would be very useful to the observations of the micro-and macro-physical details of nonprecipitating clouds.
Many aspects of the proposed approach still need further investigations. Testing the model with other cloud types, especially the ones primarily in the lower troposphere, is underway. Moreover, studies on precipitating clouds will be carried out, in which the main challenge is to distinguish the attenuation field induced by cloud liquid water from that induced by the rain water.