Data Acquisition of GNSS-Based InSAR: Joint Accuracy-Efficiency Optimization of 3-D Deformation Retrieval

In the Global Navigation Satellite System-based Synthetic Aperture Radar Interferometry (GNSS-based InSAR) system, the 3-D deformation retrieval accuracy and monitoring efficiency are very poor for the majority of the one repeat-pass period. It is necessary to select the data acquisition time and corresponding satellite combination to obtain better monitoring accuracy and efficiency. In this article, an experimental design algorithm for GNSS-based InSAR is proposed to achieve the highest monitoring efficiency with the restriction of the 3-D deformation retrieval accuracy. First, the joint optimization model is proposed, and based on an analysis of different satellite trajectories, it is proven that the optimal solution can be obtained. Second, the 3-D deformation retrieval accuracy model is derived based on the bistatic configuration, and the monitoring efficiency is evaluated based on the number of effective independent points of a single SAR image. The raw data are employed to indicate the effectiveness of the proposed joint accuracy-efficiency optimization algorithm in GNSS-based InSAR.

GNSS-based InSAR has a series of advantages. First, the repeat-pass period is much shorter than that of traditional Low-Earth-Orbit (LEO) SAR [16], [17], [18], [19], [20], [21], [22], Manuscript [23], [24], [25], [26]. Considering the Beidou system, which has a Medium Earth Orbit (MEO) satellite and Inclined GeoSynchronous Orbit (IGSO) satellite, as an example, their repeat-pass periods are 7 days and 1 day, respectively. Thus, the system is better for rapid deformation scenes, such as landslides and mining collapses. Second, there have been 10 on-orbit Beidou IGSO satellites and 27 on-orbit Beidou MEO satellites, which can guarantee multiangle deformation retrieval results for one target, and 3-D deformation detection can be achieved by combining these deformation results. Third, the cost of the GNSS-based system is low because it only requires a receiver to continuously complete 3-D deformation retrieval, i.e., there are no costs at the transmitter. Therefore, GNSS-based InSAR can be widely applied in 3-D deformation detection. Nevertheless, there are some problems in GNSS-based deformation retrieval experimental design. To achieve 3-D deformation retrieval, it is necessary to combine different single-angle measurement results of different satellites. Since many satellites are in orbit, there are always multiple satellites at any time that can illuminate the target area, which means that there are many choices for transmitters at any time. However, the accuracy and efficiency of 3-D deformation retrieval differ for different satellite combinations and experimental times. Thus, the experiment time and the corresponding satellite combination need to be taken into account in the experimental design to guarantee the best retrieval accuracy, highest monitoring efficiency, and shortest processing time. Specifically, the following factors need to be considered in the experimental design.
First, 2-D resolution needs to be considered. Along the range direction, based on the waveform and effective bandwidth, the resolution is approximately 10 m. Along the azimuth direction, limited by the data size, the accumulation time is 600 s, and the azimuth resolution after bistatic configuration optimization is approximately 5 m. Compared with traditional InSAR systems, the resolution of GNSS-based InSAR is low, which will blur the image. Moreover, the phase of the PS point can be affected by other weak targets in the resolution cell. The lower the resolution is, the greater the effect. Therefore, the resolution of GNSSbased InSAR needs to be restricted.
Second, in GNSS positioning, the positioning accuracy is related to the elevation angle. Similarly, in GNSS-based InSAR, the accuracy of the interferometric phase is also related to the elevation angle. The elevation angle also needs to be restricted. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Third, 3-D deformation retrieval is achieved by combining 1-D deformation retrieval results from different angles. Therefore, the 3-D deformation retrieval accuracy depends on the joint configuration of multiple transmitting satellites and the measurement accuracy of each satellite. Therefore, a quantitative evaluation method needs to be proposed to jointly optimize these two factors to obtain the best 3-D deformation retrieval accuracy.
Fourth, in permanent scatterer (PS) InSAR technology, to achieve high monitoring efficiency of the scene, the number of PS points in the scene needs to be large. However, in the experimental design stage, a long series of SAR images cannot be obtained, and PS points cannot be selected. Therefore, how to select the experiment time and satellite combination that correspond to the highest monitoring efficiency according to a single SAR image need to be solved. The experimental design for GNSS-based InSAR is a process of joint optimization of these factors. A corresponding experimental design algorithm needs to be proposed.
There have been many experiments based on the GNSS [27], [28], [29], [30], [31], [32], [33], [34], and the most widely employed for GNSS satellites is positioning. In [35], a geometry dilution of precision (GDOP) analysis was performed to obtain concise analytical expressions for several scenarios, which are generally applicable to geometries where the mobile device is surrounded by base stations. In [36], neural network-based, navigation satellite subset selection was presented. The approach is based on approximation or classification of the satellite GDOP factors that utilize the neural network approach. In [37], the approximation of GPS GDOP using statistics and machine learning methods was discussed. This study compares the performance of several methods, such as linear regression, pace regression, isotonic regression, support vector machines (SVMs), artificial neural networks, and genetic programming. The experimental results show that SVMs and genetic programming have better performance than other methods. In [38], a methodology in the context of convex geometry theory was proposed to reduce the number of satellites among all visible satellites utilized for the position computation by a navigation receiver. This reduction is based on the approximate equivalence between the geometric dilution of precision minimization and the maximization of the size of the polytope expanded by the user-to-satellite unit vector endpoints. In [39], a gradient-based control method where a position dilution of precision (PDOP) field is applied to indicate the performance of the positioning service was proposed. The gradient of the PDOP field is driven through rigorous derivation, which can be the distributivity computed by each robot using only the local information from itself and its neighbors. In [40], a unified framework was introduced to compare different system types based on geometric matrix and error modeling, and a performance analysis in terms of GDOP and positioning accuracy bounds was presented. However, for GNSS-based InSAR, the interferometric phase is employed for deformation retrieval. The signal processing method and the form of the error are different from positioning. Therefore, these algorithms cannot be utilized in GNSS-based InSAR systems.
For deformation monitoring experiments, commonly employed systems include ground-based SAR (GBSAR) and LEO SAR, and some public experimental results are obtained. For GBSAR, in [41], polarimetric GBSAR was employed for a one-year measuring campaign in the village of Sallent in northeastern Spain. The problem of extracting subsidence information from fully polarimetric SAR acquisitions for the retrieval of high-quality deformation maps is addressed. In [42], a new methodology for the nondestructive measurement of the absolute displacement of a pier during a bollard pull trial was proposed by GBSAR interferometry. This technique measures displacement patterns with submillimeter precision in any weather condition, and the operating distance can reach 4 km from the target area. In [43], conventional nonlocal methods developed for spaceborne SAR were investigated with GBSAR data for the first time. The proposed method was integrated into a complete GBSAR, small-baseline subset algorithm.
For LEO SAR, in [44], to investigate the impact of ground subsidence on the railway line, 26 TerraSAR-X SpotLight images, 23 TerraSAR-X StripMap images, and 31 Sentinel-1 Interferometric Wide images were processed to obtain the deformation of the Yizhuang area. Through the vertical and horizontal deformation of the railway using the fusion of the railway deformation results from different lines of sight directions, it is found that the deformation of the railway line mainly occurs in the vertical direction. In [45], the structure and deformation features of the Qinghai-Tibet Railway were retrieved and analyzed using time-series interferometry with Sentinel-1 A and TerraSAR-X images. The backscattering and coherence features of the Qinghai-Tibet Railway are analyzed in medium-and very high-resolution SAR images. In [46], a study that employs the corner reflector InSAR technique to monitor the deformations of ten corner reflectors installed at the Xiangtan converter station in China was reported, with seven TerraSAR-X spotlight images. Two-step phase unwrapping tactics are proposed based on the least squares ambiguity decorrelation adjustment algorithm to focus on estimating the nonlinear deformation without being affected by unwrapping errors. In [47], a method based on the minimization of the L1-norm within local spatial cubes to reconstruct 3-D displacement vectors, which can be obtained from tomographic, SAR point clouds available from at least three different viewing geometries, was proposed. The methodology is applied to two pairs of cross-heading combinations of ascending and descending TerraSAR-X spotlight image stacks over the city of Berlin. However, for the GNSS-based InSAR system, first, the bistatic configuration makes it different from GBSAR and LEO SAR, and multiple on-orbit satellites render the experimental design very diverse and complex. Second, due to the low 2-D resolution and SNR, in the experimental design, both the image quality and the interferometric phase accuracy need to be considered. Traditional deformation retrieval methods do not need to consider these issues.
GNSS-based InSAR uses different on-orbit navigation satellites as transmitters to achieve 3-D deformation retrieval. Compared with other InSAR systems, the target area can be illuminated by many navigation satellites at any time, and the accuracy and efficiency of 3-D deformation retrieval that can be achieved by different satellite combinations at different times are different. Therefore, the experiment time and the satellite combination corresponding to the best monitoring accuracy and efficiency need to be selected. The experimental design for GNSS-based InSAR is a process of joint optimization, which requires consideration of the 2-D resolution, elevation angle, joint configuration, 1-D monitoring accuracy, and monitoring efficiency. Traditional InSAR systems do not have many transmitters available, and few factors need to be considered in the experimental design process. Therefore, an experimental design algorithm for GNSS-based InSAR needs to be proposed.
In this article, an experimental design algorithm for GNSSbased InSAR is proposed. The theoretical 3-D deformation retrieval accuracy and monitoring efficiency are modeled and analyzed. The rest of the article is organized as follows. Section II introduces the experimental design model. Section III analyzes the characteristics of multiangle satellites. Section IV proposes the experimental design algorithm. Section V shows the design result of the raw data. Finally, Section VI concludes the article.

II. GNSS-BASED INSAR SYSTEM EXPERIMENTAL DESIGN OPTIMIZATION MODEL
The bistatic configuration of GNSS-based InSAR which uses east-north-up coordinates is shown in Fig. 1. The transmitters adopt Beidou satellites. The elevation angle ς is defined as the angle between the horizontal plane and the upward direction is positive, while the azimuth angle ϑ is defined as the angle between the east direction and the counterclockwise direction is positive. The receiver is fixed at the origin of the coordinate. Two antennas receive the echo signal and direct signal. The system parameters are shown in Table I [48].
There are already more than 100 on-orbit navigation satellites. At any time, the GNSS can guarantee that many satellites illuminate the target area. In addition, 3-D deformation retrieval requires combining the measurement results of multiple satellites. However, the monitoring accuracy and efficiency of different satellite combinations at different times are different. Therefore, an optimization model is required to obtain the best experiment time and corresponding satellite combination.
For the deformation retrieval system, the highest monitoring efficiency is desired. For the PS-InSAR technology adopted in GNSS-based systems, the number of PS points represents the monitoring efficiency of the scene. However, the number of PS points available for different satellite combinations at different times is different. Therefore, the experiment time and satellite combination that correspond to the greatest number of PS points, i.e., highest monitoring efficiency, need to be obtained.
At the same time, there are some restrictions on the system. First, the 2-D resolution determines the quality of the image. Moreover, the phase of the PS point is affected by other weak targets in the resolution cell. The lower the resolution is, the larger the number of interfering targets and the lower the phase accuracy. Second, the elevation angle reflects the influence of the atmosphere on the phase, which also needs to be considered. Third, the theoretical 3-D deformation retrieval accuracy, which depends on the joint configuration and 1-D measurement accuracy, needs to meet the requirements of deformation monitoring experiments. All these restrictions need to be considered in the optimization model.
Combining the abovementioned analysis, monitoring efficiency is selected as the optimization target, and other factors are employed as the limiting conditions. An optimization model is proposed (1) where t is the experiment time, SC is the satellite combination, and ef f is the function of the monitoring efficiency. S reso , ς, and ACC are the resolution, elevation angle, and theoretical 3-D deformation accuracy, respectively. S reso_thre , ς thre , and ACC thre are the thresholds. There have been studies on the resolution S reso and elevation angle ς [49], [50]. Therefore, this article mainly discusses the monitoring efficiency and theoretical 3-D deformation retrieval accuracy.

III. ANALYSIS OF MULTIANGLE SATELLITES
For one navigation satellite, the corresponding bistatic configuration varies greatly within a repeat-pass period. Consider the Beidou2 IGSO3 satellite as an example to show the change in the bistatic configurations. The target scene is located in Changshu, Suzhou. In one repeat-pass period, the change in elevation and azimuth angle is shown in Fig. 2(a). The azimuth angle ranges from 124 • to 260 • , and the elevation angle ranges from −4 • to 76 • . The threshold of the elevation angle is set to 15 • . The change in resolution [50] is shown in Fig. 2(b), where the range of range resolution is 8.9-25.2 m, and the range of azimuth resolution is 3.24-18.8 m. The azimuth threshold is 8 m, and the range threshold is 10 m. The elevation angle should be greater than the threshold, and the resolution should be less than the threshold. The results indicate that both conditions cannot be satisfied for most of the repeat-pass period.
In addition, the monitoring efficiency is analyzed. From Fig. 2, it can be seen that the illumination angle of one satellite changes a lot in one repeat-pass period. The range cross section (RCS) of the scene varies with the illumination angle. It means that the number of the PS point and the corresponding monitoring efficiency also varies with the illumination angle. What's more, for different satellite combinations at different times, the effect of deformation monitoring will have greater changes. Therefore, an experimental design algorithm is needed to obtain the satellite combination and the experimental time that corresponds to the best accuracy and highest efficiency. In the experiment, based on the prior satellite trajectory, the resolution, elevation angle and theoretical monitoring accuracy can be calculated. Since the scattering characteristics of the PS points are constant, the number of PS points, i.e., monitoring efficiency, can be predicted. Therefore, the highest monitoring efficiency within the restriction can be achieved, i.e., (1) is solvable based on the information in the experimental design stage.

IV. GNSS-BASED INSAR EXPERIMENTAL DESIGN ALGORITHM
In this article, an experimental design algorithm is proposed for GNSS-based InSAR. First, based on the satellite trajectories and corresponding bistatic configurations, the restriction of the experiment time and satellite combinations can be obtained. Second, according to the single image result, the effective independent point (EIP) is introduced here to describe the monitoring efficiency. The flowchart is shown in Fig. 3.

A. Restriction of the Experimental Design
As discussed in Section II, the restriction of the experiment includes resolution, elevation angle, and theoretical 3-D deformation retrieval accuracy. Based on the bistatic configuration, the resolution and elevation angle can be obtained [49], [50]. The theoretical 3-D deformation retrieval accuracy is then discussed. As mentioned earlier, 3-D deformation retrieval accuracy depends on the joint configuration and 1-D measurement accuracy. This section discusses the optimization of these two factors.
PS-InSAR is employed in the GNSS-based InSAR system. Based on this, the interferometric phase model is established. To eliminate the synchronization error between the receiver and the transmitter, the echo phase is subtracted from the direct signal phase. For a PS point Q, the phase in the kth SAR image can be expressed as where t k is the central time of kth synthetic aperture time, P S is the position of the satellite, P Q is the position of the target, and P E is the position of the receiver. ϕ noise is the phase noise, and λ is the wavelength. Let t n be the reference time, and the interferometric phase can be expressed as From (3), the interferometric phase can be simplified as In (4), ϕ topo (t k , t n ; P Q ) is the topography phase and Δϕ defo (t k , t n ; P Q ) is the deformation phase The topography phase can be compensated for through digital elevation matrix (DEM) information, and the accuracy of the deformation retrieval is determined by Δϕ defo (t k , t n ; P Q ). According to (6), the relationship between the actual deformation and the interferometric phase can be expressed as follows: where Φ is the observation result of M satellites, D is the actual deformation of the target, and n is the noise within the observation of M satellites.
The variance of the deformation D can be obtained as From (9), the optimization of the joint configuration can be achieved by H. Next, the effect of the 1-D measurement accuracy is discussed. Similarly, based on (9), the measurement accuracy is reflected by E(nn T ). n needs to be expressed based on the bistatic configuration.
For the deformation retrieval system, only the interferometric phase is utilized. The accuracy of the interferometric phase is related to the coherence coefficient, which can be expressed as where γ tem is the temporal coherence coefficient, γ noi is the noise coherence coefficient, and γ spa is the spatial coherence coefficient. In PS processing, γ tem is approximately equal to 1. In addition, the noise coherence coefficient is also approximately equal to 1 when the SNR is 25 dB [51]. Therefore, the total coherence coefficient is approximately equal to the spatial coherence coefficient, i.e., γ ≈ γ spa . Different from the traditional PS-InSAR technique, as mentioned in Section I, the resolution of GNSS-based InSAR is low, approximately 10 m × 5 m (r × a). The phase of the resolution cell is contributed by the PS point and other weak targets in the resolution cell around the PS point. The spatial coherence of the PS point can be disregarded, but the spatial coherence of other targets must be considered [15]. The lower the spatial coherence coefficient is, the greater the influence of the weak targets on the interferometric phase of the PS point. Therefore, the spatial coherence coefficient of the scene can be applied as a basis for calculating the interferometric phase noise. The expression of the spatial coherence coefficient is whereW (x, y) is the normalized point spread function (PSF). x and y represent the target position in the imaging plane. r(x, y; t) is the distance from the satellite to the target at time t. The measurement accuracy can then be expressed as where σ ϕ is the standard deviation of the interferometric phase and σ w is the standard deviation of deformation detection. The monitoring accuracy of each angle, i.e., E(nn T ), can be expressed as Next, the 3-D deformation retrieval accuracy can be calculated. Equation (9) can be rewritten as The total accuracy is defined as The range of the spatial coherence coefficient of Beidou IGSO satellites in one repeat-pass period is calculated; the results are shown in Fig. 4. Fig. 4(a) shows the result of the spatial coherence coefficient, and Fig. 4(b) shows the result of the corresponding measurement accuracy. It can be found that some satellites have high spatial coherence coefficients most of the time, while the spatial coherence coefficients of some satellites vary greatly. This is because different satellite trajectories have different angle changes relative to the monitoring scene. In 3-D deformation retrieval processing, the total accuracy is reduced when 1-D measurements with low spatial coherence coefficients are introduced. For different target positions, the selected satellite combinations are different. Therefore, the optimization of the satellite combination is necessary. Thus far, combining the results of the resolution, elevation angle, and theoretical 3-D deformation retrieval accuracy, several sets of experimental times and satellite combinations that meet the requirements can be obtained. Let T time points {t 1 · · · t T } that correspond to T satellite combinations {SC 1 · · · SC T } meet the requirements. Next, the combination corresponding to the highest monitoring efficiency needs to be further selected.

B. Monitoring Efficiency
In PS processing, the number of PS points represents the monitoring efficiency. However, in the experimental design stage, long time series images cannot be obtained, so the PS points cannot be extracted. The concepts of independent point (IP) and EIP are introduced here to replace the PS point to describe the monitoring efficiency. For the PS points, its RCS should be higher than other targets to ensure that its phase is not affected by other targets and noise. The amplitude and phase of the resolution cell should almost depend on the PS point, and its resolution cell should be similar to the theoretical PSF. A point needs a complete PSF to become a PS point. These points, which have complete resolution cells, are referred to as IPs.
The selection of IPs can be based on the coherence between the theoretical resolution cell and the actual resolution cell. For the theoretical resolution cell, the ambiguity function between the target vector Q and its neighboring target vector B can be expressed as where Φ T Q and Φ RQ are the unit vectors from the transmitter and receiver to the target Q, respectively. β is the bistatic angle. Θ is the direction along the bisector of β. ω E and Ξ are the equivalent angular speed and equivalent motion direction. p is the result of range pulse compression. m Q is the result of azimuth pulse compression. λ is the wavelength, and c is the speed of light. The theoretical resolution cell of target Q is defined as Next, the actual resolution cell is discussed. In GNSS-based InSAR systems, backprojection algorithm is used to get the SAR images [14], [52]. The actual resolution cell refers to the resolution cell extracted from the real image. In the real image results, the resolution cell of one target is influenced by noise and other weak targets, so it is different from the theoretical resolution cell. However, for PS points, the extracted resolution cell should be similar to the theoretical resolution cell. The extraction method is, first, find target Q when Q is the local maximum point in S(Q) P SF . Second, the actual resolution cell can be expressed as

S(Q) ACT = {I(B)|I(B) > I(Q) − 3 dB} (18)
where I is the real SAR image amplitude. Let where cohe is the coherence coefficient. The points that have high cof f are considered IPs. IPs that are subject to a large amount of interference need to be removed. For IPs, there are two main sources of phase error: 1) noise and 2) the sidelobe of the adjacent point. The diagram is shown in Fig. 5. The blue point is the target, and the red point is the strong adjacent point, whose sidelobe will have an impact on other targets. The orange point is noise. When the interference and noise cannot affect the phase of the target, this IP can be utilized as an alternative to the PS point, which is referred to as EIP. Since the number of PS points cannot be obtained, it is suggested that the larger the number of EIPs, the larger the number of PS points, and the higher the monitoring efficiency.
EIPs are selected based on the results of IPs. The interference of noise can be expressed through the SNR. The 5-mm deformation retrieval accuracy corresponds to a 17-dB SNR. IPs with insufficient SNR will be eliminated.
The influence of the adjacent points is analyzed. For target Q, the influence of adjacent point A can be expressed as where σ A is the amplitude of the SAR image. The total interference value of L adjacent points to Q is If Q can satisfy (22) where Amp thre is the threshold, which is also 17 dB corresponding to the 5 mm accuracy of the deformation retrieval; Q is considered the EIP. The experiment time and corresponding satellite combination can be obtained. Consider t 1 in {t 1 · · · t T } as an example. The number of EIPs are {M 1 · · · M N }, which correspond to satellite {S 1 · · · S N }. The average number of EIPs for t 1 is The satellite combination that corresponds to t can be obtained according to the results in Section IV-A.

V. RAW DATA PROCESSING
In this section, a GNSS-based experiment is conducted. Based on the proposed algorithm, the best experiment time and satellite combination are obtained.

A. Experimental Scene
The experimental scene, which is located in Changshu city, Jiangsu Province, China (31.7582N, 120.9323E), is shown in Fig. 6. The scene includes a lake, a road surrounding the lake, and a highway farther away. The remainder of the scene is grass, whose RCS is low. The receiver contains two independent channels, including an omnidirectional antenna for receiving the direct signal and a horn antenna for receiving the scene echoes. The receiver is located on the roof near the lake, and the grazing angle is around 2.7 • . The two antennas are shown in Fig. 7. The left antenna is the scene echo antenna, and the right antenna is the direct signal antenna.

B. GNSS-Based InSAR Experimental Design
The experiment started on 2021.05.03. The number of selection satellites was 10, including 7 Beidou2 IGSO satellites and 3 Beidou3 IGSO satellites. The repeat-pass time of the IGSO satellite is 1 d, which is divided into 48 periods, from which the optimal acquisition time is selected.
Based on the trajectories of the satellites, the 2-D resolution and elevation angle are analyzed. Selecting Beidou3 IGSO1 as an example, the resolution cell area and elevation angle in one repeat-pass period are calculated, which is shown in Fig. 8. Fig. 8(a) shows the results of the resolution cell area, and Fig. 8(b) shows the results of the elevation angle. The two factors vary greatly during one repeat-pass period.
The restriction of the resolution cell and the elevation angle are calculated first. The area of the resolution cell is required to be less than 80 m 2 , and the elevation angle is required to exceed 15 • . Satellites that meet these two conditions are considered available in this experiment. The same process is performed on all 10 IGSO satellites, and the period that satisfies the conditions is selected. The result is shown in Fig. 9, where yellow means that the satellite meets resolution and elevation angle requirements, and blue means that the satellite is not available at that moment. To achieve 3-D deformation retrieval, the required number of available satellites is at least three so that the periods that contain at least three available satellites are selected. Fig. 9. Selection results of 10 Beidou IGSO satellites. Yellow means that the satellite is available at that moment, and blue means that the satellite is not available at that moment.  Then, the theoretical 3-D deformation retrieval accuracy is calculated in each period. First, the spatial coherence coefficients that correspond to each satellite in the available periods are shown in Fig. 10, which are served as the error input of ACC. For the selected periods, calculate the optimal ACC, and the result is shown in Fig. 11. The upper limit shown in the figure is set to 10 mm. Based on the ACC results, there are three time points whose ACC is less than 5 mm, and the corresponding selected satellites are shown in Table II. To date, the restriction   of the experimental design has been discussed, and the final experimental time and satellite combination will be determined based on the monitoring efficiency.
The monitoring efficiency is optimized based on the number of EIPs. Considering the first time point 2:00 as an example, the selected satellites are Beidou2 IGSO1, Beidou2 IGSO4, and Beidou2 IGSO6. The 3-dB resolution cell S(Q) P SF of each satellite is shown in Fig. 12. Fig. 12(a), (b), and (c) are the results of Beidou2 IGSO1, Beidou2 IGSO4, and Beidou2 IGSO6, respectively.
The IPs are then selected. The selection results of Beidou2 IGSO1, which are shown in Fig. 13, are applied as an example. Both the thresholds of the SNR and adjacent point interference are 17 dB. For the local maximum points, the distribution of the coherence coefficient is shown in Fig. 13(a). There are 320 points with a coherence coefficient above 0.8, which can be considered IP points. The distribution of the IPs is shown in Fig. 13(b). Next, the EIPs are selected. Fig. 14(a) shows the    selected results based on the SNR threshold. The distribution of strong points is the same as the scene information. Fig. 14(b) is the enlarged result of the black box. The marked point is near the adjacent strong point, and its amplitude is much smaller than that of the adjacent strong point. Therefore, the phase of the marked point is unreliable, and the point cannot be implemented as an EIP. Fig. 14(c) is the result of removing the points with strong adjacent interference, and the enlarged results are shown in Fig. 14(d). The marked point in Fig. 14(b) is eliminated based on the proposed algorithm. The points in Fig. 14(c) are the EIP selection results, which can represent the monitoring efficiency.
The EIP selection results of Beidou2 IGSO4 and Beidou2 IGSO6 are shown in Fig. 15. Their distribution is similar to that of the scene information. The numbers of IPs and EIPs are shown in Table III. Through the proposed algorithm, only 35% of the IPs are implemented as EIPs.
According to the previous analysis, there are three periods to meet the requirements of deformation retrieval accuracy. The same algorithm is used to obtain the number of EIPs of different SAR images, and the results are shown in Table IV. According to the results, at 2:00, the number of average EIPs is the largest, which means that it corresponds to the largest monitoring efficiency. The satellites employed for transmitting images are Beidou2 IGSO1, Beidou2 IGSO4, and Beidou2 IGSO6. To date, the experimental time and corresponding satellite combination have been determined.  From the entire process, the experimental design is necessary for a GNSS-based InSAR system. First, in one repeat-pass period, the resolution and elevation angle change greatly, so these two factors must be selected to ensure the image quality and phase accuracy. Second, according to the theoretical 3-D deformation retrieval accuracy, in all 48 periods, only three periods meet the accuracy requirements. These factors must be considered to limit the experimental design.
Based on the SAR images from different satellites, EIPs are used to analyze the monitoring efficiency. Because the scattering characteristics of the scene under different angles vary greatly, the number of strong scattering points, i.e., the number of EIPs, varies greatly. From the results in Table IV, it can be seen that the number of EIPs for different satellites changes a lot. So, EIP analysis based on SAR images is necessary in GNSS-based InSAR system.

C. Comparison Results
To indicate the effectiveness of the proposed algorithm, 22 days of data were collected at different periods for deformation retrieval. The raw data were collected twice in a repeat-pass period. For the first set of data, according to the experimental design results, the collection time point was 2:00 on 2021.05.03. The second set of data is employed as the comparison result, and the collection time point is 14:00 on 2021.05.03. The experimental design results at the two time points are shown in Table V.
Based on the algorithm proposed in [53], the PS points can be selected. The algorithm proposed in [15] is selected to obtain the deformation retrieval results. The number of PS points and 3-D deformation retrieval accuracy are shown in Table VI. It can be seen from the experimental results that the deformation retrieval accuracy at 2:00 is better than that at 14:00. The number of PS points is also greater than that at 14:00. The results indicate that the analysis of the theoretical deformation retrieval accuracy and EIPs are effective in GNSS-based InSAR systems.
However, for the 3-D deformation retrieval results, the actual accuracy and theoretical accuracy are quite different, because the interferometric phase error, including the channel phase error and atmospheric phase error, cannot be fully compensated. What's more, the scene DEM error will also introduce a deviation in the interferometric phase. These factors cause a reduction in the deformation retrieval accuracy. However, under the same processing flow, the 3-D deformation retrieval accuracy at 2:00 is better than that at 14:00, which indicates that the theoretical deformation retrieval can represent the actual deformation accuracy.
For the number of PS points, the final PS point number is less than the number of EIP points because some EIPs are affected by time decoherence and spatial decoherence, and the scattering characteristics change during the experiment time. However, the experimental results indicate that when the number of EIP points is large, the number of final PS points is also large. The result shows that analyzing the number of PS points, i.e., monitoring efficiency through EIP points is effective. The results of the comparative experiments demonstrate the effectiveness of the proposed experimental design algorithm in GNSS-based InSAR systems.

D. Additional Result
To further prove the difference of multiangle image results, three Beidou satellites are used to image the mining area, and the results are shown in Fig. 16. The azimuth angle and elevation angle are (176 • , 45 • ), (141 • , 33 • ), and (152 • , 39 • ), respectively. From the results, two phenomena can be observed.
First, the image intensity is uneven in one satellite image result. As shown in Fig. 16, the intensity in the red box is much larger than that in the green box. Based on the bistatic configuration, the average angles of three satellites between the bistatic angle direction and the scene normal vector direction can be calculated, which are 18.92 • and 86.52 • of the red and green boxes, respectively. It means that the red box is closer to the specular reflection than the green box, so the intensity of the former is higher.
Second, for the SAR images of different satellites, the image intensity also varies a lot. To explain this phenomenon, the range of the bistatic angle are calculated, which are [37.05, 59.80], [21.69, 35.13], and [24.66, 39.42], respectively. It is suggested that the change of the image intensity is corresponding to the range of the bistatic angle. Moreover, Bragg Scattering phenomena can also cause this difference [54], which requires further study. These imaging results further demonstrate that the imaging results from different angles vary a lot. Therefore, the proposed algorithm is necessary in GNSS-based InSAR systems.

VI. CONCLUSION
In this article, an experimental design algorithm is proposed for GNSS-based InSAR to achieve the highest monitoring efficiency with the restriction of resolution, elevation angle, and theoretical 3-D deformation retrieval accuracy. First, the joint optimization model is proposed. Second, the trajectories are analyzed to indicate that the optimal solution can be obtained. Finally, the theoretical 3-D deformation retrieval accuracy is calculated based on the bistatic configuration, and the monitoring efficiency is evaluated based on the EIPs of the signal SAR image. The raw data are applied to indicate the effectiveness of the proposed algorithm. Based on the design results, the resolution, elevation angle, theoretical 3-D deformation retrieval accuracy, and monitoring efficiency change greatly in one repeat-pass period. Therefore, the proposed experimental design algorithm is necessary in GNSS-based InSAR systems.