{'fraction_non_alphanumeric': 0.025707692807696057, 'fraction_numerical': 0.009717572881796614, 'mean_word_length': 4.771900206340274, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Terrain radiometric calibration of airborne UAVSAR for forested area

In the field of biomass estimation, terrain radiometric calibration of airborne polarimetric SAR data for forested areas is an urgent problem. Illuminated area correction of σ -naught could not completely remove terrain features. Inspired by Small and Shimada, this paper tested gamma-naught on one mountainous forested area using airborne Uninhabited Aerial Vehicle Synthetic Aperture Radar data and found it could remove most terrain features. However, a systematic increasing trend from far range to near range is found in airborne SAR cases. This paper made an attempt to use the relationship between distance to SAR sensor and γ-naught to calibrate γ -naught. Two quantitative evaluation methods are proposed. Experimental results demonstrate that variation of γ -naught can be constrained to a limited extent from near range to far range. Since this method is based on ground range images, it avoids complicated orthorectification.


Introduction
Forests are one of the most important components of terrain ecosystems and play a vital role in maintaining the ecosystem and carbon dioxide balance. Using the SAR backscattering coefficient to estimate the biomass and carbon storage is one popular research field of remote sensing (1)(2)(3)(4).
NASA's Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR) is a fully polarimetric L-band instrument designed and built for acquiring airborne repeat-track interferometry SAR data. Its applications include monitoring ground deformations caused by earthquakes, volcano dynamics, ice dynamics, local sea ice dynamics, time-varying evaporation and hydraulic properties of soils, and aboveground biomass (5). This radar was designed to be part of an Uninhabited Aerial Vehicle, but is currently being tested on a NASA Gulfstream III (6). Its swath width in range direction is about 16 km and the nominal flying height is about 13.8 km. The accuracy of flying a predefined path is extremely high with the assistance of high precision real-time global positioning system (GPS) and a sensor-controlled flight management system (7).
UAVSAR's received signal is inevitably influenced by many factors in the scanned area including terrain, roughness, soil moisture, vegetation type, and so on (8)(9)(10). As a side-looking system, it is most sensitive to terrain (9). In applications that focus on nonterrain features such as land use, biomass, and soil moisture, these terrain-induced distortions must be calibrated to allow researchers to identify robust relationship between backscatter coefficients and these features (10,11).
During the processing of polarimetric SAR data, a flat ground or an ellipsoid earth is usually assumed (12)(13)(14). However, study sites such as forest stands are sometimes located in hilly or even mountainous regions thus violating the flat earth assumption. van Zyl et al. (15) found that ignoring terrain will cause more than one decibel error in the illuminated area correction and a large look angle estimation error in antenna pattern correction. Thus, it has been suggested that digital elevation model (DEM) should be used to calibrate the terrain influence on SAR data (12).

A review of current terrain calibration techniques
At present, researchers mainly implement three types of terrain correction: illuminated area correction, antenna pattern correction, and incidence angle correction (12,15,16).

Illuminated area correction
Illuminated area correction is replacing areas computed under no terrain assumptions with real illuminated areas. Its methods are relatively mature and can be classified into two categories. One is orthorectification-based area correction method and the other is ground range imagebased area correction. The first category, mainly proposed by , utilizes the corresponding relationship between slant range pixels and ground range pixels to compute the real illuminated area. In spite of the strictness of these methods, orthorectification is computationally intensive and requires highprecision orbit data and fine-resolution DEM. As to second category methods, several projection factors were raised (8,(19)(20)(21)(22) to obtain the real illuminated area.

Antenna pattern correction
Antenna pattern correction replaces the elevation angle under no terrain assumptions with real values from a DEM to get more precise antenna gain (15). Since antenna pattern correction is already performed in UAVSAR product, we will not address this issue here.

Incidence angle correction
The two corrections described in Sections 1.1 and 1.2 cannot totally remove the influence of terrain. The change of local incidence angle across the whole scene causes large fluctuations in the backscattering coefficient. Some researchers tried to correct the influence of incidence angle, and advanced a number of semi-empirical and empirical models (11,16,(23)(24)(25)(26) that can be considered part of scattering process correction. Although these model-based correction methods can partly solve problems in certain cases, to what degree they are correct is still uncertain.

Calibration for airborne and spaceborne systems
Most of current calibration methods were tested on spaceborne SAR data while only a few papers were on airborne SAR data. As is well known, the high orbit of a spaceborne system makes the estimation of local incidence angle relatively reliable. The variation of incidence angle across the whole scene in spaceborne systems is just a few degrees and much less than that of airborne systems (15), hence the error associated with illuminated area correction will not vary much across the swath width. In airborne case such errors vary considerably (15). As pointed out by van Zyl (22), error in antenna pattern correction without terrain can be neglected but for airborne systems this error may be fairly large. In airborne SAR, the backscatter coefficient contrast for foreslope and backslope is more obvious. Therefore, methods suitable for spaceborne SAR may not work for airborne SAR. Unfortunately, most of the airborne SAR papers were devoted to retired AIRSAR system. Only one paper by Zhang et al. (11) used the newest UAVSAR data for a flat area, so the assumption was made that the incidence angle equals the look angle. It remains unknown if the conclusions in that paper can be used for hilly or mountainous areas as well. For these reasons, it is critical to find a calibration method suitable for airborne SAR or at least UAVSAR.
The aims of this paper are as follows: (1) To propose quantitative methods to evaluate to what degree the image is terrain-flattened. (2) To compare the illuminated area corrected sigmanaught and gamma-naught and determine which one is better for terrain calibration purposes for forested areas.
(3) To explore if there is a systematic trend of backscatter coefficients for different channels along range direction and propose an empirical method to correct this trend.

Data
All original data was acquired in PolSAR mode by UAVSAR in August 2009 and were processed in 2010. The processed products were downloaded from the Alaska Satellite Facility DAAC (https://ursa.asfdaac. alaska.edu). Nevertheless, only part of the data was used because the proposed method is purely based on ground range image. The data used in this study includes: SRTM DEM in the WGS84 coordinate system (HGT file); illuminated area divided complex cross products of HH, HV, VH, and VV projected to SRTM DEM (GRD file); and a text file containing metadata (annotation file).

Preprocessing
(1) Based on the metadata information provided in the annotation file, SRTM DEM were projected to the UTM coordinate system with 5 m resolution using cubic interpolation; (2) Original GRD files corresponding to S hh S hh ⁄ , S hv S hv ⁄ , and S vv S vv ⁄ were converted to sigma-naught (often denoted as r 0 ) which is widely used for the characterization of distributed targets using equation (1). r 0 hh , r 0 hv , and r 0 vv files were also projected to UTM coordinate system with the same resolution as DEM using nearest-neighbor interpolation. Gaussian filtering using a 7 Â 7 window was employed to remove the speckles from the sigma-naught images but retaining the main features. The resulting images are referred to as r 0 orig in this paper.

Forest pixel extraction
In the field of radar remote sensing, forest scattering is widely considered to be mainly volume scattering in contrast with building characterized by dihedral scattering and water, barren land characterized by single scattering (35). Among many available polarimetric decomposition methods, Freeman-Durden's three component method (FDD) is one of the most famous methods and can estimate the powers of volume scattering, double-bounce scattering, and surface scattering (36). Due to the fact that the decomposition results are powers corresponding to each scattering type, the absolute powers have large dependency upon local incidence angle. Hence, the proportion of volume scattering power to the total power R vol is used instead. Forest pixels are extracted using the following steps: (1) R vol is computed using FDD in 7 Â 7 window; (2) r 0 hv is classified into more than five classes using unsupervised classification. This is performed to keep the boundary precise and avoid boundary positioning blur induced by 7 Â 7 window processing; (3) connected regions are acquired using four direct neighborhood for all classes that possibly contain forests; (4) the mean R vol for each region is used to distinguish forest regions from non-forest regions. An empirical threshold of 0.8 is set and region with mean R vol higher than this threshold is considered a possible forest region; (5) small regions containing too little pixels are excluded from possible forest regions.All pixels left in the forest regions are regarded as forest pixels.

Incidence angle computation
The computation of local incidence angle involves the estimation of SAR sensor position S and also the local normal vector L ! at pixel P. Simple linear interpolation was utilized based on sensor position at the center of the flight track and heading angle to estimate the SAR position S corresponding to each image pixels. Since the ability of UAVSAR flying a predefined straight path is excellent, the interpolated track is precise enough for this study. L ! is computed from SRTM DEM using four direct neighborhoods of P. Local incidence angle h incid was calculated using where acos is inverse of cosine function.

Illuminated area correction for sigma-naught
The illuminated area correction factor proposed by Ulander (20) was employed due to its use across a broad range of applications (19,37,38). An illuminated area calibration factor f r was developed for applying a DEMbased area calibration and removing the imprecise calibration performed in the UAVSAR GRD product. f r was computed by where u is complementary to the smallest angle between the surface normal and the image plane, cos / is the calibration factor given by Ulander's method (20), r e is the local radius of curvature of earth, h is global average sensor altitude, h t is the global average terrain height, and h l is the corrected angle between look vector and height axis (7). Illuminated area-corrected sigma-naught r 0 correct was computed using

γ-naught computation
Small (19) argued that r 0 correct fails to account for radar backscatter changes in regions with significant topographic variation and that γ-naught normalization is terrain-flattened. Shimada and Ohtaki (33) also hold that using γ-naught for diffuse reflection can reduce the incidence angle dependency. Why γ-naught is better than σ-naught in eliminating terrain features still remains unknown. One possible reason is that the decimeter-level wavelength of the L band wave determines that extremely rough land covers like forests can be considered as a diffuse scatterer. For the lambert scatterer, an ideal isotropic diffuse scatterer, it is assumed that the radiant intensity observed by it is directly proportional to the cosine of h incid (denoted as cos incid ), so if radiant intensity is divided by cos incid , the result would remove the dependency upon h incid . σ-naught and γ-naught also have the following relationship.
Therefore, if σ-naught cannot remove influence of terrain or incidence angle, an improvement may be obtained by using gamma-naught. However, γ-naught has only been successfully demonstrated using spaceborne SAR data (19,25). In order to test its applicability on ground range images from UAVSAR, this paper calculated gamma-naught. Here, a factor f c given in Equation (6) is directly employed for sake of its simplicity and high efficiency.
Pixels with negative cos incid are assigned a fixed small factor 0.001 since they are located in the shadow area. This small factor will not significantly change the original low value. f c is multiplied with r 0 correct to get c 0

Gamma-naught empirical calibration
In σ-naught images, the most prominent pixel value change in SAR image scene is in range direction. Generally, the near range has much larger σ-naught than the far range even after illuminated area correction. Many researchers relate it with incidence angle and experiments have shown that there is a negative relationship between σ-naught and incidence angle (11,16,(23)(24)(25)(26).
That relationship is stronger for smoother surfaces than for vegetated surfaces (39). γ-naught calibration may not be able to completely remove the brightness change in range direction. The ultimate goal of calibration in forested area is to make the same forest type share nearly the same backscatter coefficient no matter they are located in near range, middle range, or far range. To achieve this goal, an empirical model is needed to find and calibrate this trend.
(1) Grid distribution The survey area was equally divided into n groups in range direction and m nonoverlap segments in azimuth direction.
(2) Computation of mean backscatter coefficient of forest pixels in each group The mean backscatter coefficient values of forest pixels in each grid were calculated. Then each group's mean value was computed. In ideal case, forests with the same type, age, disturbance, and soil moisture must be employed to find the range direction trend. However, that is impossible to achieve since detailed information is unavailable. However, the mean value of one group with millions of samples could remove the influence of different forest parameters to a large degree. For this reason, the mean backscatter coefficient of forest pixels in each group is comparable.
(3) Computation of mean distance of each group to SAR sensor Scatter plots were made between the mean distances of each group to SAR sensor and backscatter coefficients. In this way, backscatter coefficients having a systematic change in range direction can be clearly seen in scatter plots. If a trend does exist, the mathematical relationship between these two variables can be estimated and a model fitted to obtain an empirical function.
(4) Correction based on an empirical model To apply an empirical model to c 0 , Bayer's method (23), originally based on local 1n angle was slightly modified. Suppose the correction function is F(dis), where dis is the distance from pixel to sensor, the reference distance is dis ref , then the correction factor f dis is Distance-corrected c 0 dis is This method is another form of choosing the reference incidence angle for calibration. Researchers often relate σ-naught or γ-naught to incidence angle instead of distance to sensor. In fact, when the terrain is flat, there is an approximate relationship between distance to sensor and incidence angle as follows where H F is the plane height, H T is local terrain height at pixel P, dis is distance from P to S. From Equation (10), h incid has positive relationship with dis. If necessary, we could roughly translate the distance to incidence angle and express Equation (8) as a function of incidence angle. Remember that in a mountainous area, dis is much more stable than h incid . Therefore, dis is more suitable as a reference.

Evaluation method
The most popular method to evaluate calibration results is examining the reduction in variance (8,23,26). But one single indicator cannot reflect some aspects of local dependency on terrain in the corrected backscattering coefficient. The other widely used method is using scatter plots of previous backscattering coefficients against calibrated coefficients or the incidence angle against the calibrated coefficient along a profile (8,9,18,19,25,33). But generally, the profile number and sample number per profile is limited. One visual and two quantitative evaluation methods were used to assess the effect of different terrain calibration methods.
(1) Visual check A visual evaluation of the calibration results checks whether the obvious terrain-induced radiometric distortions still exist especially in the foreslope and backslope of mountainous areas.
(2) Linear regression analysis Appropriate calibration will nearly remove all the dependence of backscatter coefficient on cos incid with a low R 2 and low absolute value slope in a linear regression. Each segment is divided into many groups at an interval of 60 m in the range direction. Mean backscattering coefficient and cos incid values of these groups were used in simple linear regression as dependent variable and independent variable separately. One benefit of this method is that it makes full use of the dramatic change of cos incid along range direction hence leading to more robust regression.
(3) ANOVA for backscattering coefficients The second quantitative method is using ANOVA to test if there are significant differences between mean backscattering coefficient values of different groups in the range direction. When the sample number and group number are fixed, a lower F value or higher p value indicates a smaller difference between the mean values of different groups. Consequently, different calibration results could be validly compared within the same channel. ANOVA was performed on two levels: the segment level and whole data level. In a range direction, each segment is divided equally into three groups at the segment level or five groups for the whole data level.

Results
Data were processed and analyzed using Matlab, IDL, and ArcMap. Part of the ANOVA results including the mean of all segments' p values for the segment level and the F value in whole data level are listed in Table 1.
Mean R 2 along with the slopes from all segments' linear regression (some outlier are excluded in regression) are also given.
The HV channel for the mountainous Hubbard brook forest is used to illustrate the calibration results, considering that the HV channel is the most sensitive channel for forest biomass (2). The flying direction is from northeast to southwest. The forest mask is shown in Figure 1. Figure 2(a) and (b) clearly shows that height fluctuates sharply in Hubbard brook forest, especially in the west and southwest part of the forest as well as in a small part in the northeast. The calculated cos incid in Figure 2 (c) exhibits the influence of terrain on local incidence angle. The incidence angle at foreslope (slopes facing Radar) is obviously smaller than that at backslope which has a distinct correspondence in r 0 orig of HV as shown in Figure 2(d). In Table 1, large R 2 and slope shows that r 0 orig strongly depends on cos incid for all three channels. However, the relationship seems stronger in HH and HV than VV. Visual check of the scatter plots of cos incid vs. r 0 orig reveal the existence of heteroscedasticity which makes it nearly impossible to model the relationship with a simple linear model.
The image of f r displayed in Figure 2(e) indicates that r 0 orig is enhanced in backslope and decreased in foreslope to get r 0 correct . Therefore, the contrast in Figure 2 (f) is less obvious. Nonetheless, as shown in Table 1, the large R 2 and slope still show a high dependence of r 0 correct on incidence angle, while a large F and small p shows the unevenness of r 0 correct in all three channels. Heteroscedasticity still exists. Hence, it could be stated that only performing illuminated area correction for sigma-naught is not enough. This conclusion is consistent with the argument of Small (19).
It seems f c is a much stronger compensator than f r as depicted in Figure 2(g). Note that the low cos incid area in lower left part of Figure 2(c) receives a much larger f c than image center. The derived c 0 is shown in Figure 2 (h). It seems all obvious terrain features have been flattened in this image and the grey value seems fairly uniform. The distinction between foreslope and backslope is nearly invisible. The existence of a small area with rough features in the lower left part may be explained by the fact the small cos incid leads to invalid r 0 orig . According to Table 1, the R 2 and slope have become much smaller, as is the F value at the level of the whole data while the average p at the segment level is many times larger. Therefore it is reasonable to say that the dependence of c 0 on cos incid is much less than the previous backscatter coefficients while achieving more uniform results, however, as shown in Table 1, in every channel it is possible to have an average R 2 not significantly close to zero. The near range has larger c 0 , so further calibration is needed to remove this range direction trend.
Due to the improved performance of c 0 in eliminating terrain features as compared with r 0 correct , it was selected as the method to determine the systematic change in the backscatter coefficient for range direction. Seven groups at the whole data level were employed to find the relationship between c 0 and dis. The reason seven groups were used was because seven is an odd number and suitable for ordinary model fitting. The scatter plots for c 0 and distances are given in Figure 3(a).
In the polynomial fitting process, to avoid the inconvenience caused by using large distance measurements in matrix operations, the distance of the central group is used as reference. Here the reference distance is 19.811 km. and third-order fitting is used. Table 2 gives fitting results for different channels. Figure 3(a) shows results for HV. The R 2 values exceed 0.97 and the small RMSE demonstrate that the fitting results were excellent. HH and HV fit better than VV. Based on the fitting models and correction method given in previous section, the final correction functions F(dis) are computed and given in Figure 3(b). It can be seen that when the distance to sensor is small, i.e. in near range, the correction factors for three channels are nearly the same, especially HH and HV, however, in the far range, their correct factors diverse and HH is the highest, HV is the second, and VV is the lowest. According to mathematic analysis, VV has two stationary points while HH and HV have none, so the curve of HH and HV only show a monotone increase trend but VV has one local maximum at the middle range and a local minimum at far range. Figure 2(j) shows the result after distance correction. The lower right part of the image will be compensated and upper left part will be suppressed as shown in Figure 2(i). Compared with c 0 the image of c 0 dis in Figure 2(j) seems more uniform, while the R 2 histogram centered near zero and most R 2 values were below 0.3. This time, about half of the slope is negative and significantly different from previous positive slopes. As depicted in Table 1, the decreasing F value at the whole data level and increasing mean P at the segment level were observed in all channels. That demonstrates the usefulness of distance correction. Figure 4 compares box plots of different groups of c 0 and c 0 dis for different channels. In (a), (c), and (e), the mean values of five groups are significantly different. However, in (b), (d), and (f), the significant overlap of the middle two quartiles of each group as well as the valid data range evidently indicate that there is no distinct difference in c 0 dis from the near range to the far range. Group 5 has the widest valid data range and indicates that the change of c 0 dis in the near range is larger than other groups. At the same time, some extremely high c 0 dis was found in a few portion of mountainous areas in the near range. A manual check shows that all have incidence angle close to zero, however, not all areas with incidence angle close to zero have abnormally high c 0 dis . To demonstrate that the proposed method is applicable to different forest types and different areas, the r 0 orig and c 0 dis of two sites located in North Carolina (NC) and Maine (ME) of USA are illustrated in Figure 5. Clearly, c 0 dis works well, looks fairly uniform, and retains features.

Discussion
In the creation of a forest mask, FDD is employed. A well-known drawback of this method is that the volume scattering power is overestimated because all HV power is considered to be from volume scattering which often does not hold true (40). Hence, Freeman-Durden decomposition may be insufficient when describing forests. Other decompositions, like Yamaguchi decomposition and van Zyl decomposition, need to be tested.
Although f r can compensate the r 0 in the backslope and suppress the r 0 in the foreslope, the small range of its value means that f r cannot completely make the backslope and foreslope share the same grey value, since the contrast in r 0 for two kinds of slope is usually big. This confirms the argument in Small (19). But c 0 can remove terrain features and reduce dependence on incidence angle. It can be stated with confidence that γ-naught is at least suitable for diffuse scatterer and can roughly account for the influence of local incidence angle.
However, c 0 still shows a systematic increasing from far range to near range. This may be because forests are not strictly isometric and have stronger c 0 at smaller incidence angles. Small and Shimada's method works on spaceborne data and no distance correction is needed. The distances of all pixels in the same scene do not differ too much. However, in airborne cases, the largest distance maybe two times that of the smallest one, so the influence of distance to sensor must be taken into account. Distance-correction functions derived from fitting a third-order model could minimize the variation of c 0 and calibrate most sites, including challenging terrain.
Even in the best calibration results, the p value at the segment level is not so large. In one segment, there usually exists certain change in soil moisture, forest type, age, and species; making the backscatter coefficient not totally the same. These local features also determine the existence of a positive and negative slope. In addition, the p value only shows the similarity of the mean backscatter coefficient. However, the overlap of valid ranges of each group after distance correction still proves the usefulness of this calibration method.
Although F value of distance-corrected γ-naught at the whole data level was much smaller than that of previous backscattering coefficients, its corresponding p is still small. That maybe contributed to the incorrect forest mask and unevenness of forest closure, species, and ages in inner composition of each group. Some low canopy-closure forests have lower c 0 than that of high canopy-closure forest. The authors believe if exactly same forest type is used in mask, the correction function would be more valid and evaluating result would be more visually appealing.
A few high backscatter coefficient forest pixels in the near range still exist after distance correction. The reason still remains unknown perhaps due to inaccurate estimation of look angle in the antenna pattern correction and use of a higher gain factor. Another explanation is some low-closure forest may reflect more energy back to SAR sensor when incidence angle is close to zero, while the ground serves as a good specular scatterer.
Three channels, HV, HH, and VV seem not to be able to be calibrated using the same distance correction function. HH and HV behave with a similar, increasing trend in the correction function. Nonetheless, the VV correction function is highest in middle range and furthest range. A similar phenomenon was reported by Freeman (41). Although no many researchers are interested in the VV channel, it is still worthwhile to investigate why it behaves that way.
It is hoped a higher resolution DEM could be employed by UAVSAR data processing teams instead of coarse SRTM DEM. Interpolation from 30 to 5 m certainly would not add more information to DEM. Inaccurate DEM causes the loss of many detailed features especially in hilly or mountainous areas and would lead to the wrong estimation of look angles and imprecise orthocorrection and low positioning accuracy  of features (15,28). With the development of newest measuring technology like airborne LiDAR and InSAR, the availability of high-resolution and high-accuracy DEM has greatly improved, so it is natural to use the latest and most accurate DEM in polarimetric data processing. The National Elevation Dataset provide nationally seamless elevation data at resolution of 1 arc-second (about 30 m) and 1/3 arc-second (about 10 m), and in limited areas at 1/9 arc-second (about 3 m) (42). The given information for the UAVSAR flight path is too limited. More flight path information would be desirable.

Conclusion
The terrain radiometric calibration of airborne polarimetric SAR for forested areas is an urgent problem. Illuminated area correction of r 0 cannot remove terrain features. However, c 0 can remove most terrain features for diffuse scatterers like forests but still shows a systematic increase from far range to near range. Using the relationship between distance to SAR sensor and c 0 results in correction functions, so the variation of c 0 is constrained in a quite limited range from near range to far range. Hence, distance-corrected γ-naught is better than pure γ-naught as used by Small and Shimada. Since this method is based on ground range images, it avoids complicated orthorectification. Future work includes: determining if distance-corrected γ-naught can remain constant when data are acquired using different viewing geometries. If it remains constant, shadow area have a potential to be filled by data from another viewing geometry (43). It seems c 0 is at least suitable for the research of diffuse scatterers like forests. Whether it could be used for specular scatterer remains to be investigated.