High-Resolution L-Band TomoSAR Imaging on Forest Canopies with UAV Swarm to Detect Dielectric Constant Anomaly

A rigorous TomoSAR imaging procedure is proposed to acquire high-resolution L-band images of a forest in a local area of interest. A focusing function is derived to relate the backscattered signals to the reflectivity function of the forest canopies without resorting to calibration. A forest voxel model is compiled to simulate different tree species, with the dielectric constant modeled with the Maxwell-Garnett mixing formula. Five different inverse methods are applied on two forest scenarios under three signal-to-noise ratios in the simulations to validate the efficacy of the proposed procedure. The dielectric-constant profile of trees can be used to monitor the moisture content of the forest. The use of a swarm of unmanned aerial vehicles (UAVs) is feasible to carry out TomoSAR imaging over a specific area to pinpoint potential spots of wildfire hazards.


Introduction
The morphology of a forest, including canopy height, vertical structure and spatial distribution, provides useful information for ecological protection [1], biomass estimation [2], the monitoring of biodiversity [1] and the global carbon cycle [3]. More incidences of forest fires in recent years may be correlated to the trend of global warming [4]. Regular monitoring of forest morphology may help in predicting potential wildfire spots for emergency preparation.
Wildfire risk increases as the soil dehydrates [4] or more dry bushes pile up [5]. L-band and C-band radars have been used to estimate the moisture content in the leaves, branches and trunks of trees [6], by detecting their dielectric constants [7].
In [8], an area-based approach was proposed to estimate different forms of deadwood. Light detection and ranging (LiDAR) has been used to retrieve 3D data of the forest structure, surface fuel depth, coverage and canopy density [9]. Synthetic aperture radar tomography (TomoSAR) in L-and P-bands has been used to reconstruct the 3D forest morphology via volume backscattering [10].
The soil and trunks are dominant scatterers in the P-band, while small branches are dominant scatterers in the L-band [11]. The backscattered signals in the L-band were reported to have a more uniform distribution in elevation [12]. In [13], an airborne Pand L-band TomoSAR technique using a Capon estimator was proposed to reconstruct a tropical forest, and the results were validated with small-footprint LiDAR (SFL) data. In [1], an airborne L-band TomoSAR technique was proposed to retrieve the vertical profile of a forest. In [14], an airborne L-band TomoSAR technique using wavelets was applied to survey the forest morphology, and its performance in conjunction with inverse algorithms of Fourier beamforming, Capon and compressive sensing, respectively, was analyzed.
To acquire TomoSAR images with high fidelity, the baseline aperture, baseline spacing, number of tracks and elevation ambiguity must be well orchestrated [15]. A method of coordinating baselines at nonuniform and uniform spacings was proposed by judging the corresponding point spread function [16]. As nonuniform baseline spacings may yield poor elevation resolution, a baseline aperture interpolator was proposed [2] to suppress sidelobes of the point spread function.
Inverse algorithms used for TomoSAR imaging of the reflectivity profile in the target volume can be roughly categorized into model-based (parametric), model-free (nonparametric) and hybrid [14]. Model-based algorithms include multiple signal classification (MUSIC) and covariance matching estimator (COMET). Model-free algorithms include Fourier beamforming (FB) and Capon beamforming (CB) [17]. In [18], a model-free iterative adaptive approach was shown to deliver finer elevation resolution than Fourier beamforming and Capon beamforming.
Compressive sensing (CS) techniques have been proposed to enable sparse SAR imaging methods [19,20]. In [21], a sparse SAR imaging method based on 2,1 -norm regularization was proposed to reduce the ghost-target phenomenon caused by azimuth periodic block sampling. In [22], an incremental SAR imaging method was proposed to improve the image quality and computational efficiency by reducing the azimuth ambiguity. In [23], a weighted sparse reconstruction method was proposed for single-channel SAR imaging on maritime targets to mitigate the azimuth ambiguities and improve the resolution.
In [16][17][18], several TomoSAR imaging methods were developed in terms of a normalized impulse response [24,25]. In [16], the effects of track number and track configuration on the acquired forest images were studied. It was reported that at least seven tracks were required to suppress the sidelobe level of the impulse response below −6 dB, resulting in a vertical resolution of 6 m. In [11], a weighted covariance fitting-based iterative spectral estimator (WISE) was proposed to retrieve the canopy height models and the above-ground biomass, achieving finer resolution than Capon and weighted-CS methods could deliver. In [14], the pros and cons of Capon beamforming, Fourier beamforming and CS algorithm were compared. In [18], a robust iterative adaptive approach (RIAA) in conjunction with a weighted least-squares criterion was proposed to conduct TomoSAR imaging with a small aperture of 5-30 m.
A vivid tree model is useful to validate a TomoSAR imaging method with fine spatial resolution. In [12,26], a forest was modeled as a reflectivity profile composed of a ground layer and a canopy layer. Three-dimensional tree generator software such as Arbaro (1995) [27], OnyxTREE (2020) [28] and AMAPstudio (2014) [29] were commonly used in LiDAR researches [30]. Arbaro is a rule-based growth algorithm [31], with which the trunks, branches and leaves of trees are transformed to a voxel model. A voxel model with a resolution of the order of a decimeter will be useful for pinpointing potential spots of wildfire hazards.
In this work, a rigorous L-band TomoSAR imaging method is proposed to reconstruct the dielectric-constant profile of a forest. A focusing function is derived to relate the backscattered signals to the reflectivity function along an elevation line segment at a specific azimuth and depth, without resorting to calibration. The profile of the reflectivity function along the elevation line segment is reconstructed with the inverse methods of compressive sensing (CS), Fourier beamforming (FB), multiple signal classification (MUSIC), amplitude and phase estimation (APES), and Capon, respectively. The simulated TomoSAR images acquired with these five inverse methods are analyzed in different scenarios under different signal-to-noise ratios.
SAR imaging tasks carried out on an unmanned aerial vehicle (UAV) were reported to achieve a fine spatial resolution of 10 cm [32]. The UAV navigation and trajectory for TomoSAR imaging has been studied [33]. A swarm of UAVs can be flexibly deployed to acquire high-resolution TomoSAR images of forest canopies, with the proposed TomoSAR imaging method.
The rest of this work is organized as follows. The forward problem and signal model are presented in Section 2, the TomoSAR imaging methods are presented in Section 3, the simulation results of acquired images are analyzed in Section 4 and some conclusions are drawn in Section 5. Figure 1 shows the schematic of electromagnetic scattering by a dielectric object in V with distribution of the refraction index n(r ). An incident field U i (r, τ) upon the object induces a scattered field U s (r, τ), resulting in the total field U(r, τ). The origin is arbitrarily located within the volume V, the observation point P is given byr =rr p +ŷy p +ŝs p and a point in volume V is given byr =rr +ŷy +ŝs. The (r, y, s) coordinate system is related to the conventional (x, y, z) coordinate system by a rotation about the y-axis. Figure 1. Scattering of electromagnetic wave by dielectric object in V [34], P atr =rr p +ŷy p +ŝs p is the observation point, and T atr =rr +ŷy +ŝs denotes a point in volume V.

Forward Problem and Signal Model
By taking the first-order Born approximation that U(r, τ) U i (r, τ) and the weak scattering assumption that |U s (r, τ)| |U i (r, τ)|, the scattered field can be represented as [34] where τ = τ − |r −r |/c. Linear frequency modulation (LFM) pulses are periodically emitted from the radar to the object, with the waveform of P(τ) = rect(τ/T r )e j(2π f 0 τ+πK r τ 2 ) , where T r is the pulse duration, and rect(τ/T r ) is a window function of duration T r . Thus, the incident field can be represented as P τ − |r −r |/c |r −r | and the scattered field is reduced to where ζ = τ − 2|r −r |/c, is the reflectivity function [34], and f 0 = f 0 + K r (τ − |r −r |/c), with τ = τ − |r −r |/c.

Range Compression
The backscattered signal in (2) is first demodulated to obtain the baseband signal U b (r, τ) = U s (r, τ)e −j2π f 0 τ . A matched filter P mf (u) = rect(u/T r )e −jπK r u 2 [35] is then convolved with the baseband received signal U b (τ,r) to give [34] where Q b (|r −r |, τ) = Q r (|r −r |, τ)e −j4π f 0 |r−r |/c [35] and is called the range factor [34]. Figure 2 shows the schematic of azimuth sampling along track n. The sampling position P nm is located atr =r nm , with −N a /2 ≤ m ≤ N a /2 and uniform spacing d. The range between P nm and a target point atr can be approximated as

Azimuth-Depth Compression
Thus, the phase term in (4) is approximated as where Ψ nm (r ) = e −j2πm 2 d 2 /(λ|r n0 |) . Let us define an azimuth reference function Φ n(k−m) = e j2π(k−m) 2 d 2 /(λ|r n0 |) , then multiply it to both sides of (4) and sum over −N a /2 ≤ m ≤ N a /2 to give where with sinc(x) = sin(πx)/(πx) and ζ nm = τ − 2|r nm −r |/c. By making the approximation ζ nm ζ n0 , (9) is further reduced to Figure 3 shows the amplitude of ξ nk (τ,r ) in (10), with −20 ≤ y ≤ 20 and −20 ≤ r ≤ 20. The magnitude of ξ nk (τ,r ) manifests a conspicuous peak at r = −10 m and y = 10 m (k = 20), similar to a conventional point spread function or a delta function. Since ξ nk (τ,r ) is concentrated about r = r , (8) is further approximated as where D r and D y are the spatial resolutions in range (r) and azimuth (y), respectively, S nk (r ) denotes a sequence of voxels along a line parallel to the s axis viewed from track n, with r = r and y = kD y . Figure 4 shows the schematic of TomoSAR imaging. Without loss of generality, all the tracks are assumed to be at the same height H. The side-looking angle from platform P 0 along the master track at η = 0 to point O is θ . The b axis is perpendicular to P 0 O and passes through P 0 . The s axis is parallel to the b axis and passes through O. . Schematic of TomoSAR imaging with platforms P 0 and P n along master track and track n, respectively [14], which are separated by perpendicular baseline b n⊥ , parallel baseline b n and parallel distance B n .

TomoSAR Imaging Methods
The backscattered signal from the qth voxel is received along the nth track and demodulated to the baseband as where ∆V q is the volume of the qth voxel, and η m = md/v p . The received baseband signals from all the Q voxels are given by By applying range compression and azimuth-depth compression as prescribed in the last Section, we have which is the acquired single-look complex (SLC) image along track n at y = kd, where (11) is used and τ n (r ) = 2|r n0 −r |/c (15) is the round-trip time betweenr n0 andr .
In conventional TomoSAR imaging methods, the original 3D inverse problem is decomposed into multiple inverse problems of solving a 1D SLC image as [24,25] where g n is the SLC image along track n, s is the coordinate along an elevation line segment at specific azimuth y and depth r , γ(s) is the reflectivity function to be reconstructed and φ(s) is a phase function pertinent to track location and coordinate s.
In [24], a spectral estimation for coherent adaptive nulling (SPECAN) algorithm was proposed to acquire SAR images at medium-to-low resolution. The received signal was modeled as a convolution of reflectivity function and a deramped phase, usually involving some calibration or normalization. In [25], a 2D point spread function (PSF) was proposed, which was later approximated as a 2D Dirac delta function to simplify the formulation. The SLC image acquired along track n was represented as a convolution of the 2D PSF and the reflectivity function multiplied with a phase function accounting for the round-trip range between the observation point and target point.
In our procedure, a rigorous scattering theory is used to derive the backscattered signals from the target volume. A focusing function at specific azimuth y and depth r is derived to relate the backscattered signals to the reflectivity function, without resorting to calibration or normalization. The focusing function, ξ nk (τ n (r ),r ), manifests a similar feature of a point spread function in the azimuth-depth plane, as in [25].

Covariance of Deramped Images
Next, multiply g nk (τ n (r )) with e j4π f 0 |r n0 |/c to derive a deramped imageg nk (τ n (r )), then compute the covariance between the two SLC deramped images along tracks n and n , respectively, as [36] R knn (r ) = E{g nk (τ n (r ))g * n k (τ n (r )) = D 2 r D 2 y S nk (r ) ds n S n k (r ) ds n E{γ(r n )γ * (r n )} e −j4π f 0 (|r n0 −r n |−|r n0 |)/c e j4π f 0 (|r n 0 −r n |−|r n 0 |)/c ξ nk (τ n (r ),r n )ξ * n k (τ n (r ),r n ) (16) wherer n ∈ S nk (r ) andr n ∈ S n k (r ). Then, the autocorrelation of the reflectivity function is approximated as [36] E{γ where δ(r n −r n ) is the Dirac delta function. In practice, S nk (r ) and S n k (r ) almost overlap with each other, and the difference between the look angles, ∆θ nn = θ n − θ n , cannot be ignored. Thus, (16) is reduced to [14,37] R knn (r ) = S nk (r ) F knn (s n , r )e −jκ nn s n ds n (18) where F knn (s n , r ) = D 2 r D 2 y |γ(r n )| 2 ξ nk (τ n (r ),r n )ξ * n k (τ n (r ),r n ) is called the vertical reflectivity profile,r n =ŝs n +ŷkd +rr , and κ nn = 4π∆θ nn /λ 0 is the vertical wavenumber, with 0 ≤ n, n ≤ N p . The covariances between all possible pairs of SLC deramped images are compiled into a covariance matrix as [38] [14]. Similarly, all the vertical wavenumbers are compiled into a vertical wavenumber matrix as The value of F knn (s n , r ) is insensitive to n or n ; hence, we select a set of grid points (voxels) along a common s axis in a vectors = [s 1 , s 2 , · · · , s L ] t , independent of n or n . The values of F knn (s n , r ) at the grid points {s n } are put in a vectorF k (r ) = [F k (s 1 , r ), F k (s 2 , r ), · · · , F k (s L , r )] t L×1 , which is also independent of n or n . Define a steering matrix on the grid points {s n } asĀ(s) = e −jKs 1 , e −jKs 2 , · · · , e −jKs L (N p +1) 2 ×L [14].
The backscattered signals are related to the vertical reflectivity profile, F knn (s n , r ), which is the counterpart of γ(s) in [39] and σ(x, γ) in [19]. The focusing function proposed in this work is the counterpart of the point spread function or delta function in conventional TomoSAR imaging methods. Since the focusing function is derived from the scattering theory, no calibration or normalization on the reconstructed reflectivity function is required.
The main contribution of our work is to establish a rigorous formulation of backscattered signals in terms of the reflectivity function via the scattering theory, without resorting to the empirical point spread function or delta function. The focusing function manifests a similar feature of PSF yet retains all of the original information. After range compression, azimuth-depth compression, co-registration and covariance among SLC images, an inverse problem in (21) is formed, which is solved by applying various inverse methods, including compressive sensing [40], Fourier beamforming [14,37], multiple signal classification (MUSIC) [19], amplitude and phase estimation (APES) [41] and Capon [14]. The acquired distributions of the dielectric constant in the target volume is analyzed to validate our procedure and compare the pros and cons of different inverse methods.

Compressive Sensing Method
The unknownsF k (r ) in (21) at L grid points along s direction are relabeled asū(s) = [u 1 , u 2 , · · · , u L ] t . Equation (21) is then solved with the compressive-sensing method as [40] u(s) = arg min whereĀ(s) plays the role of a dictionary. By substituting F knn (s n , r ) in (19) withũ(s), the reflectivity |γ(r n )| is obtained without resorting to calibration or normalization. Then, the dielectric constant can be estimated by using (3). Equation (22) is equivalent to [14] u(s) = arg min where f 1 and f 2 are L 1 -norm and L 2 -norm, respectively, off , and λ s is a hyperparameter [40]. The CVX tool can be applied to solve (22) or (23). The reconstructed reflectivity profile is sensitive to the chosen value of λ s , which is varied between 0.1 and 1. Hence, (22) is adopted in this work.

Fourier Beamforming (FB) Method
The SLC imagesg nk (τ n (r ))'s acquired along all the tracks are co-registered at y = kd and r = R 0 + r to constitute a measurement vector where τ 0 (r ) is the round-trip time between the master track andr , as defined in (15). The phase terms in (14) along all the tracks constitute a steering vector as Without loss of generality, the baseline of the master track is set to b 0 = 0. The steering vectors at all the grid points along the s axis constitute a steering matrix as A = [ā(s 1 ),ā(s 2 ), · · · ,ā(s L )] (N p +1)×L By using (24) and (25), the signal model is rephrased as g =Ā(s) ·ū(s) (26) whereū(s) = [u 1 , u 2 , · · · , u L ] t is the vertical reflectivity profile sampled at L grid points along the s axis. The image acquired with the FB method is [37] u F (s) is the covariance matrix of SLC images and

Multiple Signal Classification (MUSIC) Method
The covariance matrixR among SLC images can be decomposed asR =Ū ·D ·Ū † , whereD is the eigenvalue matrix sorted in descending order, namely, λ 1 ≥ λ 2 ≥ · · · ≥ λ N p −1 ≥ λ N p . The eigenvectors corresponding to the largest K eigenvalues are compiled into a signal matrixS, and the remaining N p − K eigenvectors are compiled into a noise matrixN, where the choice of K is contingent upon the SLC images. The image acquired with the MUSIC method is [37] Empirical parameter K is related to the actual number of phase centers in the complex forest structures. For the proposed simulation scenario, K is varied from 1 to 5 and the resulting images suggest that K = 2 is the most suitable choice. It is observed that as K increases, the reconstructed image tends to display ghost targets.

Amplitude and Phase Estimation (APES) Method
The amplitude and phase estimation (APES) method is implemented by optimizing a filter with finite impulse response (FIR) which is applied to the SLC imageḡ in (24) to obtain an output =h † ·ḡ (28) of which the power is given by whereR = E{ḡḡ † } is the covariance matrix of SLC images. By substituting (26) into (28), we obtain The vertical reflectivity profile is acquired by solving the following optimization problem If the filterh is designed to retrieve specific u at s while suppressing the other u components, namely, h † · [ā(s 1 )u 1 + · · · +ā(s −1 )u −1 +ā(s +1 )u +1 + · · · +ā(s L )u L ] → 0 then (31) is reduced to an optimization problem onh as [37] h A (s) = min where s is relabeled as s. The optimal solution of (34) is given by [37] h A (s) = where λ > 0 is an empirical parameter used to avoid possible singularity in matrix inversion [37]. Next, the optimal solution of u(s) is determined as then the reflectivity function is estimated as under the approximation that u A (s) is insensitive to n or n . By using (3), the dielectric constant is estimated as

Capon Method
The Capon algorithm has been used to detect the forest-height profile from a power spectrum [14]. It is implemented by using the same filter as derived in (27)-(33), The power spectrum is obtained by substituting (39) into (29) to give [14] Figure 5 shows the flowchart of the proposed TomoSAR procedure, including generation of the forest model, SAR tomography and validation of results.  Table 1 lists the parameters used in the simulations. The carrier frequency is 1.25 GHz, at which the backscattered signal power was claimed to be more uniform in the elevation (s) direction [12]. UAVs have been used to carry SAR, with a typical velocity of 0.1 m/s [32] to 5 m/s [42], moving at a height of 10 m [33] to 200 m [43]. In the simulations, we assume the UAVs fly at a height of H = 150 m, with a velocity of v p = 5 m/s, and the look angle is θ = 50 • . A range bandwidth of 500 MHz [32] or 600 MHz [42] has been reported in UAV-based SAR. The range bandwidth of 300 MHz is chosen to achieve the range resolution of 0.5 m. A pulse repetition frequency of 50 Hz [44] to 1 kHz [33] has been reported. By setting F a = 50 Hz, the synthetic aperture length is L a = 51.2 m and the azimuth resolution is ∆y = 0.54 m. A pulse duration of 10 µs [44] to 100 µs [33] has been used. By setting the pulse duration to T r = 1 µs, the corresponding range chirp rate is K r = 300 THz/s. The elevation resolution is given by ∆s = λR 0 /(2B n ). By specifying ∆s = 0.5 m, compatible with the range resolution and azimuth resolution, the baseline aperture is determined as B n = 57.22 m. The maximum canopy height in Arbaro is H c = 22 m [27], which corresponds to the maximum elevation of H c / sin θ = 28.7 m. The elevation ambiguity is given by s amb = λR 0 /(2d b ), which must be larger than H c (s) to avoid elevation ambiguity. By setting s amb = 28.7 m, the baseline spacing is d b = 0.49 m. Finally, the number of baselines is determined from B n and d b as N p = 115.

Parameters for TomoSAR Simulations
The effects of noise are also simulated under the signal-to-noise ratio (SNR) of 0 to −10 dB, where SNR = 10 log(U 2 s /U 2 n ), and U s and U n are the average amplitudes of signal and noise, respectively. Table 2 lists the parameters of a healthy taramack, a healthy quaking aspen and a dry (dead) tamarack, respectively, to be used in the simulations. The wood moisture and dielectric constant of a deadwood trunk are different from their healthy counterparts. Lower wood moisture implies a lower dielectric constant. A voxel model of the forest with a resolution of 0.25 m was generated by applying an open-source Binvox [45] on the airborne laser scanning (ALS) point clouds [30]. In this work, the Arbaro [27] is used to generate a 3D tree model, composed of trunk, branches and leaves. The data of the tree model are then mapped to point clouds to form a fine-resolution voxel model by using an open-source CloudCompare software v2.12.4 [46].
The dielectric constants of tree parts can be estimated with empirical formulas [7,[47][48][49][50]. The effective dielectric constant of each voxel can be determined by applying proper mixing formulas to the constituent matters within the voxel. A trunk with a diameter larger than the voxel size is modeled as a bunch of voxels filled with trunk matter. The dielectric constant of a voxel composed of leaves and air is estimated by using the Maxwell-Garnett mixing formula [51]. The dielectric constant of a voxel composed of air and branches is estimated by using the Maxwell-Garnett mixing formula for randomly oriented, highly prolate and rotation-symmetric ellipsoidal inclusions [52][53][54]. Figure 6a shows an example of point clouds of a tree generated with Arbaro [27], and Figure 6b shows the voxel model mapped from the point clouds with Binvox [45]. Each voxel is a cubic of size 0.5 m, the same as the range and azimuth resolutions specified in Table 1. The effective dielectric constant in each voxel is estimated by applying the Maxwell-Garnett mixing formula to the constituent leaves, branches and trunks. Figure 6c shows the LAD model of leaves generated with leafR [55], and Figure 6d shows the distribution of the dielectric constant in voxels of the tree model.

Simulations on Virtual Forest
Inspired by the tamarack-ignited wildfire that burned 68,637 acres in 2021 [56], Figure 7 shows the dielectric-constant profile of two virtual forests, composed of tamaracks and quaking aspen, generated with Arbaro [27]. There are 85 trees randomly dispersed in an area of 100 m × 100 m, with a tree spacing of 2.2-5 m [57]. Figure 7a shows a scenario with a deadwood tamarack at the center, close to the surrounding trees. Figure 7b shows another scenario with the deadwood at the center, but farther away from the surrounding trees. The proposed TomoSAR procedure is aimed to pinpoint this deadwood from the acquired dielectric-constant profile of the forest.
(a) (b) Figure 7. Dielectric-constant profile of virtual forest [27], (a) deadwood taramack at center, surrounded by healthy quaking aspen and taramack, (b) same as (a) except the deadwood taramack is farther away from the surrounding trees.
The fidelity of an acquired image a against the ground-truth image b can be evaluated with a structural similarity index (SSIM) defined as [58] where µ p and σ p are the mean and standard deviation, respectively, of image p, with p = a, b; σ ab is the covariance between images a and b; and c 1 and c 2 are stability constants. The SSIM index lies between 0 and 1; a higher index means the two images are more similar to each other. In the simulations, each image pixel is stored in 8 bits, corresponding to an integer between 0 and L = 255. The stability constants are set to c 1 = (0.01L) 2 and c 2 = (0.03L 2 ). The fidelity of an acquired image a against the ground-truth image b can also be evaluated with a root-mean-square error (RMSE) defined as [59] where a p and b p are the values of the pth pixel in images a and b, respectively, and P is the number of pixels in each image. Figure 8 shows the ground-truth dielectric-constant profile in Figure 7a at y = y c and the reconstructed images with the methods of CS, FB, MUSIC, APES and Capon, respectively. Around x = x c of the reconstructed images, the deadwood is slightly obscured by the surrounding healthy quaking aspen and taramack, which have a higher dielectric constant than the deadwood due to the difference in moisture content. Figure 8c,d show a clear background, while Figure 8b,e,f manifest artifact streaks in the background. The streaks in the acquired images have similar patterns but different magnitudes, which may be attributed to numerical errors in the covariance matrix. Figure 8b,c restore the tree shape and dielectric constant more accurately than the other three reconstructed images. The contrast of the dielectric constant in the canopy is less obvious in Figure 8b than in Figure 8c. The SSIM index of each image against the ground-truth image is listed in Table 3, which is consistent with the visual inspection. The ground is blurred and only partially revealed, possibly because the incident wave and backscattered wave are attenuated by the canopies. Since the deadwood is very close to the surrounding healthy trees in this scenario, the former is not well discernible in all the reconstructed images. Figure 9 shows the ground-truth dielectric-constant profile in Figure 7b at y = y c and the reconstructed images with the methods of CS, FB, MUSIC, APES and Capon, respectively. In this scenario, the deadwood is farther away from the surrounding healthy trees than in Figure 7a. In general, Figure 9b,c match the ground-truth image in Figure 9a better than the other three images. Figure 9b shows a higher contrast of the dielectric constant than the other reconstructed images, and Figure 9c shows the sharpest tree shape and the cleanest background compared with the others. Figure 9b,e,f manifest similar patterns of artifact streaks, which may be related to numerical errors in the covariance matrix of SLC images. Figure 9d displays some speckles in the background, which are not observed in its counterpart of Figure 8d. The deadwood part in Figure 9b,c matches well to that in the ground-truth image.
No calibration is needed to acquire the image in Figure 9b with the CS method, but some artifact streaks appear on the background. The image in Figure 9c acquired with the FB method vividly imitates the true forest and the background is clean.
The proposed approach can be used to reconstruct the 3D dielectric-constant profile of a target region at a fine spatial resolution by deploying a swarm of UAVs as platforms. The time decorrelation issue of multi-pass operations with airplanes or satellites is also avoided. Table 3. SSIM indices of acquired images.

Effects of Noise
The fidelity of the acquired TomoSAR images may be deteriorated by noises from the atmosphere and other sources, which is quantified against the backscattered signal in terms of the signal-to-noise ratio (SNR). Next, simulations on the two scenarios in Figure 7a,b are conducted at SNR = 0 dB, −5 dB and −10 dB, respectively. Figure 10 shows the TomoSAR images acquired on the scenario in Figure 7a under SNR = 0 dB. Compared with their counterparts in Figure 8, the images in Figure 10a with the CS method and Figure 10d with the APES method look blurrier, while the images in Figure 10b,c,e are barely affected, consistent with the SSIM indices listed in Table 3. The SSIM index suggests that the acquired image with the MUSIC method matches the best to the ground-truth image.  Figure 7a under SNR = −5 dB and SNR = −10 dB, respectively. The FB and the MUSIC methods outperform the other three methods in terms of forest configuration and background cleanliness. The artifact streaks in the images with the APES and Capon methods are strong enough to engulf the trees, especially under SNR = −10 dB. The performance of these two methods relies on the accurate estimation of the covariance matrix from the SLC images. Strong speckles appear in the images acquired with the CS method, but the forest configurations are still recognizable. Among these five methods, the MUSIC method acquires the cleanest background and the FB method best preserves the forest configuration. Figure 13 shows the TomoSAR images acquired for the scenario in Figure 7b under SNR = 0 dB. In this scenario, the deadwood is slightly separated from the surrounding trees. The image acquired with the CS method manifests the highest contrast of dielectric constants compared with the other methods. Compared with its counterpart in Figure 9b, background speckles grow stronger. The background in the images acquired with the FB and MUSIC methods are almost unaffected by noise. The forest configuration is also well preserved with these two methods. The deadwood in the center is much better reconstructed in the images as compared to their counterparts acquired for the scenario in Figure 7a. The artifact streaks grow strong in the images acquired with the APES and Capon methods, respectively, similar to their counterparts in Figure 10.  Table 3 lists the SSIM index of the images in Figures 8-15. In general, the FB method achieves the highest SSIM index under noise-free conditions, and the MUSIC method achieves the highest SSIM index under noisy conditions. Table 4 lists the RMSE indices of the images in Figures 8-15. As the noise level is raised, the Capon method and the APES method tend to generate more stripes in the reconstructed images; hence, their RMSE indices are relatively higher than the other methods. The FB and MUSIC methods yield relatively lower RMSE indices, and the reconstructed images better resemble the ground-truth images. In general, the MUSIC method or the CS method yields the lowest RMSE index across the whole area, while the FB method yields the lowest RMSE index within the square area.

Performance Comparison of Imaging Methods
In general, the FB method and the MUSIC method deliver better TomoSAR images of the simulated forest scenarios than the other three imaging methods. The compressivesensing (CS) method has the flexibility of reconstructing TomoSAR images with sparse baselines and does not require normalization, but it takes more computational loading than the other methods.
The hyperparameter K in the MUSIC method is sensitive to the number of phase centers, which can be empirically adjusted to enhance the acquired images under different noise conditions. Under noise-free conditions, the FB method reconstructs clear images and is less prone to developing ghost targets. However, it may not effectively identify low dielectric-constant areas containing deadwood. Under severe noise, the FB method may render a noisy background image. The FB method takes the lowest computational load among these methods.
The Capon method and the APES method follow the same process of deriving the optimal filters. Both methods tend to produce stripes along the elevation direction in the images. The APES method was originally tried to avoid normalization, but it did not work out as expected.
The Capon and the FB methods proved to be the most computationally efficient. The quality of reconstructed images with the FB method excels under noise-free conditions, and that with the MUSIC method excels under noisy conditions. The CS method does not require normalization of the reflectivity function as the other four methods do. Future work may be developed to reduce the computational load of the CS method and remove the constraint of normalization on the FB, Capon, APES and MUSIC methods.

Highlighted Contributions
In conventional TomoSAR imaging methods, the backscattered signals received along multiple tracks are first compressed and co-registered along a specific elevation line segment in the target volume by using an empirical point spread function. Then, an inverse algorithm is applied to reconstruct the image along the elevation line segment. However, the dielectric-constant profile in the target volume cannot be accurately estimated without calibration or normalization. In our TomoSAR procedure, the backscattered signals are related to the reflectivity function via the scattering theory, and the image at specific azimuth y and specific depth r is reconstructed via a focusing function ξ nk , without resorting to calibration or normalization. Five different inverse methods are used to acquire the reflectivity function along the specified elevation line segment. The acquired images are analyzed to validate the proposed procedure and compare the performance of these five inverse methods in TomoSAR imaging.
The novelties and contributions of this work are summarized as follows.

1.
A complete and rigorous TomoSAR imaging procedure is developed by integrating scattering theory and signal processing models to acquire a 3D dielectric-constant profile of a forest at decimeter resolution. A focusing function is derived to upgrade the point spread function adopted in conventional TomoSAR methods, without resorting to empirical assumption or calibration.

2.
The imaging task can be carried out in near real time with a swarm of UAVs to scout hot-spots of wildfire hazards.

3.
A vivid forest voxel model is developed by integrating the Arbaro tree generator, Maxwell-Garnett mixing formulas and leaf area index (LAD) for calculating the dielectric constant of a leaf, branch and trunk.

Conclusions
A complete and rigorous TomoSAR imaging procedure is developed by integrating scattering theory and signal processing models to acquire a dielectric-constant profile of a forest at decimeter resolution. The imaging task can be carried out with a swarm of UAVs to scout hot-spots of wildfire hazards in a forest. A vivid forest voxel model is developed by integrating the Arbaro tree generator, Maxwell-Garnett mixing formulas and leaf area index (LAD) for calculating the dielectric constant of a leaf, branch and trunk. Simulation results validate the efficacy of the proposed procedure. Five different inverse methods are used to acquire and compare the dielectric-constant profiles in two forest scenarios under three different signal-to-noise ratios.
The simulation results indicate that a deadwood close to the surrounding trees is difficult to pinpoint. Future efforts are needed to enhance the spatial resolution and the contrast of the reflectivity function in the acquired images. The five imaging methods are primarily used to compare their merits and shortcomings on the proposed scenario. Future efforts are needed to update or tailor these methods for better TomoSAR imaging on forests.