A Multi-Angle Image Correlation Method Based on a Bistatic Interferometric SAR System for Navigation Satellites

Abstract

This invention provides a multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites, belonging to the field of bistatic interferometric SAR technology. The specific process of this method is as follows: First, calculate the set of actual control regions of all PS points on any satellite image, wherein the actual control region of a PS point is the region formed by superimposing the control regions with the same deformation as the PS and the set of 3dB resolution units centered on the boundary of the control regions with the same deformation; Second, superimpose the actual control regions corresponding to all satellites, filter out the associated regions, wherein the associated regions are the actual control regions belonging to three or more satellites, and associate the satellite images corresponding to the associated regions.

Description

Technical Field

[0001] This invention belongs to the field of bistatic interferometric SAR technology, specifically relating to a multi-angle image association method based on a navigation satellite bistatic interferometric SAR system. Background Technology

[0002] Bistatic Interferometric Synthetic Aperture Radar (GNSS-InBSAR) is a bistatic SAR system that uses navigation satellites as external radiation sources and deploys receivers on the ground or near the ground to receive echoes from the target scene. With the continuous improvement and expansion of navigation satellite systems, the number of navigation satellites is constantly increasing. At least 16 navigation satellites at different angles can be observed over any location on the Earth's surface. By selecting different satellites, 3D deformation monitoring can be achieved. Furthermore, compared to traditional low-Earth orbit interferometric SAR, navigation satellite-based interferometric SAR systems have significant advantages such as shorter re-orbit time, wider coverage, continuous spatiotemporal monitoring, and lower cost.

[0003] GNSS-InBSAR systems employ PS (Polarization Point) technology. In 3D deformation inversion processing, it's necessary to obtain the deformation of a PS point in at least three directions, which requires correlating one-dimensional deformation measurements from different angles. However, random errors in the system and terrain information errors can cause shifts in the focus position of the PS point. Furthermore, system parameters dictate that the resolution of navigation satellite images is lower than that of traditional interferometric SAR, resulting in blurred image texture information and hindering multi-angle image correlation. Therefore, the problem of multi-angle image correlation based on navigation satellite bistatic SAR systems needs to be addressed. Summary of the Invention

[0004] To address the problems existing in the three-dimensional deformation inversion processing of the GNSS-InBSAR system, this invention provides a multi-angle image association method based on a bistatic interferometric SAR system of navigation satellites.

[0005] The technical solution for implementing the present invention is as follows:

[0006] A multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites, the specific process of which is as follows:

[0007] Calculate the set of actual control regions of all PS points on any satellite image, where the actual control region of a PS point is the region formed by the superposition of the control regions with the same deformation as the PS and the set of 3dB resolution cells centered on the boundary of the control regions with the same deformation.

[0008] The actual control areas corresponding to all satellites are overlaid, and related areas are filtered out. The related areas are the actual control areas belonging to three or more satellites. The satellite images corresponding to the related areas are then associated.

[0009] Furthermore, the process for obtaining the actual control area of ​​each PS point in this invention is as follows:

[0010] Using distance r as a constraint, the deformation of any point within distance r is the same as that of point A. Based on this constraint, the control area of ​​point PS is determined.

[0011] Determine the boundary of the PS point control region and calculate the set of 3dB resolution cells centered on the boundary;

[0012] The actual control area of ​​the PS point is obtained by superimposing the set of 3dB resolution cells centered on the boundary of the control area.

[0013] Furthermore, the specific process for filtering out the associated regions described in this invention is as follows:

[0014] For each satellite, the pixel value within the actual control area is set to 1, and the pixel value at the location of the area is set to 0. Then, the actual control areas of all satellites are superimposed, and the areas with pixel values ​​greater than or equal to 3 in the superimposed image are filtered out. The filtered areas are the areas where the image is related, and three-dimensional deformation inversion can be achieved based on these related areas.

[0015] Furthermore, the present invention further includes normalizing the associated regions after screening them out, obtaining the normalized overlap of each associated region, and removing associated regions whose normalized overlap is less than a certain threshold.

[0016] Furthermore, the normalization process for the associated region described in this invention is as follows:

[0017] There are E satellites in the associated region F1. Calculate the maximum overlapping area. If the set of actual control regions of E stars is given, then the overlap degree of the associated region F1 is normalized to obtain:

[0018]

[0019] Furthermore, the present invention also includes correcting the position of the PS point using elevation information.

[0020] Furthermore, the location of the PS point in this invention is the average of the locations of multiple images acquired by the same satellite as the final PS point location.

[0021] Furthermore, the process of correcting the position of the PS point using elevation information as described in this invention is as follows:

[0022] Let the step size of the DEM error, dDEM, range from 1 to 5. Calculate the corresponding associated regions for each region and select the largest one from the 5 associated regions as the final associated region.

[0023] Beneficial effects:

[0024] First, this invention defines the control region of PS points and, based on the overlay of time-series image control region maps, initially realizes multi-angle image association. This method solves the problem that in GNSS-InBSAR systems, due to the bistatic configuration changes and low two-dimensional resolution, it is difficult to achieve PS point association at different angles, and plays an important role in the practical application of GNSS-InBSAR systems.

[0025] Second, the present invention calculates and utilizes normalized overlap and removes associated regions with normalized overlap less than a set threshold, thereby compensating for the impact of random errors.

[0026] Third, this invention compensates for the position of the PS point by utilizing elevation measurement error, thereby improving the accuracy of PS point position measurement, achieving maximum correlation between images from different angles, and maximizing the efficiency of three-dimensional deformation monitoring. Attached Figure Description

[0027] Figure 1 This is a GNSS-InBSAR system configuration described in the embodiments of the present invention.

[0028] Figure 2 This is a flowchart illustrating an embodiment of the present invention.

[0029] Figure 3 This is a schematic diagram of the control area of ​​the PS point in the embodiment of the present invention.

[0030] Figure 4 This refers to the positional relationship between the phase center of the direct wave antenna and the target position in the embodiments of the present invention.

[0031] Figure 5 The MEO3 satellite imaging results and PS point locations are shown in the embodiments of this invention.

[0032] Figure 6 This is a schematic diagram of the control area of ​​a GNSS satellite in an embodiment of the present invention.

[0033] Figure 7 This is the association result of multi-angle images in the embodiments of the present invention. Detailed Implementation

[0034] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0035] It should be noted that, in the absence of conflict, the following embodiments and features can be combined with each other; and, based on the embodiments of this disclosure, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this disclosure.

[0036] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this disclosure, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0037] like Figure 1 The illustrated GNSS-InBSAR system architecture uses an in-orbit GNSS satellite as the transmitter, acquires interferometric phase through the satellite's repeating orbits, and obtains high-precision interferometric phase using PS points. The PS-InSAR method is based on InSAR phase analysis of SAR stacked datasets, extracting points with low noise levels from the SAR image set, called permanent scatterer (PS) points. These points maintain stable scattering characteristics across multiple satellite periods. This application uses the actual control area of ​​the PS points to obtain the correlation of satellite images.

[0038] This application proposes a multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites. The flowchart of this invention is shown below. Figure 2 As shown, the steps are as follows:

[0039] Calculate the set of actual control regions of all PS points on any satellite image, where the actual control region of a PS point is the region formed by the superposition of the control regions with the same deformation as the PS and the set of 3dB resolution cells centered on the boundary of the control regions with the same deformation.

[0040] The actual control areas corresponding to all satellites are overlaid, and related areas are filtered out. The related areas are the actual control areas belonging to three or more satellites. The satellite images corresponding to the related areas are then associated.

[0041] In this embodiment, the actual control area of ​​a PS point is defined as the area formed by the superposition of the control area with the same deformation as the PS and the set of 3dB resolution units centered on the boundary of the control area with the same deformation. Based on the overlap of the actual control areas, related satellite images can be obtained, and thus the correlation of satellite images can be quickly realized using this method.

[0042] In another embodiment of this application, the process of obtaining the actual control area of ​​each PS point is as follows:

[0043] Using distance r as a constraint, the deformation of any point within distance r is the same as that of point A. Based on this constraint, the control area of ​​point PS is determined.

[0044] Determine the boundary of the PS point control region and calculate the set of 3dB resolution cells centered on the boundary;

[0045] The actual control area of ​​the PS point is obtained by superimposing the set of 3dB resolution cells centered on the boundary of the control area.

[0046] In another embodiment of this application, for each satellite, the pixel value within the actual control area is set to 1, and the pixel value at the location of the area is set to 0; then the actual control areas of all S satellites are superimposed, and the areas with pixel values ​​greater than or equal to 3 in the superimposed image are filtered out. The filtered areas are the areas where the image is related, and three-dimensional deformation inversion can be achieved based on the related areas.

[0047] This embodiment uses binarization and overlay of images to quickly filter out related regions based on the pixel values ​​of the overlay image.

[0048] This embodiment uses S satellites as an example to provide a detailed description of the specific method for multi-angle image association using the present invention:

[0049] Step 1: Assume there are a total of S satellites, and the number of PS points corresponding to the image of each satellite is [M1, M2, ..., M]. S If ], then the set of points to be registered is:

[0050]

[0051] For each PS point, the actual position of the PS point can be obtained in the following way:

[0052] First, for a single satellite, N time-series images can be obtained, which are [I1, I2, ..., I...]. N For the target point P, due to random errors, its position differs in each image. Therefore, its imaging result position in each image is denoted as...

[0053] Secondly, by registering images from a single angle, the imaging results can be corrected to the same position. For example, using the first image as the main image, the target points in N images can be calibrated to the position of target point P in the first image using the calibration value `regis`, that is:

[0054]

[0055] Finally, the location can be identified using the PS point selection method. Then, based on the offset of each image, we obtain... The position error of the lost machine is compensated by averaging it, and the true position of P is assumed to be:

[0056]

[0057] Step 2: Calculate the set of actual control areas of all PS points on each of the S satellites.

[0058] The actual control area of ​​all PS points on each satellite is calculated by superimposing the actual control areas of each individual PS point. The calculation process for the actual control area of ​​a specific PS point will be analyzed and explained in detail below:

[0059] The actual control area of ​​a PS point consists of two parts: the first part is the area with the same deformation as the PS point, which is called the control area; the second part is the set of 3dB resolution cell areas centered on the boundary of the control area.

[0060] Part One:

[0061] In geology, deformation is considered to be consistent within a certain small range; for example, the deformation at a certain point can be considered as the deformation within a small range.

[0062] Using distance r as a constraint, the deformation of any point within distance r is the same as the deformation of point A, and the control range of point A is D(A). R :

[0063] D(A) R ={B||A,B|<r} (3)

[0064] Part Two:

[0065] The theoretical resolution can be calculated based on the bistatic structure. Let the fuzzy function expression between the position vector A of PS point A and the position vector B of its neighboring PS point B be:

[0066]

[0067]

[0068] Where, Φ TA and Φ RA Let be the unit vectors from the transmitter and receiver to point PS A, respectively; β be the bibase angle; Θ be the direction along the bisector of β; and ω be the angle bisector of β. E Ξ and Ξ are the equivalent angular velocity and equivalent motion direction, p is the range pulse compression result, and m A It is the result of azimuth pulse compression, where λ is the wavelength and c is the speed of light.

[0069] The theoretical area of ​​a 3dB resolution cell is defined as:

[0070] D(A) P ={B|χ(A,B)>-3dB} (2)

[0071] Therefore, the actual control area D(A) of point A can be obtained. fus For D(A) R And D(A) R The set of 3dB resolution cells centered on the boundary, such as Figure 3 As shown, the region within the dashed line is D(A). R The 3dB region centered on the dashed line is D(A). P The final synthesized region is D(A). fus .

[0072] Therefore, the actual control area of ​​point A can be obtained as follows:

[0073]

[0074] Step 3: For all PS points of any star image, the complete set of actual control regions can be obtained. Set the pixel value of the actual control region to 1 and the pixel value of the region location to 0 to achieve image binarization. Then, the actual control regions of all S stars are superimposed, and the regions with pixel values ​​greater than or equal to 3 in the superimposed image are filtered out. The filtered regions are the regions in the image that are related. Based on the related regions, three-dimensional deformation inversion can be achieved.

[0075] The specific process for this step is as follows:

[0076] For any satellite image, all PS points can be used to obtain the complete set of actual control regions. For example, for satellite S1, if the total number of PS points is M1, then the actual control region D(S1) of S1 is:

[0077]

[0078]

[0079] Binarizing D(S1) yields A value of 1 is assigned to areas covered by the controlled area; otherwise, it is assigned a value of 0.

[0080] For all S satellites, the overall overlay image is as follows:

[0081]

[0082] Furthermore, it is possible to use images The median value, which is the superposition of multiple images, is used to determine the correlation results of various control points. Interferometric radar can only detect displacement along the radar line-of-sight (LOS), that is, the projection of ground deformation in various directions onto the LOS. To achieve three-dimensional deformation inversion, at least three angles of input are required, therefore... A value greater than or equal to 3 indicates a location where three-dimensional deformation inversion can be performed.

[0083]

[0084] F represents the regions that can be associated, and f represents the number of those regions.

[0085] Based on the process described in steps one through three above, the image association was completed.

[0086] Considering the influence of elevation information on the position of PS point, a fourth step is introduced based on the above steps one to three. In the fourth step, the position of PS point is corrected using elevation information.

[0087] Because the maximum overlap varies depending on the orientation or number of images, a simple determination of correlation based on the area under the vector F is insufficient. For images corresponding to two satellites with similar configurations, their D... fus The overlap is very high; even if two PSs are far apart, their control areas will have a large overlapping area. However, when the configurations of two satellites differ significantly, even if the actual positions of the PS points are exactly the same, their overlapping area will not be large. Therefore, step four is needed to normalize the overlapping area.

[0088] Step 4: Normalize the associated regions to obtain the normalized overlap of each associated region. When the normalized overlap is less than the set threshold, the associated region is removed.

[0089] In the specific implementation of this step, it is assumed that F1 has E satellites. The maximum overlapping area, when points PS are exactly the same, can be considered as the set of regions controlled by stars E. Therefore, the maximum overlapping area can be expressed as:

[0090]

[0091] Then, by normalizing the overlap of the F1-related regions, we obtain:

[0092]

[0093] The final normalized overlap of the associative regions is:

[0094]

[0095] A normalized overlap threshold is set, and regions with an overlap less than the threshold are removed. The remaining associated regions are considered to be regions with actual association, and satellite image association is achieved based on the retained associated regions.

[0096] In practice, since elevation information can affect the actual location of the PS point, it is necessary to consider the impact of elevation information in step five.

[0097] Step 5: Correct the position of PS points based on elevation information. For the corrected positions, calculate the associated region of each PS point according to the process of steps 1 to 3 above, and select the optimal region from the associated regions calculated after multiple corrections.

[0098] The formula for correcting the position of the PS point using elevation information is as follows:

[0099]

[0100] like Figure 4 As shown, let the origin O be the location of the phase center of the direct wave antenna, and its height above the ground be H. A spatial rectangular coordinate system is established with east, north, and vertically upward as the X, Y, and Z axes, respectively. Here, D and E represent the positions of the direct wave and echo antennas, respectively, and Q is the actual position of a target at a height h above the ground. Using the plane of the direct wave antenna as the imaging plane, the imaged position of target Q can be denoted as Q′. Let the satellite position be S and the velocity vector be V. s .

[0101]

[0102] Among them, P S (t k ) for t k The position of the satellite at that time, P Q (t k ) represents the three-dimensional coordinates of Q, P E This represents the three-dimensional coordinates of the echo antenna. The subscripts x, y, and z represent the horizontal and vertical components of the vector, respectively.

[0103] dDEM represents the step size of the DEM error, with a value ranging from 1 to 5. Therefore, five different position corrections can be calculated. Each position correction can be used to calculate the corresponding associated region. The maximum value among these is then selected as the optimal value.

[0104]

[0105] Specifically, in a GNSS system, the focusing position is also affected by DEM information due to a certain directional shift caused by the system configuration. Based on the GNSS system configuration, an expression for the specific offset of the error can be obtained.

[0106] Let the imaging result of target point Q be Q′. Under different planes, the bistatic distance and Doppler frequency of Q and Q′ are equal, that is:

[0107]

[0108] Among them, P s (t k ) for t k The position of the satellite at that time, P Q (t k ) represents the three-dimensional coordinates of Q, P E This represents the three-dimensional coordinates of the echo antenna.

[0109] It can be simplified to:

[0110]

[0111] Based on this, the relationship between the xy-plane offset and the elevation error can be obtained:

[0112]

[0113] in,

[0114]

[0115] The relationship between horizontal offset and vertical error can be calculated based on different satellite configurations. In actual data processing, DEM measurement errors are consistent within a certain range; therefore, we assume that the DEM error remains approximately constant within an R×R region, where R is much larger than r and the length of the resolution cell.

[0116] Within each region, the variation of each satellite image with DEM error can be obtained.

[0117]

[0118] Where dDEM is the step size of the DEM error.

[0119] In practice, the case with the highest correlation count cannot be assumed to correspond to the true DEM error. If the control regions of different images overlap under a given DEM error, it only indicates that these two points have a probability of overlapping in reality. Here, we aim to retain as many correlated points as possible; the removal of correlated points will be performed in the subsequent deformation phase inversion.

[0120]

[0121] The final result is in the overlapping graph In the equation, the normalized degree of overlap for each overlapping area is: A normalized overlap threshold is set. If the threshold is met, the associated region is retained. By selecting the normalized overlap, the associated regions can be obtained. Based on the control areas corresponding to the associated regions, the final multi-angle image association result is obtained.

[0122] Example:

[0123] In this embodiment, an experimental simulation is conducted on a scene located in Changshu City, Jiangsu Province, China (31.7582N, 120.9323E). SAR images are acquired using eight MEO satellites. Taking satellite 1 (BeiDou II MEO 3) as an example, the scene and image results are as follows: Figure 5 As shown. After random error compensation, the position of the PS point is marked in red.

[0124] First, calculate the control region. Let r = 3. Based on the bistatic configuration, the control region of a PS point can be obtained, as follows: Figure 6 As shown in (a), the control area of ​​the entire scene is as follows: Figure 6 As shown in (b). To demonstrate the variation of RCS in different directions, the control area of ​​Beidou III MEO7 is as follows. Figure 6 As shown in (c), the number and distribution of PS points vary greatly. Therefore, the correlation method based on scene information is no longer applicable. The overlap results of the control areas of each satellite image are as follows: Figure 6 As shown in (d), the highly overlapping regions are mainly distributed along the lakeside road and on both sides of the road, consistent with the scattering characteristics. Then, based on the normalized F, the correlation pairs of PS points are obtained.

[0125] Next, the focus position error was compensated. Based on the measuring equipment, the accuracy of the DEM information was 5m, and the DEM error step size was 1m. Then, the maximum normalized overlap area for each pair of PS points was obtained. In this experiment, the threshold for the normalized overlap area was 30%. The final PS point correlation results are as follows: Figure 7As shown in Table 1, the number of PS point association pairs is a total of 118 points that can be processed in three dimensions. Among the 8 satellites, 2 points can be observed simultaneously by 6 satellites. The results demonstrate that the GNSS-InBSAR system has strong three-dimensional interferometric processing capabilities.

[0126] Table 1

[0127] Number of associated PS points 3 4 5 6 The correlation logarithm obtained by ignoring DEM error 71 17 4 2 The correlation logarithm obtained considering DEM error 88 24 4 2

[0128] The focal position shift caused by DEM error disrupts some correlation pairs. The results ignoring DEM error are shown in Table 1. The total number of correlation pairs is 94, a 20% reduction compared to the result considering DEM error. The results demonstrate that DEM error is non-negligible in GNSS-InBSAR systems. Furthermore, this proves that the multi-angle image correlation problem based on a bistatic interferometric SAR system using navigation satellites is effectively solved.

[0129] 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 multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites, characterized in that, Includes the following steps: The set of actual control regions for all PS points on any satellite image is calculated. The actual control region of a PS point is the region formed by superimposing the control regions with the same deformation as the PS and the set of 3dB resolution cells centered on the boundary of the control regions with the same deformation. The process of obtaining the actual control region of each PS point is as follows: using distance r as a constraint, the deformation of any point within the distance r is the same as that of point A. Based on this constraint, the control region of the PS point is determined. The boundary of the control region of the PS point is determined, and the set of 3dB resolution cells centered on the boundary is calculated. The actual control region of the PS point and the set of 3dB resolution cells centered on the boundary of the control region are superimposed to obtain the actual control region of the PS point. The actual control areas corresponding to all satellites are superimposed, and the associated areas are filtered out. The associated areas are the actual control areas belonging to three or more satellites. The associated regions are normalized to obtain the normalized overlap degree of each associated region, and the associated regions with the normalized overlap degree less than the threshold are removed. Associate the satellite images corresponding to the associated regions; The normalization process for the associated region is as follows: Assume associated region F1 has E satellites in... Calculate the maximum overlapping area. If the set of actual control regions of E stars is given, then the overlap degree of the associated region F1 is normalized to obtain:

2. The multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites according to claim 1, characterized in that, The specific process for filtering out the associated regions is as follows: For each satellite, the pixel value within the actual control area is set to 1, and the pixel value at the location of the area is set to 0. Then, the actual control areas of all satellites are superimposed, and the areas with pixel values ​​greater than or equal to 3 in the superimposed image are filtered out. The filtered areas are the areas where the image is related, and three-dimensional deformation inversion is achieved based on the related areas.

3. The multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites according to any one of claims 1-2, characterized in that, It also includes using elevation information to correct the position of the PS point.

4. A multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites according to any one of claims 1-2, characterized in that, The location of the PS point is determined by averaging the locations of multiple images acquired by the same satellite.

5. The multi-angle image association method based on a bistatic interferometric SAR system for navigation satellites according to claim 3, characterized in that, The process of correcting the position of PS point using elevation information is as follows: Let the step size of the DEM error, dDEM, range from 1 to 5. Calculate the corresponding associated regions for each region and select the largest one from the 5 associated regions as the final associated region.

Images