GNSS Ground deformation observation network optimization assisted using prior InSAR-derived ground surface deformation and multiscale iteration estimation

ABSTRACT

Ground surface deformation must be monitored to understand and prevent geological hazards. Global navigation satellite system (GNSS) technology has been widely used in deformation monitoring; however, current GNSS station layout strategies are often subjective and lack theoretical guidance. Therefore, in this study, a GNSS deformation observation network optimization method was developed based on interferometric synthetic aperture radar (InSAR) prior to deformation assistance. This method considers the deformation distribution, cost, accuracy, topography, and actual monitoring requirements and uses kriging interpolation and multiscale iterative optimization to obtain the reference number, spatial distribution, and deformation capture accuracy of GNSS stations in areas experiencing deformation. Simulation experiments showed that the proposed method strengthened the spatial deformation monitoring ability of GNSS by 46%–68% compared with the traditional station layout methods, which ignore the distribution of the deformation. Then, based on the small baseline subset (SBAS) InSAR technology and applying the new method, we evaluated the number and location of the most suitable GNSS stations for a mining area and found that 89.6% of the error was less than 1.5 cm, which further illustrates the practicability and reliability of the proposed method. This method provides an effective supplement to station layouts for monitoring GNSS deformation.

KEYWORDS:

• Ground deformation monitoring
• GNSS station layout
• InSAR
• multiscale optimization

1. Introduction

The Global Navigation Satellite System (GNSS), as one of the most mature precision positioning technologies, provides all-weather, automated, long-term, and continuous observations, so it has been widely used in engineering surveys and for natural disaster monitoring (Petropoulos and Srivastava 2021). However, the layout of stations for deformation monitoring using GNSS technology is mainly determined with empirical guidance; moreover, the spatial density of GNSS stations is limited by the cost of the ground equipment (Zulkifli, Din, and Som 2018), so the results of GNSS deformation monitoring do not accurately reflect the spatial distribution of deformation in an area. Interferometric synthetic aperture radar (InSAR) (Massonnet and Feigl 1998; Zhu, Li, and Hu 2017) technology, with its high spatial resolution, high precision, and large-scale rapid monitoring, after decades of development has made considerable breakthroughs in various fields such as the monitoring of volcanoes, earthquakes, landslides, mining areas, glaciers, and urban infrastructure, and for atmospheric correction (Soldato et al. 2021; Xue et al. 2020). InSAR is one of the effective methods used for monitoring geological disasters. However, the temporal resolution of InSAR, when used for monitoring deformation, is limited by the revisit period of the sensors, which is 6–24 days for most satellites at present (Macchiarulo et al. 2022). When monitoring deformation in mining areas (Yang et al. 2020), for landslides (Rosi et al. 2018; Zheng et al. 2023), and groundwater (Motagh et al. 2017), which have nonlinear characteristics and large deformation gradients, fully reflecting the temporal characteristics of their deformation is not yet impossible. Therefore, InSAR deformation results cannot be easily used for timely and accurate judgments in some disaster emergency monitoring and engineering applications.

Given these problems, many researchers have focused on the integration of InSAR and GNSS technologies by, for example, using GNSS to correct the atmospheric delay of the InSAR troposphere (Houlie, Funning, and Bürgmann 2016; Li et al. 2006; Reuveni et al. 2015), combining GNSS and InSAR to retrieve three-dimensional surface deformation (Gudmundsson, Sigmundsson, and Carstensen 2002; Samsonov et al. 2007), and fusing InSAR and GNSS deformation results to obtain high spatial and temporal resolutions (Ge, Han, and Rizos 2000; Liu et al. 2019; Yan et al. 2022). These technologies and their applications have made substantial progress. In these studies, the results of spatial monitoring obtained with GNSS were used, but few studies have considered the required layout of GNSS observation stations to optimize InSAR and GNSS fusion or how to balance cost and monitoring accuracy. The accuracy of deformation monitoring increases with the number of GNSS points. However, increasing the number of GNSS points raises the costs. Therefore, the amount of GNSS equipment used and the accuracy of the deformation monitoring must be balanced. Relevant research in this area remains relatively scarce. Zou et al. (2019) designed a layout strategy for GNSS stations for strike–slip faults based on physical models and parameters such as site spacing and minimum effective distance, and they discussed the influence of the GNSS station distribution on the slip rate and locking depth. Sarsito et al. (2019) designed an ideal GNSS network by measuring the three-dimensional deformation of a volcanic surface using physical equations, such as the Mogi and Okada models. These model-driven methods require prior knowledge of the physical characteristics or equations of the monitored object that cannot be obtained through general monitoring, hindering their generalization and applicability.

The layout of traditional GNSS observation station for deformation monitoring has relied mostly on experience, lacking a theoretical basis. This GNSS observation station layout strategy is not a sufficient scheme: it cannot reflect the spatial distribution characteristics of the deformation. Although adding additional equipment can reduce the estimation error of the parameters, it increases the overall system cost. Some individual methods rely on physical parameters, which are difficult to obtain and thus these methods lack generalizability. This is primarily because of the past lack of technical methods for understanding and evaluating large-scale deformations, resulting in the selection of GNSS station locations by mainly relying on experience and data surveys. In recent years, the launching of various synthetic aperture radar (SAR) satellites (Xue et al. 2020) and the provision of free data (Torres et al. 2012) have increased the efficiency and precision of InSAR technology while decreasing its cost when determining surface deformation, so the GNSS deformation monitoring network can now be deployed through a data-driven method with the assistance of prior InSAR deformation information. However, to the best of our knowledge, no such study has yet been conducted. Therefore, a system is needed for determining the number and spatial distribution of GNSS deformation monitoring stations with the support of an a priori InSAR high-spatial-resolution deformation field.

In this study, we developed a GNSS observation network Layout optimization method based on InSAR Deformation Assistance (IDAL). Based on prior InSAR deformation results, and considering the actual monitoring requirements, cost, and spatial deformation recovery accuracy, this method accurately obtains the reference number, spatial distribution, and the corresponding spatial deformation of the GNSS station layout in areas experiencing deformation through multiscale iterative optimization. The remainder of this paper is organized as follows: Section 2 introduces the InSAR deformation retrieval method and the GNSS deformation observation network optimization method, which uses prior InSAR deformation information. Section 3 presents a verification of the reliability and practicability of the proposed method, which was achieved through simulations and real data experiments. Section 4 describes the evaluation of the influence of certain factors on the proposed method. Finally, Section 5 concludes the paper.

2. Methods

2.1. Prior deformation retrieval based on SBAS–InSAR technology

Through decades of development, InSAR technology has gradually matured. Various advanced InSAR techniques (Xue et al. 2020) can be used to determine the surface deformation rate and time-series deformation. In this study, the SBAS–InSAR (Berardino et al. 2002; Lanari et al. 2004) method was used to obtain deformation, as only interferograms with short spatiotemporal baselines are required for calculation, which can suppress the influence of spatiotemporal decorrelation to a certain extent.

Assuming $K+1$ SAR images, $E$ interferograms are generated according to the spatiotemporal baseline threshold, with only one small baseline set. (1) $\left[\begin{array}{c}\mathrm{\Delta }{\phi }_1\\ \mathrm{\Delta }{\phi }_2\\ ⋮\\ \mathrm{\Delta }{\phi }_j\\ ⋮\\ \mathrm{\Delta }{\phi }_E\end{array}\right]={\displaystyle \frac{4\pi }{\lambda }\cdot \left[\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\\ 1& 0& -1& \cdots & 0& 0\\ 0& 1& -1& \cdots & 0& 0\\ & & & ⋮& & \\ 0& 0& 0& \cdots & 1& -1\end{array}\right]\times \left[\begin{array}{c}{d}_1\\ d_2\\ ⋮\\ d_i\\ ⋮\\ d_K\end{array}\right]}$(1) The solution of Equation (1(1) $\left[\begin{array}{c}\mathrm{\Delta }{\phi }_1\\ \mathrm{\Delta }{\phi }_2\\ ⋮\\ \mathrm{\Delta }{\phi }_j\\ ⋮\\ \mathrm{\Delta }{\phi }_E\end{array}\right]={\displaystyle \frac{4\pi }{\lambda }\cdot \left[\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\\ 1& 0& -1& \cdots & 0& 0\\ 0& 1& -1& \cdots & 0& 0\\ & & & ⋮& & \\ 0& 0& 0& \cdots & 1& -1\end{array}\right]\times \left[\begin{array}{c}{d}_1\\ d_2\\ ⋮\\ d_i\\ ⋮\\ d_K\end{array}\right]}$(1) ) can be obtained using the least-squares method, which is expressed as: (2) $\hat{x}=(A^T\mathrm{PA}{)}^{-1}{A}^T\mathrm{Pb}$(2) where $\hat{x}=[d_1,d_2,\cdots ,d_i,\cdots ,d_K{]}^T$ is the final deformation result after the adjustment. $A={\displaystyle \frac{4\pi }{\lambda }\cdot \left[\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\\ 1& 0& -1& \cdots & 0& 0\\ 0& 1& -1& \cdots & 0& 0\\ & & & ⋮& & \\ 0& 0& 0& \cdots & 1& -1\end{array}\right]}$ is the coefficient matrix, $b=[\mathrm{\Delta }{\phi }_1,\mathrm{\Delta }{\phi }_2,\cdots ,$ $\mathrm{\Delta }{\phi }_j,\cdots ,\mathrm{\Delta }{\phi }_E{]}^T$ is the observation matrix, and$P$ is the weight matrix, which can be determined from the coherence of the interferogram.

After obtaining the InSAR deformation results, the GNSS station layout was determined based on the mean InSAR deformation rate because the mean deformation rate comprehensively reflects the deformation over a time period, which shows not only the spatial distribution of deformation but also the temporal accumulation of deformation.

InSAR technology can only be used to obtain one-dimensional line-of-sight (LOS) deformation. According to the SAR data coverage, multitrack InSAR deformation results can be obtained to perform a three-dimensional surface deformation reconstruction.

2.2. Establishment of optimization objective function using prior InSAR deformation

The purpose in this study was to enable GNSS monitoring to more accurately reflect the spatial characteristics of deformation. Therefore, InSAR deformation monitoring results were used as a reference; the deformation at a station monitored with GNSS was interpolated and compared with the InSAR deformation results, and the minimum root mean square error (RMSE) was taken as the optimization goal.

Assuming that set $D_n$ indicates that $n$ GNSS station candidate points (CPs) were obtained in study area D, the selection of $N$ points from $D_n$ for the layout of GNSS stations is required. The following symbols are defined:

• $P_N$ represents the set of the distribution of these $N$ GNSS stations;

• $(x_i,y_i,d_i)$ denote the position and deformation of the $i$th point in the current point set, where $i=1,2,\dots N$, which can be simplified as $P_N^i$;

• $S_M$ represents a set of $P_N$ with the same number of GNSS stations but different distributions, where $M$ denotes $M$ sets $P_N$ in the set, then $S_M=\{{P_k}_N|k=1,2,\dots M\}$, and $S_M^k={P_k}_N$, which all represent the $k$th station distribution set.

Now all $P_N$ in $S_M$ are interpolated to recover the spatial deformation field. As Kriging interpolation (Shahbeik et al. 2014) can be used to perform an unbiased optimal estimation of the regionalized variables, all $P_N$ sets in $S_M$ are subjected to Kriging interpolation (the process of recovering the spatial deformation field through GNSS stations via Kriging interpolation is referred to as the deformation recovered by GNSS stations in this paper). Prior to this, the degree of variation in the spatial deformation field must be estimated. (3) $\hat{\gamma }(h)={\displaystyle \frac{1}{2N(h)}\sum _{i=1}^{N(h)}\{z(l_i)-z(l_i+h){\}}^2}$(3) where$h$ is the distance between the InSAR points, $z(l_i)$ represents the deformation value at $l_i$; $N(h)$ is the number of pairs of all observation points, with $h$ being the interval; and $\hat{\gamma }(h)$ is the semivariance of the interval $h$, which increases with the increase in $h$ within a certain range.

An appropriate variogram model (Knospe and Jonsson 2010) was used to interpolate the deformation field. The recovered deformation field was then compared with the InSAR deformation field, and the RMSE of each GNSS distribution set was calculated as follows: (4) $\mathrm{RMS}{E}_k=\sqrt{{\displaystyle \frac{1}{N(D)}\sum _{i=1}^{N(D)}(z_i-{\hat{z}}_i)}}$(4) where $z_i$ represents the deformation monitored using InSAR, ${\hat{z}}_i$ represents the recovered deformation from the GNSS stations, $N(D)$ denotes the number of InSAR points in monitoring area D, and $\mathrm{RMS}{E}_k$ represents the RMSE of the $k$th point set in $S_M$.

The objective function is a function of the position of each station. Therefore, using the RMSE calculated using Equation (4(4) $\mathrm{RMS}{E}_k=\sqrt{{\displaystyle \frac{1}{N(D)}\sum _{i=1}^{N(D)}(z_i-{\hat{z}}_i)}}$(4) ) as an indicator, the problem of solving the position of the GNSS stations can be converted into the problem of minimizing the RMSE: (5) $\begin{array}{r}\mathit{X}\end{array}=\{(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\}=\arg \min J[(x_1,y_1,v_1),(x_2,y_2,v_2),\dots ,(x_N,y_N,v_N)]=\arg \underset{R}{min}\{W\cdot R(x_1,y_1,v_1),(x_2,y_2,v_2),\dots ,(x_N,y_N,v_N)\}$(5) where $X$ represents the optimal solution of the station position, $J$ is the objective function with the coordinates of each GNSS point as the parameters to be solved, $W$ is the weight or scale factor, and $R$ represents the RMSE between the deformation recovered from each potential GNSS station distribution and the InSAR deformation.

2.3. Determination of GNSS CP set

Before solving the model in Equation (5(5) $\begin{array}{r}\mathit{X}\end{array}=\{(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\}=\arg \min J[(x_1,y_1,v_1),(x_2,y_2,v_2),\dots ,(x_N,y_N,v_N)]=\arg \underset{R}{min}\{W\cdot R(x_1,y_1,v_1),(x_2,y_2,v_2),\dots ,(x_N,y_N,v_N)\}$(5) ), CPs set of the GNSS stations must be determined. As the computational burden of using all InSAR deformation points as CPs would be huge, quadtree downsampling is required (Decriem et al. 2010) for the InSAR deformation points, which enables the selection of GNSS station CPs that reflect the spatial characteristics of deformation and reasonably reduces the number and density of InSAR points.

Specifically, the InSAR deformation field is continuously subdivided into squares according to the quadtree recursive algorithm, and the variance in the InSAR deformation in each square is calculated as follows: (6) ${\sigma }^2={\displaystyle \frac{{\sum }_{k=1}^q{(z_k-\overline{z})}^2}{q}}$(6) where ${\sigma }^2$ is the variance of the InSAR points in the square, $q$ is the number of InSAR points in the square, $z_k$ denotes the deformation of the $k$th point in the square, and $\overline{z}$ is the mean value of $q$ InSAR deformation points in the square.

Set the segmentation threshold δ. If the variance is greater than δ, continue to divide until the variance of the deformation in each square is less than or equal to δ. The points after quadtree downsampling are taken as the CPs.

Not every point is suitable for a GNSS station. Simultaneously, contemplating future equipment maintenance, the available geological structure, topography, and other external data must be integrated. Subsequently, areas with excessive terrain undulation and overly lush vegetation were masked.

2.4. Multiscale iterative solution

Even with the use of quadtree downsampling to reduce the number and density of candidate GNSS points, attempting to directly generate different distributions of GNSS stations through permutation and combination and subsequently applying Kriging interpolation to recover deformations for solving the models in Equation (5(5) $\begin{array}{r}\mathit{X}\end{array}=\{(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\}=\arg \min J[(x_1,y_1,v_1),(x_2,y_2,v_2),\dots ,(x_N,y_N,v_N)]=\arg \underset{R}{min}\{W\cdot R(x_1,y_1,v_1),(x_2,y_2,v_2),\dots ,(x_N,y_N,v_N)\}$(5) ) would become challenging. As the number of solutions exponentially increases, directly traversing all possible combinations become impractical. Some heuristic searches, such as the greedy algorithm (DeVore and Temlyakov 1996) and simulated annealing (Bertsimas and Tsitsiklis 1993) fall into local solutions, and the solution is unstable. Therefore, considering the efficiency and accuracy of the solution, a multiscale iterative optimization method was developed to solve the proposed model.

The initial station distribution must be determined before solving the model. This is recorded as $Q_N$, where $Q_N^i$ represents the $i$th point at the initial station. There are two types of stations in $Q_N$, one of which must be monitored according to actual requirements, and these stations will not change in the future; the other includes stations that require iterative optimization, and the initial positions of these stations can be roughly selected according to the InSAR deformation.

Next, a multiscale iterative optimization method is applied to solve the objective function model according to the InSAR deformation, $Q_N$, CPs, and variogram. First, the initial station points are iteratively solved on a large scale (i.e. with a small number of CPs). One method of obtaining large-scale CPs is to increase the variance threshold of the quadtree downsampling, as described in Section 2.3, to reduce the number of CPs. Another method involves downsampling the InSAR deformation field by reducing the resolution to reduce the number of CPs, as shown in Figure 1.

Figure 1. CPs downscaling using (a) different quadtree downsampling thresholds and (b) different downsampling resolutions.

Assume that the large-scale CP set $D_n^{{}^{\mathrm{\prime }}}$ is obtained through the first method and that its density and number are much smaller than those of $D_n$. Now, the optimal solution $X_N^{{}^{\mathrm{\prime }}}$ of $Q_N$ for large-scale $D_n^{\mathrm{\prime }}$ can be iteratively solved.

First, the RMSE between the spatial deformation recovered from the $Q_N$ and the InSAR deformations is calculated and denoted as $r_0$. Second, for each point in $Q_N$, excluding the fixed point, a search is conducted in $D_n^{\mathrm{\prime }}$ for points with a distance less than $T_1$, forming a new point set $C_{n1}$, where $n1$ indicates $n1$ points in $C_{n1}$. The current station in $Q_N$ is replaced by each point in $C_{n1}$ in turn, and the RMSE $r_n$ between the spatial deformation recovered from the formed $n1$ station distributions and the InSAR deformation is calculated. Finally, we obtain (7) $r=r_n-r_0.$(7) The minimum value $r_{\min }$ in $r$ and its corresponding station distribution $C_{n1}^j$, $j$ are recorded, indicating that the minimum $r_{\min }$ corresponds to the $j$th station distribution in $C_{n1}$. (8) $Q_N^{\mathrm{\prime }}=\{\begin{array}{c}{C}_{n1}^j{r}_{\min }<0\\ Q_N{r}_{\min }\ge 0\end{array}$(8) As shown in Equation (8(8) $Q_N^{\mathrm{\prime }}=\{\begin{array}{c}{C}_{n1}^j{r}_{\min }<0\\ Q_N{r}_{\min }\ge 0\end{array}$(8) ), when $r_{\min }$<0, the station distribution $C_{n1}^j$ is more appropriate than the current station distribution. Thus, $C_{n1}^j$ is used to replace the current station distribution; otherwise, the current station remains unchanged. This process is repeated for all nonfixed stations until the distance between the corresponding positions of each point in the two iteration stations is less than a given threshold. The final solution can be expressed as (9) $\{\begin{array}{c}{\mathit{X}}_{\mathit{N}}^{\mathrm{\prime }}=\arg minr\{(x_1^{\mathrm{\prime }},y_1^{\mathrm{\prime }}),(x_2^{\mathrm{\prime }},y_2^{\mathrm{\prime }}),\dots ,(x_N^{\mathrm{\prime }},y_N^{{}^{\mathrm{\prime }}})\}\\ (x_1^{{}^{\mathrm{\prime }}},y_1^{{}^{\mathrm{\prime }}}),(x_2^{{}^{\mathrm{\prime }}},y_2^{{}^{\mathrm{\prime }}}),\dots ,(x_N^{{}^{\mathrm{\prime }}},y_N^{{}^{\mathrm{\prime }}})\in {\mathit{D}}_{\mathit{n}}^{{}^{\mathrm{\prime }}}\end{array}.$(9) The final station distribution ${\mathit{X}}_{\mathit{N}}^{{}^{\mathrm{\prime }}}$ is the optimal solution for the initial station $Q_N$ for large-scale ${\mathit{D}}_{\mathit{n}}^{\mathbf{\prime }}$.

${\mathit{X}}_{\mathit{N}}^{{}^{\mathrm{\prime }}}$ is still a rough solution relative to the final station, but its approximate position is correct. Therefore, only a smaller search threshold needs to be used to iterate further on a small-scale point with denser CPs. Similar to the process of solving for station ${\mathit{X}}_{\mathit{N}}^{{}^{\mathrm{\prime }}}$, the initial station $Q_N$ needs to be replaced with ${\mathit{X}}_{\mathit{N}}^{{}^{\mathrm{\prime }}}$, where the large-scale CPs set ${\mathit{D}}_{\mathit{n}}^{\mathbf{\prime }}$ is replaced with the small-scale candidate point set ${\mathit{D}}_{\mathit{n}}$, and the larger distance threshold $T_1$ is replaced with the smaller distance threshold $T_2$. Then, the optimization of large-scale stations is repeated to obtain the final GNSS station optimization result ${\mathit{X}}_{\mathit{N}}$. (10) $\{\begin{array}{c}{\mathit{X}}_{\mathit{N}}=\arg minr\{(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\}\\ (x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\in {\mathit{D}}_{\mathit{n}}\end{array}$(10) Then, the IDAL method is used to calculate the RMSE of the optimal station distribution corresponding to different numbers of stations and to determine the corresponding relationship between the number of stations and the RMSE. The inflection point of the curve represents the optimal number of GNSS stations in the area experiencing deformation. Figure 2 presents a flowchart of the IDAL method.

Figure 2. Flow chart of GNSS station layout optimization based on InSAR deformation.

3. Results

3.1. Simulation experiments

A simulation was conducted to verify the effectiveness of the IDAL method. Figure 3 shows a simulated deformation field (200 $\times$ 300 pixels; resolution 10 $\times$ 10 m). To make the deformation more realistic, local subsidence and Gaussian white noise were added.

Figure 3. Simulated deformation.

3.1.1. Optimal distribution of GNSS stations based on IDAL method

First, the stable area was masked according to the InSAR deformation field to reduce unnecessary calculations. Second, employing the IDAL method, the station layout results of the deformation area were obtained. Finally, taking the deployment of nine stations in this area as an example, traditional methods, such as the empirical layout (EL) and uniform layout (UL) methods, were compared with the IDAL method (as shown in Figure 4). The EL method is an empirical GNSS station layout method that initially considers the deformation distribution. The UL method evenly lays out stations in areas with deformation, and their number and spacing are determined using human experience.

Figure 4. Accuracy comparison of different station layout methods (black triangles: fixed stations; red triangles: final optimized stations). (a–c) Deformation results recovered using EL, UL, and IDAL methods; (d–f) the between (a–c) and the simulated deformation; (g–i) error histograms.

Figure 4(a–i) show the deformation fields and their error distribution results recovered using the three methods based on the corresponding nine stations. The deformation obtained using these three methods was similar to the simulated deformation, indicating that they could recover the general characteristics of the deformation. However, their error distributions and histograms significantly differed. The IDAL method fully considers the distribution of deformation; therefore, it accurately captured the spatial characteristics of deformation in the entire region using only nine stations, and the error was small and almost obeyed a normal distribution. The results of the EL method were primarily arranged near the two central axes, which was reasonable. However, compared with the IDAL method, the accuracy on the axis was not significantly improved, and the deformation error of some areas exceeded 4 cm. Similarly, the local error produced with the UL method was also large. In addition, the errors of the EL and UL methods were more dispersed.

Table 1 lists the mean absolute error (MAE), standard deviation (STD), and RMSE of the three methods. Overall, the IDAL method exhibits the smallest errors. Compared with the EL method, its MAE, STD and RMSE increased by 54.96%, 46.32 and 55.28% respectively; compared with the UL method, its performance in MAE, STD and RMSE increased by 66.67%, 61.17 and 68.42% respectively. These results demonstrated the rationality and reliability of the proposed method.

Table 1. Statistical error analysis of EL, UL and IDAL methods from Figure 4(g–i).

Method	MAE (cm)	STD (cm)	RMSE (cm)
EL	1.31	1.36	1.61
UL	1.77	1.88	2.28
IDAL	0.59	0.73	0.72

3.1.2. Optimal number of GNSS stations

To evaluate the most appropriate number of stations for an area experiencing deformation, the variation in the RMSE for different numbers of stations was obtained using the IDAL method. As shown in Figure 5, the shape of the curve approximates that of an L curve (Hansen and O’Leary 1993). With an increase in the number of stations, the increase in deformation accuracy was strong at the beginning but eventually tended to flatten.

Figure 5. RMSE corresponding to different number of stations (blue dotted line: RMSE with 9 stations; green ellipse: reference range of the number of stations in deformation area).

Because a large number of stations cannot be laid, regardless of cost, a compromise is required between the number of stations and accuracy. The green ellipse in Figure 5 presents the reference range of the number of stations for the simulated deformation given by the new method. Within this range, a quality standard for deformation monitoring can be designed, and the number of stations can be reasonably determined according to the actual monitoring requirements and budget. Generally, when the deformation accuracy was higher than 1 cm, most of the monitoring requirements were achieved; therefore, only nine stations were deployed in this simulation experiment (Figure 4 g).

The above simulation experiments and analyses demonstrate that the traditional EL and UL methods, which ignore the deformation distribution, are unreasonable. The proposed IDAL method fully considers the distribution of the deformation field, which maximizes the ability of GNSS stations to capture spatial deformation.

3.2. Real test case: mining deformation monitoring in Datong

3.2.1. SBAS–InSAR ground surface deformation retrieval and analysis

Datong is the northernmost city in Shanxi province, China. To promote urban energy conservation, emission reductions, and local economic development, Datong has built many photovoltaic facilities, including photovoltaic panels, booster stations, and transmission towers, for power generation, and has achieved good results. However, owing to the ongoing coal mining activities in Datong, ground subsidence has occurred in many places (Jia, Ma, and Zhao 2008; Zhao et al. 2014), which seriously threatens the safety of these photovoltaic facilities. Numerous researchers (Palamà et al. 2022; Pawluszek-Filipiak and Borkowski 2020; Pawłuszek-Filipiak et al. 2023) have employed InSAR technology to monitor the ground deformation in mining areas. These studies have considerably contributed to reducing temporal and spatial decoherence, eliminating atmospheric interference, restoring complete deformation information, and improving the overall accuracy of deformation monitoring in mining areas. However, the revisit period of InSAR data is limited. Therefore, in order to grasp the safety and deformation of photovoltaic equipment in the region in a timelier manner, GNSS stations must be deployed in key locations to enable collaborative monitoring on the basis of InSAR deformation monitoring.

Figure 6(a) shows an optical image of the mining area in Datong, where several photovoltaic devices were located. DT05 was located at a booster station; therefore, a GNSS station was deployed to monitor the local deformation. DT03 is another GNSS station located far from DT05, in an area with factories, which monitored the local deformation.

Figure 6. (a) Optical image of the study area (red box: coverage of ALOS2 image; yellow polygon: factory area with installed solar photovoltaic panels; white stars: transmission towers; black triangles: GNSS stations that were deployed; white dotted line: boundary of mining area), and (b) annual mean deformation rate obtained with InSAR (black box: selected deformation region).

In the study area, 11 ALOS2 ascending SAR images from June 19, 2020 to April 23, 2021 and the Shuttle Radar Topography Mission (SRTM) 30 m DEM data were obtained. After registration, a baseline set was formed, with a temporal baseline threshold of 126 days and a spatial perpendicular baseline threshold of 270.8 meters. Then, a multi-look of 1:2 in the range and azimuth directions was used for differential interference, and the obtained differential interferograms were filtered using the Goldstein method (Goldstein and Werner 1998). In general, the phase quality in the center of a mining area is relatively low, which causes difficulties in phase unwrapping. Therefore, the hierarchical networking and constrained adjustment (HNCA) phase unwrapping method (Mao et al. 2021) was used to recover the complete phase with low coherence at the center of the mining area. The orbital error was removed using polynomial fitting (Ferretti et al. 2007). The least-squares method was used to obtain the final mean deformation rate. For the atmosphere in nonlinear deformation, spatiotemporal filtering (Hooper et al. 2012; Osmanoğlu et al. 2016) was used to separate the atmosphere; finally, the complete time series deformation was obtained. Figure 6(b) shows the annual mean deformation rate in the LOS direction obtained using SBAS–InSAR technology.

Subsequently, to compare the InSAR and GNSS results, the GNSS deformation was projected in the LOS direction using Equation (11(11) $D_{\mathrm{gnss}}^{\mathrm{los}}=D_U\cdot \mathrm{cos\theta }-D_E\cdot \mathrm{sin\theta }\cdot \sin \left(\alpha -{\displaystyle \frac{3\pi }{2}}\right)-D_N\cdot \mathrm{sin\theta }\cdot \cos \left(\alpha -{\displaystyle \frac{3\pi }{2}}\right)$(11) ) (Hanssen 2001); DT06, a station far from the deformation area, was employed as a reference station to correct the InSAR deformation results. (11) $D_{\mathrm{gnss}}^{\mathrm{los}}=D_U\cdot \mathrm{cos\theta }-D_E\cdot \mathrm{sin\theta }\cdot \sin \left(\alpha -{\displaystyle \frac{3\pi }{2}}\right)-D_N\cdot \mathrm{sin\theta }\cdot \cos \left(\alpha -{\displaystyle \frac{3\pi }{2}}\right)$(11) where $D_{\mathrm{gnss}}^{\mathrm{los}}$ represents the deformation obtained with the GNSS transformed into the LOS direction; $D_E$, $D_N$ and $D_U$ represent the east–west, north–south, and vertical deformation of the GNSS points, respectively; $\theta$ represents the incident radar angle of the satellite; and $\alpha$ represents the radar azimuth angle.

In Figure 7, the InSAR time-series deformation is compared with that of the two GNSS stations in the study area, and the corresponding RMSE was calculated. The RMSE of DT03 and DT05 were found to be 0.97 and 0.54 cm, respectively, which demonstrate the reliability of the deformation results obtained with InSAR.

Figure 7. Comparison of the time-series deformation accuracy of the InSAR and GNSS results at DT03 and DT05.

Figure 6(b) shows a large subsidence basin in the study area, with the maximum subsidence rate exceeding –15 cm/year. Simultaneously, some transmission towers were at the edge of the subsiding basin, which may cause disasters such as ground fissures. Figure 8 shows that on June 25, 2021, a wall that cracked and sank by 15 cm was observed near DT05, indicating that the subsidence was serious. The ALOS2 data have a spatial resolution of 6 m, but their temporal resolution is 14 days. As such, the satellite sometimes fails to capture data, resulting in only one SAR scene per month. Although Sentinel-1 data can be obtained for free, their spatial resolution is too low (approximately 20 m) to meet monitoring requirements. Therefore, to more accurately identify the subsidence in the area and ensure the safety of the photovoltaic power generation equipment, more GNSS stations are required for joint monitoring.

Figure 8. A wall sinking 15 cm near station DT05 (Level 1 represents the original position of the wall corner, while Level 2 indicates the corner position after sinking 15 cm).

3.2.2. GNSS observation station optimization using IDAL method

Because the spatial distribution of the deformation was ignored when the GNSS stations were laid out, some of the GNSS stations indicated in Figure 6(a) are not reasonable, so they could not effectively monitor the deformation. To enable the new GNSS stations to more accurately capture the spatial characteristics of the deformation, the IDAL method was used for evaluating the number and locations of the most suitable GNSS stations in this area. Figure 9(a) shows the clipped deformation area (black rectangular area in Figure 6(b)); the stable area without deformation and the area away from the plots (purple polygon region) are masked out. Because DT05 was located near an important facility booster station at the center of the local subsidence, a GNSS station is needed here. DTM1 and DTM2 are two local subsidence centers, where the annual mean deformation rate also exceeds –16 cm/year, so GNSS stations must also be deployed at these two locations.

Figure 9. (a) Deformation region of interest, and (b) RMSE corresponding to the number of different GNSS stations (blue dotted line: deformation recovery accuracy with 12 stations; green ellipse: recommended range for the number of stations in this area).

Considering the above three stations as constraints, the IDAL method was used to determine the change in the RMSE of the deformation recovery accuracy with the number of stations in the area. As shown in Figure 9(b), the result is similar to that of the simulation experiment. Compared with the numbers found in the simulation, more stations are required to achieve the same accuracy in this deformation area because of the three local settlement centers in this area, and the deformation is more complicated. However, this reflects the universality of the IDAL method: it can handle complex scenes. The green ellipse in Figure 9(b) represents the reference range for the number of stations in the region based on the curve. Considering the high cost of GNSS ground equipment and the substantial mining deformation, an RMSE of less than 1 cm can satisfy the deformation monitoring requirements. Therefore, 12 GNSS stations are required in this area, which would produce an RMSE of 0.95 cm at this time.

As the InSAR-based deformation results were previously verified, these results for these stations were used to recover the deformation field. Figure 10(a) shows the deformation field recovered using these 12 GNSS stations (three black triangles and nine red triangles). The black triangles denote the constrained stations, while the red triangles represent station positions optimized by the IDAL method. The deformation field retrieved from the GNSS stations aligns with the deformation characteristics observed in the InSAR deformation field depicted in Figure 9(a). Figures 10(b) and 10(c) show that the error obeys an approximately normal distribution, which indicates that the spatial distribution of these stations is reasonable. The statistical results of the absolute error percentage in Table 2 show that 89.6% of the point errors were less than 1.5 cm, which indicates that the deformation recovered using the 12 stations is reliable and can be used to determine the deformation in the entire study area. Of the point errors, 4.8% were greater than 2 cm because the deformation in this area was more complicated. However, given the few photovoltaic devices in these areas, this error is acceptable.

Figure 10. (a) The station-recovered deformation (red triangles: stations optimized using IDAL method), (b) the difference between the deformation recovered using the GNSS stations and the InSAR deformation, and (c) distribution histogram of the error in b.

Table 2. Statistics of the absolute error from Figure 10(b).

Absolute error (cm)	<1.0	1∼1.5	1.5∼2	>2
Percentage (%)	76.0	13.6	5.6	4.8

Figure 10(a) also shows that although DT03 was not selected as a necessary station, it almost coincides with the optimized station DTB2, which indicates that the position of station DT03 is reasonable. Five stations, DTB1, DTB2, DTB6, DTB7, and DTB8, are located inside or near the equipment; therefore, they can be used for the high-precision monitoring not only of the overall deformation but also of the deformation of the equipment. Four stations, DTB3, DTB4, DTB5, and DTB9, were located outside the plot and were used to control the spatial characteristics of the entire deformation area. In the future, these stations could be used to invert the deformation of the entire region without using InSAR data.

Generally, the IDAL method can be used to determine the number and spatial distribution of GNSS stations in a mining area experiencing deformation, considering accuracy, cost, and specific facilities. These GNSS stations captured the spatial characteristics of the deformation to the greatest extent possible, and the proposed method was used to improve the locations of GNSS stations and to provide guidance for station layout in GNSS deformation monitoring networks.

4. Discussion

4.1. Effect of initial GNSS station distribution on deformation capture ability

To explore the effect of the initial GNSS station distribution on the ability of the IDAL method to capture deformation, the simulated data and the layout of nine stations were considered as an example. The random layout (RL), uniform layout (UL), and empirical layout (EL) methods were used for the initial stations, which was followed by the application of the IDAL method to obtain a solution. The processes for solving the locations of these three different initial stations using on the IDAL method are named RL-IDAL, UL-IDAL, and EL-IDAL, respectively.

Figure 11(a–c) show the deformation recovered using the RL-IDAL, UL-IDAL, and EL-IDAL methods, respectively: the distribution of the initial stations widely varies among the methods. However, these three completely different initial station layouts can recover the spatial characteristics of the deformation after being optimized using the IDAL method. The distribution of the red triangle stations shows that the positions of the optimized final stations are not considerably different, which preliminarily shows that the IDAL method is less dependent on the initial locations of the stations. Figure 11(d–f) present the differences between Figure 11(a–c) and the simulated deformation, respectively. After the initial stations of these three methods were optimized using the IDAL method, the error distributions of the recovered deformation were similar, and the magnitudes were relatively small. The error distribution histogram in Figure 11(g–i), the error of the RL-IDAL method was slightly higher, the error of the EL-IDAL method was the smallest, and the errors of these three methods were close to a normal distribution.

Figure 11. The effect of different initial station locations on the results produced with the IDAL method (black triangles: fixed stations; white triangles: initial stations to be optimized; red triangles: stations optimized using IDAL method).

Table 3 lists the MAE, STD, and RMSE values for the three methods. The EL-IDAL method was superior to the RL-IDAL and UL-IDAL methods for all indices. This is because the EL-IDAL method considers deformation, which helps guide the reasonable distribution of initial stations; thus, the optimized deformation accuracy obtained using this method was the highest.

Table 3. Statistical error analysis of RL-IDAL, UL-IDAL and EL-IDAL methods in Figure 11(g–i).

Method	MAE (cm)	STD (cm)	RMSE (cm)
RL-IDAL	0.67	0.82	0.84
UL-IDAL	0.61	0.73	0.75
EL-IDAL	0.59	0.73	0.72

The results of the above analysis show that the different initial station distributions weakly impact the final optimization results, showing that the IDAL method is relatively robust. However, from subtle differences, it can still be observed that the more reasonable the distribution of the initial stations, the higher the final optimization accuracy. This means that before determining station layout, a somewhat reliable initial station distribution can be provided according to the actual requirements and InSAR deformation results, which can help to further strengthen the spatial deformation monitoring capability of GNSS.

4.2. Effect of search thresholds on the efficiency of the IDAL method

Figure 12 shows the variation in the RMSE of the deformation on two different scales for different numbers of iterations during the optimization process. The RMSE relatively quickly changes when the initial station is optimized on large-scale CPs, and the rough solution quickly converges. When the rough solution is optimized on small-scale CPs, the change in the RMSE is small; however, it also quickly converges to obtain a refined solution. This confirms the validity of the theoretical method developed in this study.

Figure 12. When the number of stations (Ns) is 9, the RMSE of the stations obtained with the IDAL method varies with the number of iterations.

However, the distance threshold must be set for the iterative solution because the distance threshold at different scales affects the efficiency of the solution. Therefore, to study the rule for selecting the distance threshold at different scales, the simulated data and the layout of nine stations were again taken as examples to design an experiment of the influence of the combination of different distance thresholds on the accuracy and efficiency of the solution.

As shown in Figures 13(a) and (b), as long as $T_1$ and $T_2$ are sufficiently large, the optimization accuracy is not affected; however, the solution efficiency is seriously impacted. Moreover, $T_1$ and $T_2$ cannot be too small; otherwise, the CPs in some regions cannot be searched, resulting in low accuracy.

Figure 13. Effect of distance search threshold on different scales on IDAL (a) efficiency and (b) accuracy.

The results of the above experiment demonstrate that selecting the appropriate parameters can increase the efficiency by more than two times without affecting the accuracy, and the improvement is more noticeable when the study area is larger. Here, the selection of $T_1$ can be determined based on the initial stations. If the initial stations are reliable (such as those obtained using the EL method), a smaller $T_1$ value could be used. $T_2$ can be smaller because, after large-scale optimization, the positions of the stations are relatively accurate. To be safe, a value of 200–400 m can be chosen.

5. Conclusions

In this study, a GNSS observation network layout and an optimization method based on InSAR prior deformation information were developed. For an area experience ground deformation, the user only needs to provide a quality standard for deformation monitoring to rapidly obtain the number, spatial distribution, and corresponding deformation recovery accuracy of the GNSS station layout. The main findings are as follows:

1. The IDAL method fully uses of the advantages provided by InSAR technology in deformation monitoring, such as its large scale, high precision, and high spatial resolution. Considering factors such as spatial deformation monitoring accuracy, cost, topography, and the actual monitoring requirements of GNSS stations, according to the deformation spatial correlation, a multiscale iterative optimization method was adopted to take advantage of the ability of GNSS stations to capture the spatial characteristics of regional group deformation and increase the availability of GNSS deformation monitoring results.

2. Simulation experiments demonstrated that, compared with the traditional station layout methods that scarcely considered deformation distribution, the new method improved GNSS spatial deformation monitoring ability by 46%–68%. The algorithm exhibited high stability, efficient performance, and strong applicability. Utilizing ALOS2 data and the SBAS–InSAR method, the deformation of a mining area in Datong City, Shanxi Province, China, was assessed. The findings revealed that deployed GNSS stations could effectively capture the deformation of the mining area with local subsidence. In comparison with the InSAR deformation results, 89.6% of the deformation recovered by the new method exhibited a point error of less than 1.5 cm. In situations where InSAR data is unavailable, this method can be used as a substitute to promptly understand the safety and deformation status of photovoltaic equipment in the region. The proposed method provides guidance and a reliable basis for determining the station layout of GNSS deformation monitoring networks and provides a suitable supplement for the current GNSS station layout, which is conducive to the monitoring and understanding of geological disasters given the frequency of geological disasters.

3. In the future, this method will need to be extended to different ground deformation scenarios, such as landslides, earthquakes, and groundwater as InSAR technology matures and given the widespread availability of free SAR data, which reduces the difficulty of acquiring deformation information. However, due to the high dependence on the InSAR deformation results, the method will fail in some completely decoherent regions. Additionally, in this study, only InSAR LOS deformation data were considered. However, surface deformation occurs in the three-dimensional direction. In the future, attention will be paid to incorporating three-dimensional deformation results into the IDAL model.

Acknowledgements

The authors wish to thank the Japan Aerospace Exploration Agency (JAXA) for providing the ALOS-2/PALSAR-2 data.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The ALOS-2/PALSAR-2 data used in this paper can be applied for on the https://auig2.jaxa.jp/ips/home. Other data in this paper can be obtained by contacting the authors.

Additional information

Funding

This work was partly supported by the National Natural Science Fund for Distinguished Young Scholars, grant number 41925016; the National Natural Science Foundation of China, grant number 42330112; and the Innovation Foundation for Postgraduate of Central South University under grant 2021zzts0264.