GNSS-based InBSAR system PS selection and registration method based on generalization area

Abstract

The invention discloses a new PS selection and registration method based on a generalization region for a GNSS-based InBSAR system. The method comprises the specific steps of calculating a scaling factor based on DEM precision, constructing a theoretical resolution unit according to system parameters so as to obtain the generalization region, screening a preselected PS based on coherence, and superposing the generalization region of the preselected PS in multiple days so as to obtain a preselected PS correlation sequence. And the sequence is subjected to repetition removal processing through same pixel removal and adjacent coherence criterion, and the final PS is determined based on amplitude deviation, so that high-precision and rapid PS point selection and registration are realized.

Description

Technical Field

[0001] This invention belongs to the field of bistatic synthetic aperture radar technology, specifically a PS selection and registration algorithm for a GNSS-based InBSAR system based on a generalized region. Background Technology

[0002] GNSS-based InBSAR (Global Navigation Satellite System Based Bistatic Synthetic Aperture Radar Interferometry) uses on-orbit navigation satellites as transmitters and receivers fixed on the ground to receive echoes from the target scene. Repeated orbit interferometry is employed to monitor surface deformation in the target area. Compared to traditional InSAR, this technique features shorter repetition periods, lower cost, and three-dimensional deformation detection. Furthermore, given the large number of navigation satellites (any scene can be illuminated simultaneously by eight or more BeiDou satellites), this system only requires the deployment of receivers to achieve deformation monitoring. Therefore, GNSS-based InBSAR offers advantages such as near real-time performance, high accuracy, and low cost, making it widely applicable in deformation monitoring.

[0003] GNSS-based InBSAR uses BeiDou navigation satellite signals, which are not specifically designed for InSAR deformation measurement. These signals have limited bandwidth and low transmission power, resulting in poor SAR image resolution, low signal-to-noise ratio, and increased susceptibility to interference from nearby targets. PS-InSAR is used in this system. In long-term SAR interferometric image sequences, the phase of the PS signal is relatively stable. Utilizing its phase characteristics can avoid interference from other phases and obtain high-precision deformation information. However, trajectory offsets and elevation errors of navigation satellites can lead to random focal position errors, making traditional PS selection and registration methods unsuitable. Furthermore, obtaining three-dimensional deformation requires processing deformation monitoring data from multiple satellites (at least three) on a single receiver, and traditional PS processing algorithms are time-consuming. Therefore, a fast and accurate PS selection and registration algorithm is needed. Summary of the Invention

[0004] In view of the above, this invention provides a novel PS selection and registration method for GNSS-based InBSAR systems based on generalized regions. This method replaces the single pixel in traditional PS selection methods with a generalized region, superimposing the generalized regions to obtain a sequence of pre-selected PS points. The final PS points and interferometric phase are then determined based on the amplitude distribution, achieving high-precision and fast PS selection and registration.

[0005] The technical solution of this invention is: a PS selection and registration method for a GNSS-based InBSAR system based on a generalized region, comprising:

[0006] GNSS-based InBSAR system configuration as follows Figure 1 As shown, an in-orbit navigation satellite serves as the transmitter, and a receiver fixed on the Earth's surface receives the target's reflected signal. The offset of the satellite's repeating channel (satellite baseline and positioning error vector sum) D and the DEM error Δh cause the focus position to shift. The following uses an East-North-Up coordinate system, with the receiver's position as the origin, to establish the signal model.

[0007] Under the conditions of random orbital error and elevation angle error, a phase model for target Q is established:

[0008]

[0009] In the formula, σ is the scattering coefficient, and P T For the transmitter position, P Q P is the location of the target. E Let Q be the location of the receiver. At non-target locations, the ambiguity function is used to represent the ambiguity between target vector Q and its neighboring target vector A. This ambiguity function can be expressed as:

[0010]

[0011] In the formula Φ TQ and Φ RQ Let θ be the unit vector from the transmitter and receiver to the target Q, respectively; β be the bistatic angle; Θ be the direction vector along the bisector of β; and ω be the angle bisector. E Ξ represents the equivalent angular velocity and equivalent motion direction, p represents the result of range pulse compression, and m represents the distance pulse compression. A This is the result of azimuth pulse compression, where λ is the wavelength and c is the speed of light. The target is displayed on the image in resolution units.

[0012] However, due to DEM errors, satellite positioning errors, and the influence of satellite spatial baselines, resolution cells exhibit random offsets. The projection offset, related to DEM errors and satellite positions, can be expressed as:

[0013]

[0014] In the formula () xy For the horizontal component, () z For the vertical component, V S For satellite velocity, G Q′ It can be represented as:

[0015]

[0016] A method for PS selection and registration in a GNSS-based InBSAR system based on generalized regions, flowchart as follows: Figure 2As shown, firstly, a generalization region is established to identify the associated sequences of pre-selected PS points. Then, duplicate pixels in the associated sequences of PS pre-selected points are removed, and non-PS points are removed using a coherence threshold. Finally, phase is extracted based on distance weights and phase weights.

[0017] Step 1: Select PS pre-selection points

[0018] Using radar SAR images obtained through the BP imaging algorithm on a GNSS-based InBSAR system, theoretical resolution cells are determined based on radar configuration and experimental parameters. Points with coherence coefficients greater than a threshold are selected as pre-selected PS points.

[0019] Perform PS pre-selection, and define the theoretical resolution cell area of ​​target Q as:

[0020] S PSF (Q)={A|χ(Q,A)>-3dB} (5)

[0021] In the formula, A represents the pixels surrounding Q, and the actual resolution unit area of ​​the target Q can be expressed as:

[0022] S(A) ACT ={A|I(A)>I(Q)-3dB} (6)

[0023] Where I is the amplitude of the SAR image. The correlation coefficient between the PS resolution unit extracted from the actual image and the theoretical resolution unit is:

[0024] coff(Q) = cohe(S(Q)) PSF ,S(Q) ACT (7)

[0025] Where cohe is the coherence coefficient of the two images, and points with sufficiently high coherence coefficients are selected as PS pre-selected points.

[0026] Step 2: Calculate the generalization area

[0027] To facilitate the association of pre-selected PS points, this invention employs a generalized area method. The concept of a generalized region is introduced, which replaces the pre-selected PS points.

[0028] In GNSS-based InBSAR, the SAR image is projected onto the imaging plane, and it is assumed that the phase of the PS is consistent with the phase within its generalization region. Through ambiguity function derivation, the expression for establishing the generalization region based on the theoretical resolution unit is:

[0029] S GEN (Q)={A|(σ GEN (AQ)+Q)∈S PSF (Q)} (8)

[0030] Where, σ GEN The scaling factor is determined based on the range of variation of the DEM across the entire scene, and is usually expressed as the variance of the DEM across the entire scene.

[0031] Overlapping the generalized regions of SAR images across multiple reorbiting cycles. For imaging results of a satellite over K consecutive days, find the pre-selected PS points in the daily imaging results, and denote the set of pre-selected PS points for each day as [Y1, Y2, ..., Y]. K This is used for subsequent PS selection. Taking the first image as an example, for the PS pre-selected point set Y1, the generalization area S of the entire image is... YGEN It can be represented as:

[0032]

[0033] Step 3: Extract the association sequence of PS pre-selected points

[0034] The generalized regions of PS pre-selected points over multiple days are overlaid, and the associated sequences of the pre-selected points in the image are extracted. Non-PS sequences with a small number of valid days are then removed.

[0035] If S YGEN A value of 1 indicates that the pixel can be effectively measured. The S values ​​of all images are then superimposed. YGEN Find the pixels that are continuously measured over time. Determine the association sequence of the PS preselected points based on the membership relationship between each pixel and the PS preselected points.

[0036] For the imaging results over K days, the expression after generalizing the region overlay is:

[0037] S GEN_ALL =S YGEN (Y1)+S YGEN (Y2)+…S YGEN (Y K (10)

[0038] For S GEN_ALL Threshold filtering is performed on the value of each pixel:

[0039] F GEN ={(x,y)|G GEN_ALL(x,y) ≥thre GEN} (11)

[0040] Where F GEN For the selected pixels, then GEN This represents the minimum threshold for the number of days a pixel can be detected. The association sequence of PS pre-selected points is determined based on the membership relationship between each pixel and the PS pre-selected points, as illustrated in the diagram below. Figure 3 As shown.

[0041] Introducing the coefficient ζ GENThis compensates for the direction error caused by the large resolution difference between range and azimuth in GNSS-based InBSAR.

[0042]

[0043] In the formula, R1 is the distance between the registered pixel and the PS pre-selected point, and R2 is the length of the generalization region in the R1 direction. When a pixel belongs to two different PS pre-selected points, its corresponding ζ is calculated separately. GEN Choose ζ GEN Smaller PS pre-selection points.

[0044] Step 4: Remove duplicate pixels

[0045] Since each resolution unit contains multiple pixels, there are many duplicates in the results. The PS pre-selected point association sequences are deduplicated, retaining as many longer sequences as possible. Then, a sequential coherence coefficient filtering method is used to retain sequences with higher coherence coefficients.

[0046] The PS pre-selected point association sequences in the above processing have a large number of repetitions. To ensure that the phase of each PS can only be used once, longer sequences are selected first for deduplication processing. The operation process is as follows: Figure 4 As shown.

[0047] 1) Arrange the pre-selected PS points in descending order of sequence length based on the previously obtained PS pre-selected point association sequence.

[0048] 2) Remove elements from the subsequent sequences that are identical to those in the current sequence to perform deduplication filtering. Then, rearrange the filtered results in descending order of sequence length.

[0049] 3) Repeat filtering and sorting until all duplicate pixels are removed.

[0050] This method preserves the integrity of long sequences as much as possible. A sequential coherence filtering method is employed, further filtering each sequence based on resolution unit coherence. Assume a sequence has L pairs of adjacent pre-selected PS points, each pair representing Q1 and Q... l+1 Calculate the coherence coefficient using the resolution unit:

[0051] coff(Q1,Q l+1 )=cohe(S(Q1) PSF ,S(Q l+1 ) PSF (13)

[0052] The threshold γ for the inter-resolution coherence coefficient is calculated based on the theoretical coherence coefficient γ as follows:

[0053] γ=γtem ·γ noi ·γ spa (14)

[0054] In the formula γ tem γ is the time coherence coefficient. noi γ is the noise coherence coefficient. spa Let be the spatial coherence coefficient. The scattering characteristics of PS do not change over a short period of time, meaning the temporal coherence coefficient is approximately 1. The relationship between the noise coherence coefficient and the signal-to-noise ratio can be expressed as:

[0055]

[0056] When the coherence coefficient of a set of pre-selected PS points exceeds T C At the same time, subsequent processing is performed to make the interferometric phase monitoring results of the same target more continuous.

[0057] Step 5: Select PS

[0058] After obtaining the correlation sequence of pre-selected PS points, the PS is determined using the amplitude deviation threshold method. The interferometric phase is extracted from the PS obtained by this algorithm and applied to the next step of high-precision deformation inversion.

[0059] For a PS preselected point association sequence of length N, according to equation (6), the peak amplitude value within each PS preselected point resolution cell is extracted as S(Q)=[Q1,Q2,…Q N Then, calculate its mean m. A =mean(S(Q)) and standard deviation σ Q = std(S(Q)). The expression for amplitude deviation is:

[0060]

[0061] Calculate the magnitude deviation of the associated sequence for each PS preselected point and compare it with the threshold. D If the comparison is greater than the threshold requirement, the associated sequence of the PS preselected point is considered to be PS.

[0062] The phase of the resolution cell containing the PS is extracted to reduce the impact of noise on the accuracy of the interferometric phase. Considering that for each PS, its effective phase extraction range is the intersection of its theoretical and actual resolution cells, according to equations (5) and (6), the phase extraction range S... PHA It can be represented as:

[0063] S PHA (Q1)=S PSF (Q1)∩S ACT (Q1) (17)

[0064] The phase of Q can then be expressed as:

[0065]

[0066] Where M is S PHA The number of pixels, W is the phase weight. p This represents the distance weight. The phase weight and distance weight can be expressed as:

[0067]

[0068] In the formula, P(m) represents the position of the m-th pixel. and σ P For S PHA The phase variance and distance from the peak are calculated. By using a weighted average, the influence of noise can be minimized, thus improving phase accuracy. The resulting PS will be used for multi-angle PS correlation, utilizing interferometric phase to achieve high-precision three-dimensional deformation inversion.

[0069] The present invention has the following beneficial effects:

[0070] 1. Solved the PS selection and registration problems in GNSS-based InBSAR systems caused by low signal-to-noise ratio, poor resolution, low DEM accuracy, and large satellite positioning errors.

[0071] 2. The PS point offset was corrected, and an interference phase with high phase accuracy was obtained, laying a good foundation for subsequent high-precision three-dimensional deformation measurement.

[0072] 3. It improves the robustness of the GNSS-based InBSAR system, reduces the error of the interferometric phase, and plays a significant role in the practical application of GNSS-based InBSAR deformation monitoring. Attached Figure Description

[0073] Figure 1 Configure a model for a GNSS-based InBSAR system;

[0074] Figure 2 Here is a flowchart of the algorithm proposed in this invention;

[0075] Figure 3 A schematic diagram of PS pre-selected point association sequences overlaid based on generalized regions;

[0076] Figure 4 A schematic diagram of deduplication processing for PS pre-selected point association sequences;

[0077] Figure 5 The above describes the experimental scenario and required equipment for embodiments of the present invention.

[0078] Figure 6, is the theoretical resolution unit of PS;

[0079] Figure 7 Statistical analysis of coherence coefficients for all resolution units;

[0080] Figure 8 The image shows the PS pre-selected point extraction results in the SAR image of day 1 of the example.

[0081] Figure 9 , represents the number of PS pre-selected points within 32 days in the example;

[0082] Figure 10 To introduce a resolution cell image with generalized area;

[0083] Figure 11 The image shown is an overlay of the generalization region on day 1 of the embodiment.

[0084] Figure 12 The image shows the overlay of the total generalization region over 32 days in Example 3;

[0085] Figure 13 The above is a superimposed image of the generalized area after 32 days of filtering in the example.

[0086] Figure 14 1. PS preselection point lookup table: initial sequence (left), after removing duplicate pixels (right);

[0087] Figure 15 , where is the coherence coefficient of the adjacent PS preselected point pairs in the associated sequence of the 30th PS preselected point;

[0088] Figure 16 , represents the number of associated sequences for the pre-selected PS points after filtering;

[0089] Figure 17 , represents the magnitude deviation of the PS preselected point associated sequence;

[0090] Figure 18 The above represents the Beidou3-igso3 PS selection results in the embodiment.

[0091] Figure 19 The figures show the deformation monitoring results of two points, PS1 and PS2, in the embodiment.

[0092] Figure 20 The accuracy of deformation monitoring for 8 IGSO satellites;

[0093] Figure 21 The time spent processing 8 IGSO satellites. Detailed Implementation

[0094] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0095] GNSS-based InBSAR system configuration as follows Figure 1 As shown, an in-orbit navigation satellite serves as the transmitter, and a receiver fixed on the Earth's surface receives the target's reflected signal. The offset of the satellite's repeating channel (satellite baseline and positioning error vector sum) D and the DEM error Δh cause the focus position to shift. The following uses an East-North-Up coordinate system, with the receiver's position as the origin, to establish the signal model.

[0096] Under the conditions of random orbital error and elevation angle error, a phase model for target Q is established:

[0097]

[0098] In the formula, σ is the scattering coefficient, and P T For the transmitter position, P Q P is the location of the target. E Let Q be the location of the receiver. At non-target locations, the ambiguity function is used to represent the ambiguity between target vector Q and its neighboring target vector A. This ambiguity function can be expressed as:

[0099]

[0100] In the formula Φ TQ and Φ RQ Let θ be the unit vector from the transmitter and receiver to the target Q, respectively; β be the bistatic angle; Θ be the direction vector along the bisector of β; and ω be the angle bisector. E Ξ represents the equivalent angular velocity and equivalent motion direction, p represents the result of range pulse compression, and m represents the distance pulse compression. A This is the result of azimuth pulse compression, where λ is the wavelength and c is the speed of light. The target is displayed on the image in resolution units.

[0101] However, due to DEM errors, satellite positioning errors, and the influence of satellite spatial baselines, resolution cells exhibit random offsets. The projection offset, related to DEM errors and satellite positions, can be expressed as:

[0102]

[0103] In the formula () xy For the horizontal component, () z For the vertical component, V S For satellite velocity, G Q′ It can be represented as:

[0104]

[0105] A method for PS selection and registration in a GNSS-based InBSAR system based on generalized regions, flowchart as follows: Figure 2 As shown, firstly, a generalization region is established to identify the associated sequences of pre-selected PS points. Then, duplicate pixels in the associated sequences of PS pre-selected points are removed, and non-PS points are removed using a coherence threshold. Finally, phase is extracted based on distance weights and phase weights.

[0106] Step 1: Select PS pre-selection points

[0107] Using radar SAR images obtained through the BP imaging algorithm on a GNSS-based InBSAR system, theoretical resolution cells are determined based on radar configuration and experimental parameters. Points with coherence coefficients greater than a threshold are selected as pre-selected PS points.

[0108] Perform PS pre-selection, and define the theoretical resolution cell area of ​​target Q as:

[0109] S PSF (Q)={A|χ(Q,A)>-3dB} (24)

[0110] In the formula, A represents the pixels surrounding Q, and the actual resolution unit area of ​​the target Q can be expressed as:

[0111] S(A) ACT ={A|I(A)>I(Q)-3dB} (25)

[0112] Where I is the amplitude of the SAR image. The correlation coefficient between the PS resolution unit extracted from the actual image and the theoretical resolution unit is:

[0113] coff(Q) = cohe(S(Q)) PSF ,S(Q) ACT (26)

[0114] Where cohe is the coherence coefficient of the two images, and points with sufficiently high coherence coefficients are selected as PS pre-selected points.

[0115] Step 2: Calculate the generalization area

[0116] To facilitate the association of pre-selected PS points, this invention employs a generalized area method. The concept of a generalized region is introduced, which replaces the pre-selected PS points.

[0117] In GNSS-based InBSAR, the SAR image is projected onto the imaging plane, and it is assumed that the phase of the PS is consistent with the phase within its generalization region. Through ambiguity function derivation, the expression for establishing the generalization region based on the theoretical resolution unit is:

[0118] S GEN(Q)={A|(σ GEN (AQ)+Q)∈S PSF (Q)} (27)

[0119] Where, σ GEN The scaling factor is determined based on the range of variation of the DEM across the entire scene, and is usually expressed as the variance of the DEM across the entire scene.

[0120] Overlapping the generalized regions of SAR images across multiple reorbiting cycles. For imaging results of a satellite over K consecutive days, find the pre-selected PS points in the daily imaging results, and denote the set of pre-selected PS points for each day as [Y1, Y2, ..., Y]. K This is used for subsequent PS selection. Taking the first image as an example, for the PS pre-selected point set Y1, the generalization area S of the entire image is... YGEN It can be represented as:

[0121]

[0122] Step 3: Extract the association sequence of PS pre-selected points

[0123] The generalized regions of PS pre-selected points over multiple days are overlaid, and the associated sequences of pre-selected points in the image are extracted. Non-PS sequences with a small number of valid days are then removed.

[0124] If S YGEN A value of 1 indicates that the pixel can be effectively measured. The S values ​​of all images are then superimposed. YGEN Find the pixels that are continuously measured over time. Determine the association sequence of the PS preselected points based on the membership relationship between each pixel and the PS preselected points.

[0125] For the imaging results over K days, the expression after generalizing the region overlay is:

[0126] S GEN_ALL =S YGEN (Y1)+S YGEN (Y2)+…S YGEN (Y K (29)

[0127] For S GEN_ALL Threshold filtering is performed on the value of each pixel:

[0128] F GEN ={(x,y)|G GEN_ALL(x,y) ≥thre GEN} (30)

[0129] Where F GEN For the selected pixels, then GENThis represents the minimum threshold for the number of days a pixel can be detected. The association sequence of PS pre-selected points is determined based on the membership relationship between each pixel and the PS pre-selected points, as illustrated in the diagram below. Figure 3 As shown.

[0130] Introducing the coefficient ζ GEN This compensates for the direction error caused by the large resolution difference between range and azimuth in GNSS-based InBSAR.

[0131]

[0132] In the formula, R1 is the distance between the registered pixel and the PS pre-selected point, and R2 is the length of the generalization region in the R1 direction. When a pixel belongs to two different PS pre-selected points, its corresponding ζ is calculated separately. GEN Choose ζ GEN Smaller PS pre-selection points.

[0133] Step 4: Remove duplicate pixels

[0134] Since each resolution unit contains multiple pixels, there are many duplicates in the results. The PS pre-selected point association sequences are deduplicated, retaining as many longer sequences as possible. Then, a sequential coherence coefficient filtering method is used to retain sequences with higher coherence coefficients.

[0135] The PS pre-selected point association sequences in the above processing have a large number of repetitions. To ensure that the phase of each PS can only be used once, longer sequences are selected first for deduplication processing. The operation process is as follows: Figure 4 As shown.

[0136] 1) Arrange the pre-selected PS points in descending order of sequence length based on the previously obtained PS pre-selected point association sequence.

[0137] 2) Remove elements from the subsequent sequences that are identical to those in the current sequence to perform deduplication filtering. Then, rearrange the filtered results in descending order of sequence length.

[0138] 3) Repeat filtering and sorting until all duplicate pixels are removed.

[0139] This method preserves the integrity of long sequences as much as possible. A sequential coherence filtering method is employed, further filtering each sequence based on resolution unit coherence. Assume a sequence has L pairs of adjacent pre-selected PS points, each pair representing Q1 and Q... l+1 Calculate the coherence coefficient using the resolution unit:

[0140] coff(Q1,Q l+1 )=cohe(S(Q1) PSF ,S(Q l+1 )PSF (32)

[0141] The threshold γ for the inter-resolution coherence coefficient is calculated based on the theoretical coherence coefficient γ as follows:

[0142] γ=γ tem ·γ noi ·γ spa (33)

[0143] In the formula γ tem γ is the time coherence coefficient. noi γ is the noise coherence coefficient. spa Let be the spatial coherence coefficient. The scattering characteristics of PS do not change over a short period of time, meaning the temporal coherence coefficient is approximately 1. The relationship between the noise coherence coefficient and the signal-to-noise ratio can be expressed as:

[0144]

[0145] When the coherence coefficient of a set of pre-selected PS points exceeds T C At the same time, subsequent processing is performed to make the interferometric phase monitoring results of the same target more continuous.

[0146] Step 5: Select PS

[0147] After obtaining the correlation sequence of pre-selected PS points, the PS is determined using the amplitude deviation threshold method. The interferometric phase is extracted from the PS obtained by this algorithm and applied to the next step of high-precision deformation inversion.

[0148] For a PS preselected point association sequence of length N, according to equation (6), the peak amplitude value within each PS preselected point resolution cell is extracted as S(Q)=[Q1,Q2,…Q N Then, calculate its mean m. A =mean(S(Q)) and standard deviation σ Q = std(S(Q)). The expression for amplitude deviation is:

[0149]

[0150] Calculate the magnitude deviation of the associated sequence for each PS preselected point and compare it with the threshold. D If the comparison is greater than the threshold requirement, the associated sequence of the PS preselected point is considered to be PS.

[0151] The phase of the resolution cell containing the PS is extracted to reduce the impact of noise on the accuracy of the interferometric phase. Considering that for each PS, its effective phase extraction range is the intersection of its theoretical and actual resolution cells, according to equations (5) and (6), the phase extraction range S... PHA It can be represented as:

[0152] SPHA (Q1)=S PSF (Q1)∩S ACT (Q1) (36)

[0153] The phase of Q can then be expressed as:

[0154]

[0155] Where M is S PHA The number of pixels, W is the phase weight. p This represents the distance weight. The phase weight and distance weight can be expressed as:

[0156]

[0157] In the formula, P(m) represents the position of the m-th pixel. and σ P For S PHA The phase variance and distance from the peak are calculated. By using a weighted average, the influence of noise can be minimized, thus improving phase accuracy. The resulting PS will be used for multi-angle PS correlation, utilizing interferometric phase to achieve high-precision three-dimensional deformation inversion.

[0158] The specific implementation examples are as follows:

[0159] This embodiment is located in Chongqing, China (30.6855°N, 108.2925°E), and the scene includes buildings, roads, and hillsides, using BeiDou-2 and BeiDou-3 satellites as transmitters. Actual images of the scene, equipment, and antenna in this embodiment are shown below. Figure 5 As shown.

[0160] The radar SAR image was obtained using the BP imaging algorithm, with an imaging range of 1000m (East) * 1400m (North). Figure 8 As shown. This experimental example processed Beidou2-igso2 satellite data, using a single-angle PS (photometric) selection, over a time span of 32 days.

[0161] Select pre-selected PS points from the acquired SAR images. Determine the theoretical resolution units based on radar configuration and experimental parameters, such as... Figure 6 As shown, the distance is 14m and the azimuth is 6m, forming an ellipse. Referring to the coherence coefficient, the pre-selected PS points are extracted using theoretical resolution cells, and the results are as follows. Figure 7 As shown, when the cohesion threshold is set to 70%, 76% of the pre-selected PS points have a cohesion coefficient higher than this threshold.

[0162] The PS selection results for day 1 are as follows: Figure 8As shown, there are 237 pre-selected PS points. These PS points are mainly concentrated on bare ground and roads, with almost no points in vegetated areas, accurately reflecting the characteristics of the actual scene. The number of PS pre-selected points over 32 days is statistically analyzed, as shown below. Figure 9 As shown, the scene features densely packed targets and close proximity between photocells (PSs). However, the distance sidelobes of resolution cells do not attenuate, leading to mutual interference between different PSs. The number of pre-selected PSs within 32 days ranges from 200 to 340, exhibiting significant distribution variation. Therefore, a novel PS selection and registration algorithm is needed to handle this fluctuation.

[0163] This algorithm is used to process the PS pre-selected point sequence. Based on the shape of the theoretical resolution cell and the accuracy of the DEM, the generalization factor is set to 1.2, and the generalization region is obtained as follows: Figure 10 As shown, the yellow area represents the theoretical resolution unit. A generalized region is applied to the daily PS pre-selected point results; the generalized region for day 1 is shown below. Figure 11 As shown. Using the same method, the 32-day generalization region is superimposed, as shown. Figure 12 As shown, the value of each point represents the number of days of effective registration. Setting the minimum registration threshold to 22, the associated sequences of PS pre-selected points less than 22 days are removed from the overlay image. The results are as follows. Figure 13 As shown, it can be observed that the PS pre-selected point association sequences with larger values ​​are concentrated around the ground, which is consistent with the characteristics of the actual scene.

[0164] Sort all associated sequences in descending order of length, such as Figure 14 As shown in the left image. After removing duplicate pixels and re-sorting, the result is as follows. Figure 14 As shown in the right figure, a total of 443 PS pre-selected point association sequences were extracted. Adjacent PS pre-selected sequences were filtered based on the coherence coefficient, with a threshold γ set to 0.95. Figure 15 The coherence coefficient calculation results for the PS preselected point pairs in the 30th PS preselected point sequence show that the coherence coefficients of 5 adjacent point pairs are below 0.95. PS preselected point association sequences that do not meet the coherence coefficient requirements are filtered out, and the results are as follows. Figure 16 As shown, 275 PS pre-selected point association sequences were ultimately retained.

[0165] After obtaining the association sequences of the pre-selected PS points, the magnitude deviation threshold method is used for PS selection. The magnitude deviation results of all PS pre-selected point association sequences are as follows: Figure 16 As shown in the figure. The threshold for amplitude deviation was set to 0.4. A total of 270 PS pre-selected point association sequences met the threshold, and the final distribution is shown in the figure. Figure 17 As shown, regions with high SNR have more PS (Simultaneous Predictive Values). The number of PS selected from a single satellite at a single angle is relatively small; however, the number of PS increases significantly when results from multiple satellites are combined.

[0166] The interference phase of the PS was extracted for deformation monitoring. Since this experiment lasted only 32 days, changes in the experimental setting can be considered invisible. Figure 18 Taking the red regions PS1 and PS2 as examples, after error phase compensation, the deformation results are as follows: Figure 19 As shown. The blue line represents the result without error phase compensation; the orange line represents the final deformation monitoring result. The deformation monitoring results of the two PS fluctuate around 0. Its deformation accuracy can be expressed as:

[0167] acc = RMSE(D s1 -D s2 (39)

[0168] Where D s1 and D s2 The figures show the actual deformation and the deformation monitoring results using GNSS-based InBSAR, respectively. The monitoring area is a natural landscape with virtually no deformation, so the actual deformation is 0. The monitoring accuracies of PS1 and PS2 are 9.7 mm and 10.3 mm, respectively.

[0169] The PS (Static Points) of eight IGSO satellites in the BeiDou system were processed, with an average PS count of 265 and a monitoring accuracy as follows: Figure 20 As shown, the average monitoring accuracy is 12.0 mm. Using a computer with 16GB of onboard RAM and an Intel(R) Core(TM) i9-10900 CPU @ 2.80GHz 2.81GHz processor, processing signals from 8 satellites, 3 frequency points, and each segment 600s in length, the algorithm's runtime per run is as follows: Figure 21 As shown, each satellite consumes an average of 4.9 minutes, which meets the requirement of GNSS-based InBSAR to output deformation data once every 2 hours.

[0170] In summary, this invention proposes a novel PS selection and registration algorithm for GNSS-based InBSAR systems based on generalized regions. This algorithm can quickly and accurately extract and filter suitable PSs, achieving high-precision 3D deformation inversion. It effectively solves the problems of low signal-to-noise ratio, poor resolution, low DEM accuracy, and low satellite positioning accuracy existing in GNSS-based InBSAR.

[0171] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A PS selection and registration method for a GNSS-based InBSAR system based on a generalized region, characterized in that, include: Step 1: Select PS pre-selection points: Using the radar SAR image obtained by the BP imaging algorithm through the GNSS-based InBSAR system, determine the theoretical resolution unit based on the radar configuration and experimental parameters, and select points with a coherence coefficient greater than the threshold as PS pre-selection points. Step 2: Calculate the generalization area: To facilitate the association of PS pre-selected points, this invention adopts the generalization area method; the concept of generalization region is introduced, and the generalization region is used to replace the pre-selected PS points; Step 3: Extract the associated sequence of PS pre-selected points: Overlay the generalized regions of PS pre-selected points over multiple days, extract the associated sequence of pre-selected points in the image, and remove the non-PS sequences with a small number of effective days; Step 4: Remove duplicate pixels: Since each resolution unit contains multiple pixels, there are a lot of duplicates in the result; perform duplicate removal processing on the PS pre-selected point association sequence, and keep as long sequences as possible; then use the sequential coherence coefficient screening method to keep the sequences with higher coherence coefficients. Step 5: Select PS: After obtaining the PS pre-selected point association sequence, determine PS using the amplitude deviation threshold method; extract the interferometric phase based on the PS obtained by this algorithm, and apply it to the next step of high-precision deformation inversion.

2. The method as described in claim 1, characterized in that, In step two, the expression for establishing the generalization region based on the theoretical resolution unit is: S GEN (Q)={A|(σ GEN (AQ)+Q)∈S PSF (Q)} Where Q is the target in the image, A is the pixels surrounding Q, and S... PSF (Q) represents the theoretical resolution cell area of ​​target Q, σ GEN This is the scaling factor.

3. The method as described in claim 1, characterized in that, In step two, the set of pre-selected PS points for each day is denoted as [Y1, Y2, ..., Y]. K For the preselected point set Y1 of PS, the generalization area S of the entire image YGEN Represented as: Among them, S GEN (Q) represents the area of ​​the generalized region corresponding to the target Q in the image.

4. The method as described in claim 1, characterized in that, In step three, for the imaging results after K days, the expression after generalization region overlay is: S GEN_ALL =S YGEN (Y1)+S YGEN (Y2)+…S YGEN (Y K ) Among them, S YGEN The generalized area of ​​the entire image; For S GEN_ALL Threshold filtering is performed on the value of each pixel: F GEN ={(x,y)|G GEN_ALL(x,y) ≥thre GEN } Where F GEN For the selected pixels, G GEN_ALL(x,y) The result of overlaying the generalized region at each point, thre GEN This is the minimum threshold for the number of days a pixel can detect.

5. The method as described in claim 1, characterized in that, In step four, assume a sequence has L pairs of adjacent pre-selected PS points, each pair representing Q1 and Q2 respectively. l+1 Calculate the coherence coefficient using the resolution unit: coff(Q1,Q l+1 )=cohe(S(Q1) PSF ,S(Q l+1 ) PSF ) The threshold γ for the inter-resolution coherence coefficient is calculated based on the theoretical coherence coefficient γ as follows: c = c tem ·c noi ·c spa In the formula γ tem γ is the time coherence coefficient. noi γ is the noise coherence coefficient. spa is the spatial coherence coefficient.

6. The method as described in claim 1, characterized in that, In step four, the relationship between the noise coherence coefficient and the signal-to-noise ratio is expressed as: When the coherence coefficient of a set of pre-selected PS points exceeds T C At the same time, subsequent processing is performed to make the interferometric phase monitoring results of the same target more continuous.

7. The method as described in claim 1, characterized in that, In step five, the expression for the amplitude deviation is: Calculate the magnitude deviation of the associated sequence for each PS preselected point and compare it with the threshold. D If the comparison is greater than the threshold requirement, the associated sequence of the PS preselected point is considered to be PS.

8. The method as described in claim 1, characterized in that, In step five, the phase extraction range S PHA Represented as: S PHA (Q1)=S PSF (Q1)∩S ACT (Q1) The phase of Q is then expressed as: Where M is S PHA The number of pixels, W is the phase weight. p This represents the distance weight.

9. The method as described in claim 1, characterized in that, In step five, the phase weight and distance weight are expressed as follows: In the formula, P(m) is the position of the m-th pixel; and σ P For S PHA Mid-phase variance and distance from the peak.

Images