Quantitative Analysis of SAR Image Geometric Distortion and Its Application in Deformation Rate Fusion Mapping

Interferometry synthetic aperture radar (InSAR) technology has emerged as a powerful tool for the early identification and monitoring of geological hazards in the Loess Plateau. Since the SAR satellite adopts the range-azimuth two-dimensional side-looking imaging mode, the corresponding field area of the pixel will change with slope, aspect, and incident angle, especially in areas with severe changes in terrain and landforms, which causes serious area distortion, leading to abnormal deformation rate mapping. Considering the SAR imaging geometry, terrain slope and aspect conditions, this article proposes a quantitative analysis method for accurate geometric distortion (GD) based on pixel area and corresponding field area and constructs a weighted fusion method based on distortion values for InSAR deformation rate maps in overlapping regions of adjacent tracks. Taking Lingtai County, Gansu Province as the research area and sentinel-1 data as the research data, we utilize the SBAS method to process the sentinel-1 data and estimate the deformation rate, and then use the GD value as a weight factor to estimate the fusion rate of the overlapping area of adjacent tracks. The results demonstrate that the weighted fusion method based on distortion values can achieve a high-quality mosaic effect of deformation rate maps of adjacent tracks, effectively improve the GD influence in deformation rate mapping, and the identification of geological hazards.


Quantitative Analysis of SAR Image Geometric Distortion and Its Application in Deformation
Rate Fusion Mapping Zhengyang Wang , Xu Ma, Junhuan Peng , Mengyao Shi , Yun Peng, Yuhan Su , Zhaocheng Guo, and Wenwen Wang Abstract-Interferometry synthetic aperture radar (InSAR) technology has emerged as a powerful tool for the early identification and monitoring of geological hazards in the Loess Plateau.Since the SAR satellite adopts the range-azimuth two-dimensional side-looking imaging mode, the corresponding field area of the pixel will change with slope, aspect, and incident angle, especially in areas with severe changes in terrain and landforms, which causes serious area distortion, leading to abnormal deformation rate mapping.Considering the SAR imaging geometry, terrain slope and aspect conditions, this article proposes a quantitative analysis method for accurate geometric distortion (GD) based on pixel area and corresponding field area and constructs a weighted fusion method based on distortion values for InSAR deformation rate maps in overlapping regions of adjacent tracks.Taking Lingtai County, Gansu Province as the research area and sentinel-1 data as the research data, we utilize the SBAS method to process the sentinel-1 data and estimate the deformation rate, and then use the GD value as a weight factor to estimate the fusion rate of the overlapping area of adjacent tracks.The results demonstrate that the weighted fusion method based on distortion values can achieve a high-quality mosaic effect of deformation rate maps of adjacent tracks, effectively improve the GD influence in deformation rate mapping, and the identification of geological hazards.

I. INTRODUCTION
T HE Loess plateau region exhibits diverse landforms and significant variations in topography, along with numerous geological disasters [1].Landslides in the region are highly destructive and occur frequently, causing significant economic losses and casualties.Conventional monitoring methods are based on geodetic and geophysical equipment to collect surface deformation through single-point measurement or network deployment.Despite their high monitoring accuracy, these methods cannot offer macro-scale monitoring results [2], [3].Interferometry synthetic aperture radar (InSAR) is an active remote sensing technology using microwave telemetry, which can make up for the shortcomings of conventional monitoring methods such as limited spatial scale coverage, low monitoring frequency, long period, and sparse point distribution.
Over the past two decades, time series InSAR (TS-InSAR) technology has been continuously developed and applied to long-term geological disaster monitoring [4], [5].Previous studies have illustrated that TS-InSAR can detect landslide areas and monitor the long-term cumulative deformation [6], [7].With the increase of SAR acquisition, numerous studies have analyzed the applicability of SAR data to obtain different types of landslides [8], [9], [10].In addition, numerous studies have explored the usage of multiple data or methods for joint analysis to enhance the reliability of deformation results [11], [12], [13].
In regions characterized by intricate topography, the effectiveness of InSAR technology in identifying and monitoring landslides may decrease due to the special side-looking imaging mode.This mode introduces geometric distortion (GD), leading to unreliable or even undetectable results.Despite the maturity of InSAR technology in landslide identification and monitoring, its application still has limitations.Therefore, numerous studies evaluate SAR image visibility through SAR imaging parameters and terrain parameters to determine data applicability.Notti et al. [14] proposed the R-Index method to identify the GD area by defining the ratio of the radar slant distance to the ground distance.Except for shadow area, the R-index effectively evaluates the impact of various terrains on the suitability of SAR images [15], [16].Dai et al. [17] quantitatively analyzed the correlation between line of sight (LOS) displacement and slope displacement.They generated a sensitivity map of LOS displacement using InSAR measurement to reveal the relationship among LOS displacement, actual slope displacement, and GD.To reveal the applicability of sentinel-1 data along the Qinghai Tibet Railway, they quantitatively analyzed the GD along the railway line to determine the suitability of ascending and descending track data [18].Under different geometric conditions, Chen et al. [19] summarized the influence of GD when the landslide faced the satellite and faced away from the satellite and obtained the observation effect of InSAR on slopes with different slope angles and aspect angles.In complex terrain, certain studies utilized slope, aspect, and incident angles to identify local incident angles on slop surfaces for identifying layover and shadow areas.This identification method is significant for SAR data selection [20], [21].These methods for assessing GD are all based on satellite LOS measurement and terrain parameters, without considering the changes of the entire pixel on different slopes [22], [23], [24].Changes in the slope angle and aspect angle of the ground result in variations between the actual length and image length in both the range and azimuth directions.Hence, there is a distinction between the image area and the corresponding field area.
This article proposes a method to improve the recognition of GD by introducing azimuth distortion, which can extend one-dimensional (1-D) range distortion to 2-D area distortion to analyze the degree of GD quantitatively.The method comprehensively considers the GD in both the range and azimuth directions and calculates the ratio of pixel area to the corresponding field area to obtain the distortion value.This value is then used to analyze the applicability of SAR data quantitatively.In addition, a weighted fusion method is constructed based on distortion values to fuse deformation rate maps in overlapping areas of adjacent tracks effectively and mitigate the impact of GD on deformation rate mapping.
The rest of this article is organized as follows: Section I begins by presenting the research difficulties and current situation through literature analysis.Section II introduces a novel method for quantitative analysis of GD and a fusion method for deformation data.Section III is the experimental chapter, which describes the experimental area and methods and provides the experimental results.Section IV analyzes and discusses the results of the quantitative analysis for GD and the fusion results of deformation data, as well as the validity of fusion results in identifying landslides.Finally, Section V summarizes the innovative content and experimental verification analysis presented in this article.

A. Azimuth Angle Deviation and Field Resolution Calculation
Since the SAR system employs the 2-D side-looking imaging mode, the corresponding field area of the pixel will change with slope angle, aspect angle, and incident angle, which causes serious area GD, including foreshortening, layover, and shadow, in generated SAR image [25].
Changes in the slope angle and aspect angle cause differences between the actual lengths of the ground and the corresponding lengths in the image in both the range and azimuth directions.To quantify the GD, we can then calculate the area difference based on the length differences at the corresponding position.
In the ascending track, the satellite flight azimuth is α h , and the LOS azimuth is α h − 270 • (see Fig. 1).The aspect angle of slope is β 0 , and the azimuth deviation β between the LOS azimuth and the aspect angle can be derived by According to the geometric relationship of SAR imaging, the single pixel in the image represents the approximate rectangular ABCD on the corresponding slope in the field [see Fig. 2].The lengths of the azimuth direction and the range direction of one Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
pixel on the slope are represented by l AC and l CD , respectively.After calculating the slope angle α and the azimuth deviation β, the slope angle α AC of l AC and the slope angle α CD . of l CD can be derived by When the range and azimuth resolutions of the SAR image are known, the actual lengths of the range and azimuth direction on the slope can be given by where ΔR and ΔZ are the range and azimuth resolutions, respectively.l CD and l AC are the actual lengths of the range and azimuth directions on the slope, respectively.

B. Quantitative Analysis of Geometric Distortion
To quantitatively analyze the GD, we first calculate the area of the pixel ΔR • ΔZ and the corresponding field area l CD • l AC .To accurately calculate GD, it is crucial to consider identifying shadow areas.Merely calculating distortion value based on the ratio of pixel area to corresponding field area is insufficient for this purpose.In order to address this limitation, a judgment parameter B can be incorporated into the calculation.By setting B = −1 when α CD > 90 • − θ on the slope facing away from the satellite.Then, the GD value σ is given by ( 6) as σ increases within the range of 0 ≤ σ ≤ 1, the difference between the pixel area and the corresponding field area also increases, resulting in more noticeable GD.The layover area is recognized when σ > 1, and the shadow area is recognized when σ < 0. Furthermore, the SAR satellite follows parallel tracks with similar azimuth angles but exhibits significant variations in incident angles.When the slope and aspect angles remain unchanged, the GD decreases with the increase of the incident angle.

C. Deformation Rate Fusion
The quality of data acquisition can be affected by the GD caused by the SAR side-looking imaging mode.In the same area, the measurement results of different tracks may be inconsistent or even contradictory.Therefore, the impact of GD must be considered when fusing data.In the InSAR data processing process, we also consider the benchmark deviation d O of the deformation rate between two tracks.After processing the InSAR data of two tracks within the study area using the SBAS method, we need to unify the spatial reference benchmark.To achieve this, we utilize geocoding and Kriging methods, which ensure the deformation rate maps are within the same geographic coordinate system [26].Finally, one track is selected as the master track, and the deformation rate of the slave track is projected onto the master track to unify the LOS deformation [27], [28].The prefusion deformation rate of the slave track can be calculated by where d P and d S represent the prefusion deformation rate and the original deformation rate of the slave track.θ M and θ S are the incident angles of the master track and slave track.d O represents the deviation between the deformation rate maps of the two tracks.
Considering the layover and shadow areas in one orbit, which may result in inaccurate or missing deformation information, restoring the accurate deformation by utilizing the data from another track becomes essential.In order to accomplish this, appropriate weights must be set before data fusion.We utilize the condition of σ < 0 to detect the shadow area and designate the rest area as a layover area based on σ > 1.Thus, the weight assignment method is given by ⎧ ⎨ ⎩ P M = 0 and P S = 1, 0 ≤ σ S ≤ 1 < σ M P M = P S = 0, σ S > 1 or σ M < 0 P M = 1 and Without layover or shadow areas (0 ≤ σ i ≤ 1, i = M, S), weight determination is directly based on the GD value.The weight assignment method is given by where P M and P S are the weights of the master track and the slave track.
After configuring the weights of each region, the fusion of the deformation rate is given by where d fusion represents the fusion deformation rate and d M represents the deformation rate of the master track.

A. Study Area and Datasets
Lingtai County is situated in the gully area of the Loess Plateau in the southeastern part of Gansu Province, between the Jinghe River and the Weihe River.The altitude of the region is between 900 and 1520 m, with the terrain showing a downward trend from northwest to southeast.There are three types of geography: remnant plains; hills; and tundra.The region has a typical continental climate, with concentrated rainfall from July to September.The characteristics of loess collapsibility  are evident here, and rainwater infiltration further increases the occurrence rate of geological disasters.The unique topography and climatic conditions of Lingtai County provide a natural condition necessary for the occurrence of landslides, collapses, and other geological disasters.
We collected sentinel-1 data and high spatial resolution remote sensing data from two tracks to monitor deformation in the study area [29].The local incidence angles of PATH84 FRAME110 and PATH157 FRAME107 are 33.8°and43.8°i n the study area, respectively.These different incident angles enable data acquisition under distinct observation geometries, leading to distinct observation results.The range sampling spacing and azimuth sampling spacing are 2.33 and 13.99 m, respectively.In order to ensure the consistency of the deformation rate, we selected the images of the two tracks from October 3, 2020, to May 26, 2023, totaling 113 images [see Table I].The study area and coverage of datasets are shown in Fig. 3.

B. Results of Geometric Distortion
We introduce the azimuth factor into GD analysis, increasing the traditional 1-D range GD to 2-D area distortion.The degree of GD is evaluated by comparing the difference between the pixel area and the corresponding field area.To determine the applicability of sentinel-1 data in the study area, we processed ALOS WORLD 3D 30-m DEM to acquire terrain information (the slope angle and aspect angle).According to the slope angle, aspect angle, and SAR satellite parameters, we utilized the GD quantitative analysis method to generate GD maps for Lingtai County in PATH84 FRAME110 and PATH157 FRAME107 [see Fig. 5].
Different colors are used in GD maps to distinguish the degree of GD.The distortion value determines the color of each point on the map.The red color indicates a higher distortion value, while the green color indicates a lower distortion value.Moreover, the layover area is represented by blue, and the shadow area is represented by black.

C. Results of Deformation Rate Fusion
Before time-series deformation analysis, we preprocessed SAR dataset covering the study area.Initially, one image was selected as the main image for registration and resampling of the remaining images in the frame.To suppress the influence of decorrelation noise, the multilooking operation (range × azimuth: 10 × 2) was applied.
After preprocessing, an appropriate spatiotemporal threshold was set to obtain a small baseline subset.During differential interference processing, the DEM was utilized for the removal of the terrain phase.In addition, we selected the same reference point in the overlapping area to eliminate the benchmark deviation between deformation rate maps.Maintaining consistency in reference point selection between two tracks is crucial, as inconsistent choices can result in deviations between deformation rate maps.Finally, the SBAS method was employed to process the differential interferograms and obtain the deformation rate points.
Due to variations in observation geometry, the coherent points identified in two parallel tracks may differ when monitoring the same area.To achieve spatial consistency and subsequent deformation rate fusion, we geocoded the coherence points into the WGS84 coordinate system and performed Kriging interpolation to obtain the two track deformation rate maps with evenly distributed coherence points [30], [31].Once the complete deformation rate maps were obtained, the fusion process could be achieved.In contrast to conventional data fusion methods that use equal weight, we incorporated GD values as prior weights to avoid inaccurate deformation data caused by GD.Fig. 4 illustrates the comprehensive flowchart of data process.
After correcting for the deviation, the deformation rate maps of PATH84 FRAME110 and PATH157 FRAME107 exhibit a high consistency.The fusion map clearly indicates regions with significant deformation [see Fig. 6].

A. Overall Geometric Distortion Analysis
The SAR satellite maintains a consistent flight direction in both tracks, with a local incident angle of 33.8°for PATH84 FRAME110 and 43.8°for PATH157 FRAME107 in the study area.Consequently, the overall GD is minor in PATH157 FRAME107 than in PATH84 FRAME110 due to the influence of the incident angle.We conduct statistical analysis on the GD Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 4. Flowchart of data processing.The data processing is divided into four parts.The first part is SAR dataset preprocessing, the second part is SAR data temporal processing, the third part is GD quantitative analysis, and the fourth part is deformation data fusion.II].The GD values of PATH84 FRAME110 mainly concentrate in the range of 0-0.8, amounting to 96.8%.In addition, GD is severe in the range of 0.8 to 1.0, amounting to 3.1%.In PATH157 FRAME107, most of the GD values are within the range of 0-0.8, amounting to 99.8%.Furthermore, only 0.2% of the GD values fall within the range of 0.8-1.0.The GD value greater than 1 indicates the presence of layover area.For PATH84 FRAME110, only 0.1% of the areas are identified as layover.In contrast, no layover areas are found in PATH157 FRAME107.A GD value less than 0 indicates the shadow area, but no shadow areas are found on either track.

TABLE II GEOMETRIC DISTORTION PROPORTION STATISTICS values of different intervals in two regions [see Table
Table II illustrates that the slope facing the satellite exhibits significant GD in the two tracks.Conversely, the slope facing away from the satellite exhibits relatively minor GD.Moreover, the degree of GD is influenced by the change in the incident angle of the satellite.As the incident angle increases, the degree of GD decreases.Since the incident angle of PATH157 FRAME107 is larger than that of PATH84 FRAME110, the overall GD of PATH157 FRAME107 is significantly reduced [see Fig. 5 ].
In the study area, the joint monitoring using parallel track dataset can reduce the influence of GD, which is in response to the problem that only ascending track sentinel-1 data can be acquired in most areas of the Loess plateau.According to the statistical analysis of the GD of the two tracks in Table II, the serious foreshortening areas can be reduced from 3.1% to 0.2%, and the layover areas can be reduced from 0.1% to 0% through joint monitoring.The results show that combining parallel track sentinel-1 data enhances the applicability of observation results.

B. Geometric Distortion of Different Slopes
The SBAS method is utilized to process the SAR dataset and obtain deformation points.These points are subsequently geocoded into the WGS84 coordinate system to determine their distribution.We investigate the influence of GD on the spatial distribution of deformation points by overlying them on the analysis maps for GD.PATH84 FRAME110 exhibits significant overall GD and a relatively sparse distribution of deformation points compared to PATH157 FRAME107.On the other hand, PATH157 FRAME107 exhibits a dense and uniform distribution of deformation points overall [see Figs.7 and 8].
We analyze the relationship between terrain features and GD by selecting five slopes (S1-S5) with varying slope angles and aspect angles in the white dashed lines in Figs.7 and 8.The subgraphs on the right are optical images showing the slope area with blue circle, slope direction with yellow arrow, and deformation points with pink dots.At the same time, we calculate the actual GD values by measuring the average intervals in both range and azimuth directions of deformation points on the slopes using Google earth.To assess the accuracy of the proposed GD quantitative analysis method, we compare the quantitative analysis results of GD with the actual GD results, along with the R-index results (see Table III) [14].
The optical image shows that S1 is flat ground, and the deformation points on both tracks are uniformly distributed.All three GD calculation methods achieve the same results.The other four slopes have different slope angles and aspect angles.We compare the GD results of four slopes obtained by three calculation methods and find that the GD values obtained by the proposed method are more consistent with the actual distortion values, with a maximum deviation of 0.05 and an    average deviation of −0.007.As a result, this GD quantitative analysis method can accurately assess the applicability of SAR data in complex terrains.

C. Deformation Rate Fusion Analysis
After unifying the spatiotemporal benchmark and reference point, we fused the two original deformation rate maps and the final deformation rate fusion map showed good consistency in spatial distribution with the two original deformation rate maps.To evaluate the reliability of the fusion result, we compare the deformation rates of two tracks in overlapping regions.The statistical distribution map is drawn to analyze the deformation rate differences of 140 494 samples in the overlapping area for further quantitative analysis [see Fig. 9].Before fusion, the mean value of the deformation rate differences is −0.8 mm/year [see Fig. 9(a)].After that, the mean value of the deformation rate differences changes from −0.8 mm/year to −0.1 mm/year [see   It is worth noting that the mean value after fusion approaches zero, indicating the elimination of systematic bias in the deformation rate fusion map.In addition, the standard deviation of the latter (3.3 mm/year) is smaller than the former (5.6 mm/year), suggesting an enhancement in accuracy.

D. Validity of Deformation Rate Fusion Map for Landslide Identification
The long-term deformation rate reveals the spatial and evolutionary features of ground deformation.From the deformation rate fusion map, the regions with obvious deformation are clearly visible.After delineating the deformation areas, we assess the potential landslide characteristics by considering the morphological characteristics and sliding signs observed in the optical remote sensing image.Finally, we have identified four potentially dangerous slopes [see Fig. 10].In the optical images, these four potential landslides are characterized by clear sliding boundaries, distinct trailing edges, and accompanied by local collapse [see Fig. 10(a-4)-(d-4)].
The observation effects of PATH84 FRAME110 and PATH157 FRAME107 vary due to different observation geometries in the same area, leading to inconsistent GD.To assess the impact of GD, we analyzed the deformation of these four slopes across two tracks.In PATH84 FRAME110, both slope a and slope d exceed their landslide boundaries in terms of deformation area [see Fig. 10(a-1) and (d-1)], while in PATH157 FRAME107, the deformation area of each slope aligns more closely with the landslide boundary [see Fig. 10(a-2) and (d-2)].However, the fusion map demonstrates greater consistency in both the deformation range and landslide boundary [see Fig. 10(a-3) and (d-3)].
The slope b and slope c show smaller deformation area and deformation degree in PATH84 FRAME110 [see Fig. 10(b-1) and (c-1)], but are more significant in PATH157 FRAME107 [see Fig. 10(b-2) and (c-2)].After the fusion processing, these two slopes are more prominent in the fusion map [see Fig. 10(b-3) and (c-3)].

V. CONCLUSION
Aiming at the problem of GD caused by SAR satellite sidelooking imaging mode, this article proposes a quantitative analysis method for GD based on pixel area and corresponding surface area.Using this method, we identify layover and shadow areas and quantitatively describe the extent of GD in the study region.To obtain deformation points from SAR datasets, we employ the SBAS method.In addition, we use the Kriging method to generate two complete deformation rate maps.To mitigate the influence of GD on deformation results, we fuse the two deformation rate maps by assigning weights based on the GD values.
By analyzing the deformation rate fusion map, we can effectively identify regions exhibiting substantial deformation.Coupled with the examination of optical images, we detected four potential landslides.The main conclusions are as follows.
1) The combination of GD quantitative analysis map and deformation point distribution indicates that the proposed method can accurately identify distortion areas.According to the results of statistical analysis, the degree of GD is relatively small when the slope is facing away from the satellite and relatively large when facing the satellite.At the same time, varying incident angles can alter the overall GD of the study area.2) Based on the weighted fusion method to process the In-SAR deformation rate in the overlapping area of adjacent tracks in the study area, the average value of the deformation rate differences changes from -0.8 to -0.1 mm/year, eliminating the systematic error between different tracks.The standard deviation decreased from 5.6 to 3.3 mm/year, showing an improvement in precision.
3) The weighted fusion method, utilizing GD values, effectively mitigates the influence of GD on deformation results.Using the fused deformation rate map for hidden hazard identification focuses on the issue of indistinct deformation or inaccurate deformation area caused by GD, leading to a more precise determination of the distribution and area of potential dangers.

Fig. 1 .Fig. 2 .
Fig. 1.Azimuth angle deviation between LOS direction and slope aspect.The azimuth angle deviation β is illustrated.α h represents the satellite flight azimuth.

Fig. 5 .
Fig. 5. GD analysis results of (a) PATH84 FRAME110 and (b) PATH157 FRAME107.The subgraphs depict the GD of various orbits in the same region on an enlarged scale.

Fig. 7 .
Fig. 7. GD and pixel distribution of SAR images in PATH84 FRAME110.S1-S5 represent five distinct areas with varying slope angles and aspect angles.The two subgraphs on the right are optical images, in which the yellow arrows indicate the slope direction, blue circles depict the slope area, and pink dots demonstrate the distribution of pixels on the slope.

Fig. 8 .
Fig. 8. GD and pixel distribution of SAR images in PATH157 FRAM107.S1-S5 represent five distinct areas with varying slope angles and aspect angles.The two subgraphs on the right are optical images, in which the yellow arrows indicate the slope direction, blue circles depict the slope area, and pink dots demonstrate the distribution of pixels on the slope.

Fig. 9 .
Fig. 9. Histogram of deformation rate difference in the overlapping area.(a) Before fusion.(b) After fusion.

Fig. 10 .
Fig. 10.Identification of potential landsides by deformation fusion result.Four potential landslides are indicated by (a)-(d), with the red line denoting their extent and the yellow arrow indicating their direction.The deformation of the potential landslides is represented by a series of subgraphs: (a-1)-(d-1) for PATH84 FRAME110, (a-2)-(d-2) for PATH157 FRAME107 and (a-3)-(d-3) for the fusion data.In addition, (a-4)-(d-4) are the optical images of the corresponding slopes.

Fig. 9 (
Fig.9(b)].It is worth noting that the mean value after fusion approaches zero, indicating the elimination of systematic bias in the deformation rate fusion map.In addition, the standard deviation of the latter (3.3 mm/year) is smaller than the former (5.6 mm/year), suggesting an enhancement in accuracy.

TABLE I MAIN
PARAMETERS OF THE SAR DATA USED IN THIS ARTICLE

TABLE III ANALYSIS
OF GEOMETRIC DISTORTION IN TYPICAL REGIONS