Classification of Land Cover in Complex Terrain Using Gaofen-3 SAR Ascending and Descending Orbit Data

Abstract

Synthetic aperture radar (SAR) image is an effective remote sensing data source for geographic surveys. However, accurate land cover mapping based on SAR image in areas of complex terrain has become a challenge due to serious geometric distortions and the inadequate separation ability of dual-polarization data. To address these issues, a new land cover mapping framework which is suitable for complex terrain is proposed based on Gaofen-3 data of ascending and descending orbits. Firstly, the geometric distortion area is determined according to the local incident angle, based on analysis of the SAR imaging mechanism, and the correct polarization information of the opposite track is used to compensate for the geometric distortion area, including layovers and shadows. Then, the dual orbital polarization characteristics (DOPC) and dual polarization radar vegetation index (DpRVI) of dual-pol SAR data are extracted, and the optimal feature combination is found by means of Jeffries–Matusita (J-M) distance analysis. Finally, the deep learning method 2D convolutional neural network (2D-CNN) is applied to classify the compensated images. The proposed method was applied to a mountainous region of the Danjiangkou ecological protection area in China. The accuracy and reliability of the method were experimentally compared using the uncompensated images and the images without DpRVI. Quantitative evaluation revealed that the proposed method achieved better performance in complex terrain areas, with an overall accuracy (OA) score of 0.93, and a Kappa coefficient score of 0.92. Compared with the uncompensated image, OA increased by 5% and Kappa increased by 6%. Compared with the images without DpRVI, OA increased by 4% and Kappa increased by 5%. In summary, the results demonstrate the importance of ascending and descending orbit data to compensate geometric distortion and reveal the effectiveness of optimal feature combination including DpRVI. Its simple and effective polarization information compensation capability can broaden the promising application prospects of SAR images.

1. Introduction

Accurate and timely land-cover mapping can provide reliable basis data for regional development and environmental management [1]. Many land-cover products have been developed with traditional optical images. However, high-quality traditional optical images cannot always be obtained, due to weather restrictions. Fortunately, with the development of active microwave remote sensing technology, synthetic aperture radar (SAR) is becoming a new important data source for land-cover mapping relying on its all-weather imaging capability. However, the backscattering energy of ground objects is influenced not only by the surface roughness but also by the surrounding ground objects, their inherent topography, and the incident angle of the SAR sensor. These uncertain factors limit the classification accuracy of SAR images, especially in complex terrain areas. Complex terrain refers to areas with significant variations in terrain height, shape, and detail, which can make land-cover identification and area quantification challenging. Its research has aroused widespread interest [2,3].

Recently, the effectiveness of SAR images in land-cover classification has been demonstrated [4,5]. However, most of these methods focus on the whole region without considering the geometric distortion of the complex terrain, which leads to errors in land-cover classification. Mountain areas cover about 24% of the earth’s land area. The side-view and range-imaging characteristics of SAR sensors cause geometric distortion in the images [6,7], such as layover and shadow, which may lead to scattering characteristics within the images inconsistent with the actual scene. This issue is especially acute in rolling mountains or other complicated terrain. For example, the layover phenomenon can cause the backscattering coefficient to reach the building area, and it is impossible to perceive the real types of ground objects in the layover area [8]. Moreover, the intensity and scattering angle of SAR responses are sensitive to the land surface conditions, especially in regions of complex terrain with different land-surface conditions and diversity of vegetation types. Exploring more effective feature combinations of SAR data in complex terrain regions is an important task to further improve land-cover classification accuracy [9,10].

At present, the research into land-cover classification in complex terrain areas concentrates mainly on two aspects: First, the research method based on geometric distortion regions mainly focuses on detection [11,12,13]. Second, some researchers are trying to extract efficient features to further increase the land-object discriminative capacity of SAR data.

In view of the degradation of image quality caused by the existing deformation of complex terrain, some researchers have focused on the use of interferogram analysis and geometric modeling to complete the detection of geometric distortion areas [14,15,16]. Ren et al. (2013) used correlation coefficient and interference amplitude to detect layover and shadow areas [17]. Zhang et al. (2019) determined the layover and shadow areas by using a geometric model and a morphological method of radar satellite imaging [18]. Their study showed that both the interferogram detection method and the geometric model can detect the geometric distortion region to a certain extent. However, the polarization information of geometric distortion areas has not been supplemented or corrected, and cannot properly describe the types of ground objects, nor can it be directly applied to classification. Hence, the compensation of the geometric distortion area is essential for enhancing the accuracy of land-cover classification based on SAR data. Mahdavi et al. (2019) found that the cross-track image can collect most of the information lost by the sensor in the line of sight of the opposite orbit [19]. Borlaf-Mena et al. (2021) classified the stay and shadow areas of front occlusion to avoid the influence of geometric distortion, but the real information contained was still unknown [20]. However, only a few researchers have tried to use multi-track data to compensate for the effects of geometric distortion [21,22].

An effective combination of features is one of the key factors of land-cover classification using SAR images [23,24,25,26]. Up to now, backscattering coefficient and polarization characteristics have been the main characteristics of SAR land-cover classification. Several researchers have investigated feature combinations using full-polarimetric SAR data. Guo et al. (2018) extracted twenty-two feature parameters from Radarsat-2 full-pol SAR images, and then selected the best feature-by-feature selection method to identify rice [27]. Sayedain et al. (2020) thought that better land-cover classification results could be obtained by combining NDVI with full-pol features [28]. Shen et al. (2020) used GF-3 to extract polarization decomposition features and backscattering coefficient for water extraction [29]. Compared with the full-pol modes, dual-pol modes are easier to obtain, with larger swath width and lower data volumes. However, feature separability inthe dual-polarization data of different land objects is lower and more susceptible to the impact of complex terrain. Therefore, an appropriate feature combination is important to ensure the accuracy of land-cover classification based on dual-pol modes with limited polarimetric information. Backscatter intensity from polarized channels has been widely used to identify different types of land cover [30,31,32]. Other derived information such as the H-α polarization decomposition parameters, grey level co-occurrence matrix (GLCM) [33,34,35], and co-pol ratio has been explored for classification and valuable classification results have been obtained [36]. Recently, DpRVI for dual-pol SAR data has been proposed, which has been successfully utilized in crop-growth monitoring and studies of soil moisture [37,38]. However, the effect of DpRVI and its feature combinations for land-cover classification in complex terrain areas still needs to be investigated. Moreover, these assessments have mainly been explored using VV-VH SAR data. Further experiments using the HH-HV mode are required, as the land-object responses of horizontally polarized transmitted waves and vertically transmitted waves could be different.

Based on the analysis of the above-mentioned research status, the main problems faced during the classification of ground objects in complex terrain areas can be summarized as follows: on the one hand, the current research mainly focuses on the detection of geometric distortion areas, and the real classification results of ground objects in geometric distortion areas have not been explored. On the other hand, the effect of DpRVI in land-cover classification has not been investigated. In view of the above two problems, this research aims to develop a classification method for complex terrain based on the dual-pol mode (HH-HV) SAR data from the flight path of the Gaofen-3 satellite, to improve the accuracy of land-cover classification. Firstly, by analyzing the geometric model of radar satellite imaging, the local incident angle is obtained from the radar wave side angle and slope angle. Conditional judgment is made to identify the layover and shadow area, and compensated by the data of another orbit. Then, the relevant features are extracted for each ground object. Based on the dual orbit polarization data, the DpRVI is introduced to characterize the vegetation-cover type and Jeffries–Matusita (J-M) distance analysis is performed to select the optimized feature combination. Finally, the deep learning method 2D-CNN with excellent classification ability is applied to classify the compensated images. The main contributions of this paper are as follows: designing a framework that compensates geometric distortion area by using the data of the opposite trajectory, determining the specific types of ground objects on the basis of detected geometric distortion areas. It broadens the application range of SAR remote sensing technology and has good application prospects in the field of land-cover classification. Furthermore, it demonstrates the application of DpRVI for land-cover classification using dual-pol HH-HV SAR data, especially the ability to identify vegetation types.

2. Materials and Methods

2.1. Study Area

The study area is located in Danjiangkou reservoir area (111°24′35″E~111°52′16″E, 32°31′6″N~32°50′21″N), which is the core water source area of South-to-North Water Transfer Project. It is an important ecological protection area in China, shown in Figure 1. Accurate land-cover investigation is of great significance to further protect the ecological environment.

According to official data for 2020 released by the Bureau of Land and Resources Utilization, buildings, woodlands, farmland, and water are the main types of ground elements. The slope-angle map on the right side of Figure 1b shows that the terrain fluctuates greatly. In the west and south of the study area, there are continuous mountainous areas which cause serious geometric distortion in SAR images. Major land-cover types in the eastern and northern parts of the study area are towns and farmland. The study area is characterized by its complex terrain and relatively rich land-cover types.

2.2. Data and Preprocessing

2.2.1. Gaofen-3 Data and Preprocessing

Gaofen-3 is a special remote sensing satellite of China’s High-Resolution Earth Observation System Project. It is the first C-band multi-polarization SAR imaging satellite in China and has the characteristics of multiple imaging modes, large imaging width (10~650 km), multiple resolutions (1~100 m), and multiple polarization modes. This study utilizes two images from the Gaofen-3 fine stripe (FSII) SLC product of ascending and descending orbit, with resolution of 10 m in dual polarization (HH-HV) mode. The detailed data parameters are shown in

Table 1. Parameters for Gaofen-3.

Gaofen-3 Parameter	Master Image	Slave Image
Product	SLC	SLC
Image mode	Fine stripe mode II	Fine stripe mode II
Incidence angle	31.597256	31.783311
Polarization	HH, HV	HH, HV
Pixel interval	10 × 10 m	10 × 10 m
Band	C	C
Pass direction	Ascending	Descending
Date	9 July 2017	9 July 2017

.

2.2.2. Auxiliary Data

In this study, shuttle radar topography mission (SRTM) DEM with a resolution of 30 m was selected as auxiliary data. Geometrical model construction of terrain correction and SAR imaging was completed by DEM.

The study area is mostly mountainous and it is difficult to obtain sufficient actual land-cover information from field trips. Hence, corresponding optical images from Google Earth with 1 m resolution from the same period were selected as reference annotation sources to provide samples of land cover for visual interpretation. These images were also used to evaluate qualitatively the relative accuracy of land-cover classification. The distribution of land-cover samples is shown in Figure 2. Specifically, 148 parcels were selected as samples, including 33 buildings, 58 farmland, 11 water areas, and 46 forest areas. In this study, 60% of the samples were randomly selected as training samples, and the remaining 40% were used as test samples.

2.3. Geometric Distortion Region Analysis

After standard image pre-processing, such as geocoding and topographic correction, geographic distortion remains present in a large number of SAR images, which limits their classification accuracy for complex terrain. The premise of dealing with the geometric distortion area is to identify the specific location. In this paper, the spatial geometric relationship of the synthetic aperture radar’s imaging is utilized.

The imaging method of the SAR is shown in Figure 3. Related to the projection order of the ground-target points, distant ground objects are projected on the image later. The distance to the target is directly proportional to the order of receiving information.

$$B^{\prime }-C^{\prime }=K\left(R_B-R_C\right)$$

(1)

where B’ and C’ are distance coordinates, R_B and R_C are the oblique distance of the target point, and K represents the coefficient of the actual distance and the image distance. That difference between ${\beta }_2$ and ${\beta }_3$ is very small, where $\lambda$ is the oblique irradiation angle $\lambda =\frac{\pi }{2}-{\beta }_3+\alpha$, which can be calculated as follows:

$$B^{\prime }-C^{\prime }=\left(B-C\right)\cos \left(\frac{\pi }{2}-{\beta }_3+\alpha \right)$$

(2)

It can be seen from Equation (2) that the larger the incident angle, the smaller the scale coefficient K. The phenomenon in Figure 3 shows the short-range compression characteristics of SAR imaging.

As the slope angle increases, the relief of the terrain becomes more pronounced. As shown in Figure 4a, the oblique distance from the top of the slope decreases continuously, and the local incident angle decreases continuously, until the oblique distance B2 from the top of the slope to the radar satellite equals the oblique distance A from the bottom of the slope to the satellite. At this point, the positions of the bottom of the slope and the top of the slope in the image coincide. At this time, the local incidence angle is 0°. Starting from this critical point, as the slope angle increases, the top of slope B3 is recorded before the bottom of the slope A in the SAR images, thus resulting in layover areas.

When the satellite illuminates the back of the mountain slope, the tilt angle of the back is too large, and the radar beam is blocked by the top of the mountain. Thus, it cannot illuminate the target area on the back of the mountain, resulting in no echo signal in this area, which is a shadow area. As shown in Figure 4b, when the back slope angle is ${\alpha }_2$, it is just at the critical point of the shaded area. At this time, $\lambda =2\pi -\alpha$, then $\theta$ can be expressed as:

$$\theta =\beta -2\pi -\alpha =\beta +\alpha$$

(3)

when the local incident angle $\theta$ is greater than 90$°$, the radar signal is blocked by the mountain, which is a shadow area.

In general, the geometry of SAR images is determined by the local incidence angle of the target point. When the local incidence angle is less than 0°, the radar incidence angle is smaller than the slope angle. The imaging position of the top of the slope in the SAR image is in front of the bottom of the slope, which is the layover area. When the local incidence angle is larger than 90°, the radar signal is occluded and no data are generated; this is the shaded region. Mountainous areas in complex terrain fluctuate greatly, and the problem of layover and shadow is more prominent.

2.4. Land-Cover Classification Method for Complex Terrain Area

The flowchart of our proposed method is shown in Figure 5. After data preprocessing, the process of land-cover classification can be divided into three steps. Firstly, the geometric model of SAR satellite imaging is analyzed and the geometrically distorted region is identified based on the range of the local incidence angle. Secondly, the geometric distortion area is compensated by using images of different tracks, in order to reduce the adverse effects of layover and shadow on the classification accuracy of ground objects. Thirdly, dual-orbit polarization characteristics (DOPC) and DpRVI are extracted from the compensated image, and J-M distance analysis is carried out to verify the sample separability of different feature combinations. Finally, the proposed feature combination is used for classification by 2D-CNN, which is used to generate the land-cover classification model. Its accuracy is evaluated and compared with DOPC + DpRVI feature combination based on uncompensated images and DOPC without DpRVI based on compensated images.

2.4.1. Geometric Distortion Area Detection

On the basis of the above analysis, the geometric conditions of layover area and shadow area were analyzed according to the spatial geometric relationship of SAR image imaging time. Satellite orbit data and DEM data of the study area were used for constructing the SAR imaging model, mainly including spatial coordinate transformation and morphological filtering. Firstly, the position parameters of Gaofen-3 were calculated to fit the orbit data, and the coordinate system and the three-dimensional coordinates of the study area obtained by DEM were unified to construct the imaging model. Secondly, the local incident angle was judged conditionally, and the rough geometric distortion area determined. Finally, morphological filtering was performed to reduce the influence of noise and obtain the fine geometry. The detailed implementation steps are shown in Figure 6.

In detail, satellite orbit parameters and the trinomial fitting method were used to fit the orbit of the satellite. The orbit parameter file contains the position vector (X_S, Y_S, Z_S) and velocity vector (X_V, Y_V, Z_V) of many satellite groups, as shown in

Table 2. Parameters for track data.

Parameter	Formula
XS	$X_S={\displaystyle {\displaystyle \sum }_{i=0}^n}{a}_i{t}^i$
YS	$Y_S={\displaystyle {\displaystyle \sum }_{i=0}^n}{b}_i{t}^i$
ZS	$Z_s={\displaystyle {\displaystyle \sum }_{i=0}^n}{c}_i{t}^i$
XV	$X_V={\displaystyle {\displaystyle \sum }_{i=1}^n}i{A}_{i-1}{t}^{i-1}$
YV	$Y_V={\displaystyle {\displaystyle \sum }_{i=1}^n}i{B}_{i-1}{t}^{i-1}$
ZV	$Z_V={\displaystyle {\displaystyle \sum }_{i=1}^n}i{C}_{i-1}{t}^{i-1}$

. Where $a_i,t^i,b_i,c_i,A_i,B_i,C_i$ are the coefficients to be determined. The coefficients were obtained by fitting the actual data with the trinomial ballistic equation.

Three-dimensional coordinates of the research area were obtained from DEM. DEM is an ordered numerical array that describes the surface morphology and spatial distribution of information within the whole simulated scene, reflecting the actual surface morphology. The satellite orbit coordinate system is the WGS-84 coordinate system while the coordinate system adopted by DEM is a geodetic coordinate system. Therefore, the satellite orbit coordinate (XYZ) was converted into the geodetic latitude and longitude (BLH) coordinate system to establish the unified coordinate system required for the experiment. The specific coordinate conversion expression is:

$$\mathrm{B}=\arctan \left(\frac{\mathrm{z}+{\mathrm{N}\ast \mathrm{e}}^2\sin \mathrm{B}}{\sqrt{{\mathrm{x}}^2+{\mathrm{y}}^2}}\right)$$

(4)

$$\mathrm{L}=\arctan \left(\frac{\mathrm{y}}{\mathrm{x}}\right)$$

(5)

$$\mathrm{H}=\left(\frac{\sqrt{{\mathrm{x}}^2+{\mathrm{y}}^2}}{\cos \mathrm{B}}\right)-\mathrm{N}$$

(6)

where N is the radius of a unitary circle, that is:

$$\mathrm{N}=\frac{\mathrm{a}}{\sqrt{1-{\mathrm{e}}^2{\sin }^2\mathrm{B}}}$$

(7)

where a = 6,378,137 m.

According to the actual needs, the geodetic coordinate system of known target points was converted into a Gaussian plane coordinate. Thus, the geometric relationship between radar satellite and irradiation target can be solved.

On the basis of the unified coordinate system, a geometric model was constructed to obtain the instantaneous local incident angle at this time, as shown in Figure 7.

The local incidence angle $\theta$ is expressed as:

$$\theta =\beta -\alpha$$

(8)

where $\beta$ is the radar wave side angle; $\alpha$ is the inclination angle of the target point. When $\theta <0$°, the target area produces top–bottom displacement, which is layover area. When $\theta >90$°, the target area produces signal loss, which is shadow area.

In the currently detected geometric distortion area, morphological methods were used to remove burr noise. Image corrosion and a dilation algorithm are the basis of morphological image processing, which extracts the geometric features of useful objects from images according to the logical relationship between pixels. The method used in this paper involves first expansion and then corrosion.

• Expansion operation:

  $$A⨁B=\left\{x,y|{\left(B\right)}_{xy}\cap A\ne \varnothing \right\}$$

  (9)

The rough geometric distortion region is dataset A, where 1 is the geometric distortion region and 0 is the normal region. This formula shows that the dataset A is expanded by the structure B, and the origin of the structure element B is translated to the image pixel (x, y) position. If the intersection of B and A at the image unit (x, y) is not null, the unit (x, y) corresponding to the output image is assigned as 1, otherwise it is assigned as 0.

• Corrosion operation:

  $$A\ominus B=\left\{x,y|{\left(B\right)}_{xy}⊑A\right\}$$

  (10)

This formula shows that structure B is used to corrode A. When the origin of B is translated to the pixel position (x, y) of image A, if B is at (x, y) and completely contained in the layover area of image A, then the pixel (x, y) corresponding to the output image is identified as a geometric distortion region with value 1, otherwise it is identified as a normal region with value 0.

2.4.2. Geometric Distortion Area Compensation

Detection can identify the position of the geometric distortion area, but it cannot obtain effective polarization information. Consequently, the characteristics of scatter error in layover and shadow areas still lead to incorrect classification. The images obtained using SAR sensors from different tracks in a specific area might contain different polarization information and have different geometric distortion areas because of different incident angles. That is, the geometric distorted areas in the current trajectory image may not be deformed in the opposite orbit. Therefore, most of the lost information in the line of sight of the remote sensor can be collected in the opposite orbit.

An image from any orbit can be selected as the primary image, and the image of the opposite orbit is selected as the secondary image. In this paper, the descending orbit image was set as the main image and the ascending orbit image was set as the secondary image. The geometric distortion area in the main image is replaced by the pixel value at the corresponding position in the secondary image. The image-fusion formula of ascending orbit and descending orbit is as follows:

{[M(i,j)∈A_S]∨[M(i,j)∈A_L]}→R(i,j) = S(i,j)
{[M(i,j)∈A_N]} →R(i,j) = M(i,j)

(11)

where M(i,j) is the pixel value in the main image. A_L, A_S, and A_N are the layover area, shadow area, and normal area in the main image, respectively. S(i,j) represents the corresponding image pixel value of the secondary image. R(i,j) is the pixel value in the compensation image based on the data of ascending and descending orbits.

2.4.3. Feature Extraction

The features used in the current study mainly include dual-orbit polarization characteristics (DOPC) and DpRVI. DOPC is a dual-pol feature combination including backscattering coefficient (σ₀), total scattering power (SPAN), differential intensity (DI), power ratio (PR), and H-α polarization decomposition parameters. These features are selected for their benefits to distinguish different types of surface roughness, such as water, farmland, buildings, and woodlands for land-cover classification.

The framework of feature extraction used in the current study is shown in Figure 8. After the compensation, the polarization matrix of the image was converted to generate polarization matrix C2; that is, each pixel in the image is represented by a covariance matrix, as shown in Formula (12):

$$\mathrm{Z}=\left|S_{HH}{S}_{HV}\right|\to C2=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}^{*}& C_{22}\end{array}\right]$$

(12)

where * stands for conjugate transpose, $C_{11}$ and $C_{22}$ are real numbers. Z is the target vector, S_HHS_HV defines the scattering component, and V/H represents the vertical and horizontal transmit/receive polarization.

Parameters such as eigenvalues and eigenvectors are used to represent the polarization information of ground objects in SAR images. The C2 polarization matrix in Formula (12) can also be written as:

$$C2={{\displaystyle \sum }}_{i=1}^2{\lambda }_i{e}_i{e}_i^{*}$$

(13)

where ${\lambda }_i$ and $e_i$ are the eigenvalues and eigenvectors of the C2 matrix, respectively.

SPAN is an effective parameter to represent the intensity information of scattering between pixels. The $SPAN$ value of urban areas is higher, while that of water areas is lower. It is expressed as follows:

$$SPAN=\frac{1}{2}\left({\left|S_{HH}\right|}^2+{\left|S_{HV}\right|}^2\right)$$

(14)

where |S_HH|² and |S_HV|² are the intensities of HH and HV bands, dB (decibel) is the logarithm of a given number x to base 10 or log₁₀ x.

In addition, the difference between the co-polarization and cross-polarization intensity is also reflected in DI and PR.

$$DI=\frac{1}{2}\left({\left|S_{HH}\right|}^2-{\left|S_{HV}\right|}^2\right)$$

(15)

$$PR=\frac{{\left|S_{HV}\right|}^2\left(dB\right)}{{\left|S_{HH}\right|}^2\left(dB\right)}$$

(16)

H-α polarization decomposition was applied to analyze the polarization matrix, which is mainly based on the eigenvalues and eigenvectors of C2. The physical quantities of eigenvalues and eigenvectors are defined as follows:

• Entropy indicates the randomness of target scattering:

  $$H={{\displaystyle \sum }}_{i=1}^2-P_{i}lo{g}_3{P}_i , 0\le H\le 1$$

  (17)

  where $P_i=\frac{{\lambda }_i}{{{\displaystyle \sum }}_{j=1}^2{\lambda }_j},i=1,2$, ${\lambda }_i$ is the eigenvalues of the polarization matrix C2;
• Alpha is the average scattering mechanism from surface scattering to volume scattering and then to dihedral angle scattering:

  $$\alpha =P_1{\alpha }_1+P_2{\alpha }_2$$

  (18)

  where ${\alpha }_i$ is the internal degree of freedom of the scatterer, with a value ranging from 0° to 90°;
• Anisotropy indicates the degree of anisotropy of target scattering:

  $$A=\frac{{\lambda }_2-{\lambda }_1}{{\lambda }_2+{\lambda }_1}$$

  (19)

A dual polarization radar vegetation index (DpRVI) for monitoring crop conditions was previously developed based on dul-pol SAR feature vector decomposition [38]. In this current study, we aimed to utilize DpRVI extracted from HH-HV SAR data to enhance the distinguishability between farmland and forest. Its expression is:

$$DpRVI=1-DoP\ast P$$

(20)

where $P=\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2}$ is the normalization parameter of eigenvalue, ${\lambda }_1$ represents the even and odd scattering intensity, ${\lambda }_2$ represents multiple scattering intensity. When the target is dominated by odd or even scattering, ${\lambda }_1\ge {\lambda }_2$ exists. DoP in Equation (21) is used to quantify the relative intensity between different scattering mechanisms of polarization. For dual polarization data, its expression is:

$$DoP=\sqrt{1-\frac{4C2}{Tr{\left(C2\right)}^2} }=\left|\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1+{\lambda }_2}\right|$$

(21)

where Tr () is the trace of the matrix. According to the definition of the degree of polarization, the relative intensities of different scattering mechanisms are quantified. When DoP = 1, it is considered that the target is in a completely polarized state that is dominated by a single scatter mechanism. When DoP = 0, it is considered that the target is in a completely depolarization state dominated by multiple co-existing scattering mechanisms.

2.4.4. J-M Distance

After extracting the above features from the SAR images, J-M distance was used to analyze and verify the separability of these features. This is a useful tool for evaluating the effectiveness of features. The J-M distance varies from 0 and 2. When the value is greater than 1.9, it indicates that the features are well separable. When the value is greater than 1.8 and less than 1.9, the separability of features is generally possible. When the value is less than 1, these features are inseparable. This is shown in Equation (22):

$$D_{ij}={\{{{\displaystyle \int }}_x^{}{[\sqrt{P\left(\frac{x}{w_i}\right)}-\sqrt{P\left(\frac{x}{w_j}\right)}]}^{2}dx\}}^{0.5}$$

(22)

where $D_{ij}$ represents the separability of category i and category j, and $P\left(\frac{x}{w_i}\right)$ and $P\left(\frac{x}{w_j}\right)$ indicate the conditional probability density of feature x under category i and category j, respectively.

2.4.5. 2D-CNN Classifier

In order to effectively exploit the optimal feature combination, appropriate classifiers were selected for image processing. As a widely used method in deep learning, 2D-CNN has been shown to be effective in feature mining and representation. The structure of 2D-CNN is shown in Figure 9, including three convolution layers, three aggregation layers, three batch normalization layers, and two full connection layers. The convolution layer and pool layer are used to build the feature extractor, and the full connection layer is used to build the classifier.

In order to verify the separability of feature combination and the effectiveness of the proposed land-cover classification method based on images from the ascending and descending orbits, three feature combinations of the dual-pol HH-HV mode of Gaofen-3 data were constructed and used as input files. Uncompensated_DOPC_DpRVI and Compensated_DOPC were selected to compare the Compensated_DOPC_DpRVI feature combination of the proposed method.

• Uncompensated_DOPC_DpRVI: The feature combination of DOPC (backscattering coefficient, SPAN, DI, PR, H, A, $\alpha$) and DpRVI based on the uncompensated HH and HV polarization images of Gaofen-3.
• Compensated_DOPC: The feature combination of DOPC based on compensated HH and HV polarization images of Gaofen-3.
• Compensated_DOPC_DpRVI: The feature combination of DOPC and DpRVI based on the compensated HH and HV polarization images of Gaofen-3.

2.4.6. Quantitative Analysis

In order to measure the accuracy of the classification results, quantitative and qualitative evaluations were made for the identification results. To evaluate quantitvely the accuracy of each model, various parameters such as Precision, Recall, F1_Score, OA, and Kappa coefficient are considered.

• Precision:

  $$Precision=\frac{\left(TP\right)}{\left(TP+FP\right)}$$

  (23)
• Recall:

  $$recall=\frac{\left(TP\right)}{\left(TP+FN\right)}$$

  (24)
• F1_Score:

  $$F1\_Score=\frac{\left(2\times precision\times recall\right)}{\left(precision+recall\right)}$$

  (25)
• OA:

  $$OA=\frac{\left(TP+TN\right)}{\left(TP+TN+FP+FN\right)}$$

  (26)
• Kappa coefficient:

  $$Kappa=\frac{\left(OA-p_e\right)}{\left(1-p_e\right)}$$

  (27)

  $$p_e=\frac{\left(TP+TN\right)\times \left(FP+TP\right)\times \left(FP+TN\right)\times \left(FN+TN\right)}{{\left(TP+TN+FP+FN\right)}^2}$$

  (28)

  where TP is a positive sample predicted by the model to be positive, TN is a negative sample predicted by the model to be negative, FP is a negative sample predicted by the model to be positive, and FN is a positive sample predicted by the model to be negative.

Qualitative evaluation was applied to evaluate the classification results intuitively. Compared with optical images, the classification results were interpreted and verified intuitively. The superiority of the proposed method was verified by the comparison of the three feature combinations.

3. Results and Discussion

3.1. Geometric Distortion Region Analysis

Experiments were carried out in the Danjiangkou area, where the complex terrain is mostly mountainous and widespread problems of layover and shadow are encountered. The preprocessed ascending and descending orbits of Gaofen-3 images in the Danjiangkou area are shown in Figure 10. The dimension of the study area is 4611 × 3208 pixels. The HV polarization image of the descending track was selected as the main image, while the HV image of the ascending track was selected as the sub-image, as shown in Figure 10a,b.

Figure 10b,d shows enlargements of areas A and B in red rectangular frames in Figure 10a,c. As can be seen from Figure 10, the backscattering coefficient in the dwell zone is high and the image is bright due to the superposition of the signals, while the backscattering coefficient is low in the areas of shadow which appear dark, due to the loss of signal blocked by mountains. In addition, different images of a specific area were obtained by using SAR sensors in the ascending and descending trajectories. In the same area at the same latitude, the information from another track can be collected from the cross-track images. The geometric distortion areas in the primary image are normal in the secondary image, and the information confusion in the primary image can be corrected by geometric area compensation.

3.2. Geometric Distortion Area Detection

The geometric model of satellite imaging was established by combining the three-dimensional coordinates of DEM and orbital parameters of Gaofen-3. The geometric distortion region of the main image is shown in Figure 11. Figure 11b shows the instantaneous incident angle. Subsequently, the geometric distortion area was determined by morphological processing according to the local incident angle range (Figure 11c). In Figure 11c, white areas represent is the regions of geometric distortion, and black indicates the normal region. It can be observes that the geometric distortion areas are mainly distributed near mountains, seriously affecting the accuracy of land-cover classification.

3.3. Geometric Distortion Area Compensated

The geometric distortion areas in the primary image are processed by the secondary image. Figure 12a,b shows the original polarized images and Figure 12c,d shows the compensation results. An enlarged area was selected to analyze the compensation effect, as shown in Figure 12b,d. Compared with the original polarized image, it can be seen intuitively that most of the geometric distortion areas in Figure 12b have been compensated in Figure 12d.

3.4. J-M Distance Analysis of Features Combinaion

Buildings, water, farmland, and woodland show different separability in different combinations of features. To evaluate the separability of DpRVI and DOPC after compensation, the J-M distances of the samples in three cases were calculated according to different combinations of features, as shown in Figure 13.

According to the J-M distance analysis of the feature combinations, it can be seen that compared with Uncompensated_DOPC_DpRVI and Compensated_DOPC, the Compensated_DOPC_DPRVI proposed in this paper has the highest feature separability. The J-M distance of Compensated_DOPC_DpRVI is almost always above 1.9, which indicates that the sample height can be divided. Geometric distortion compensation and the introduction of DpRVI will improve the separability of different types of ground objects.

Generally speaking, the separability of buildings, waters and other features in Compensated_DOPC is significantly higher than that in UnCompensated_DOPC_DpRVI and Compensated_DOPC_DpRVI. The separability of farmland and forest land in Compensated_DOPC_DpRVI is higher than that in Compensated_DOPC. In other words, the treatment of geometric distortion can improve the recognition accuracy for buildings and water areas, while DpRVI is more effective for woodland and farmland.

3.5. Quantitative Analysis and Discussion

3.5.1. Quantitative Evaluation

The combinations of features were input into the model, respectively, then the predicted results were compared with the actual samples to obtain the confusion matrix and the above precision parameters. The accuracy results are shown in Figure 14 and

Table 3. Quantitative evaluation of the model generated by the combinations of features.

		Building	Farmland	Woodland	Water	OA	Kappa
Uncompensated_DOPC_DpRVI	Precision	0.83	0.86	0.84	0.98	0.88	0.86
	Recall	0.73	0.94	0.86	0.99
	F1_Score	0.78	0.90	0.85	0.99
Compensated_DOPC	Precision	0.90	0.86	0.76	0.99	0.89	0.87
	Recall	0.88	0.92	0.80	0.98
	F1_Score	0.89	0.88	0.78	0.99
Compensated_DOPC_DpRVI	Precision	0.93	0.89	0.90	0.99	0.93	0.92
	Recall	0.91	0.95	0.93	0.98
	F1_Score	0.92	0.92	0.91	0.98

.

In terms of overall accuracy, the Compensated_DOPC_DpRVI feature combination has the highest accuracy, with an OA of 0.93 and Kappa score of 0.92. Compensated_DOPC also achieved good classification results, with an OA score of 0.89 and Kappa score of 0.87, but its accuracy was slightly lower than that of Uncompensated_dopc_DpRVI because DpRVI was not added. Uncompensated_DOPC_DpRVI introduces DpRVI but does not deal with the geometric distortion area, which led to the lowest accuracy: the OA score was 0.88 and Kappa score was 0.86. Experiments show that OA and Kappa can be improved by about 5% after geometric distortion processing and introduction of DpRVI.

In terms of the accuracy of different land-cover types, the accuracy of the Compensated_DOPC_DpRVI method for all land types was higher than that of Uncompensated_DOPC_DpRVI or Compensated_DOPC. In Uncompensated_DOPC_DpRVI, the value for farmland and forest land was slightly higher than Compensated_DOPC, by about 7%. This proves that DpRVI has a good effect on the identification of woodland and farmland.

3.5.2. Qualitative Evaluation

Figure 15 shows the qualitative evaluation results for the three sets of feature combinations and the optical reference images. Generally speaking, the classification results for flat areas are similar, and all the feature combinations can describe the ground objects in the study area. However, there are many geometric distortions in complex terrain areas, and the classification results for the feature combinations are quite different. Two areas in the mountainous region, A and B, were selected for further analysis and magnified for observation. Enlarged classification results are shown in Figure 16.

As shown in Figure 16b,f, due to the influence of layover and shadow, a large number of ground objects were wrongly divided into buildings or water in the classification results based on Uncompensated_DOPC_DpRVI. After the geometric distortion area was compensated, the misclassification of water and buildings improved to some extent; the classification results based on Compensated_DOPC are shown in Figure 16c,g. However, compared with Compensated_DOPC_DpRVI, the Compensated_DOPC did not introduce DpRVI, and the feature separability for different vegetation is weak, resulting in lower classification accuracy for forest land and farmland.

Based on uncompensated _DOPC_DpRVI classification results, due to the influence of layover and shadow, a large number of ground objects are wrongly divided into buildings or water (Figure 16b,f). Compared with Compensated _DOPC_DpRVI, Compensated _DOPC did not introduce DpRVI, although the geometric distortion area was compensated, which led to insufficient classification accuracy for woodland and cultivated land (Figure 16c,g).

In this study, the images of the relative track were used to compensate for the geometric distortion. The compensation of geometric distortion area can replace 86.3% error information, see

Table 4. Compensation situation.

Image	Number of Pixel	Area (Km2)
Geometric distortion region	3,460,644	346.1
Compensation area	2,986,287	298.6
Study area	14,792,088	1479.2

for details.

$$\mathrm{Compensation}\mathrm{Ratio}=\frac{\mathrm{Compensation}\_\mathrm{Num}}{\mathrm{Geometric}\_\mathrm{distorion}\_\mathrm{Num}}$$

(29)

where Compensation_Num is the number of pixels in the compensation area, and Geometric_distorion_Num is the number of pixels in the geometric distortion area.

However, due to the imaging angle of the SAR sensors, even sensors with opposite orbits cannot guarantee that all the collected information is correct. As shown in Figure 17, the image is compensated by lifting the track data, but there are still a few geometric distortion areas that cannot be compensated.

4. Discussion

In this study, the classification of ground objects in complex terrain area was studied by using ascending and descending orbit HH-HV dual-pol mode SAR data from the Gaofen-3 satellite. Through geometric distortion processing and effective feature combination, the classification confusion relating to complex terrain areas can be effectively improved.

There are vast mountains and hills on the earth’s land surface, and the SAR side-looking imaging mechanism leads to serious geometric distortion in areas of complex terrain. The experimental results show that it is difficult to identify land objects using the scattered signal of the distorted area without processing. In this study, 86.3% of the geometric distortion was compensated by the data of the lifting orbit, greatly improving the incorrect scattering information in the geometric distortion area.

Because of the imaging angle of the SAR sensors, even sensors with opposite orbits cannot guarantee that all the collected information is correct. The current data for the research area are Gaofen-3 dual-track data for compensation and fusion, but there are still a few geometric distortion areas that cannot be compensated. Orbit data from different angles should be added to ensure the integrity of collected information.

Furthermore, effective feature combination is an important guarantee of classification accuracy. In addition to enhancing the separability between land-cover types, DpRVI has been introduced and assessed. DpRVI is mainly used to reflect the growth state of vegetation and the inversion of soil moisture, while its role in land-cover classification has rarely been investigated, especially in relation to HH-HV dual-pol mode SAR data. It was found that the introduction of DpRVI in HH-HV dual-pol mode is effective and can significantly increase the J-M distance between woodland and farmland. The application of DpRVI deserves more attention in follow-up research into land-cover classification.

Finally, a 2D-CNN classification model was constructed. The results show that the classification accuracy was greatly improved by geometric distortion compensation and effective feature combination with DpRVI.

Existing research mainly focuses on detecting areas of terrain distortion. In view of the present problem of classification errors caused by geometric distortion, this study explores the identification of the correct types of objects in geometric distortion areas through compensation by reverse orbit data. This research verifies the effectiveness of DpRVI in land-cover classification, especially for farmland and forest areas.

5. Conclusions

In this study, a land-cover classification method based on Gaofen-3 lifting rail data is proposed, which determines the specific types of ground objects on the basis of detecting geometric distortion areas. Firstly, by establishing the geometric model of SAR imaging, the layover and shadow areas are detected. Then, the geometric distortion region is compensated by using different trajectory data. Finally, effective feature combinations with DpRVI are extracted from the compensated image, and a 2D-CNN model is constructed to realize classification. The effectiveness of this method is verified by qualitative and quantitative evaluation.

In the quantitative evaluation, compared with the results of other methods, F1 score, OA, and Kappa were all improved about 5%. All the evaluation indicators are above 90%. In qualitative evaluation, two typical mountain areas were selected for visual observation, and were very similar to the optical images. It is proved that this method is applicable in complex terrain.

Briefly, the land-cover classification method based on ascending and descending orbit-data compensation in complex terrain areas can improve the separability of polarization characteristics and the classification accuracy. In addition, the dual polarization data required in this method are simple and easy to obtain. Compared with full polarization, land-cover classification based on dual-pol data lacks the expression of a vegetation index. Experiments show that the introduction of DpRVI can improve the classification accuracy based on double pol data. It broadens the application scope of SAR remote sensing technology, and has good application prospects in the field of land-cover classification.