An adaptive spatial filtering method based on navigation satellite bistatic interferometric SAR system

Abstract

The present invention belongs to the field of bistatic interferometric SAR technology, and specifically relates to an adaptive spatial domain filtering method based on a navigation satellite bistatic interferometric SAR system. The invention first defines a search area for PS points and then performs non-local mean filtering on each PS point based on the search area, thereby achieving adaptive spatial domain filtering. This method solves the problem of eliminating multi-source errors through spatial domain filtering due to the system's bistatic configuration changes, low signal-to-noise ratio, and two-dimensional resolution in GNSS-InBSAR systems, and plays an important role in the practical application of GNSS-InBSAR systems.

Description

Technical Field

[0001] The present invention belongs to the technical field of bistatic interferometric SAR (bistatic SAR) and particularly relates to an adaptive spatial domain filtering method based on a navigation satellite bistatic interferometric SAR (bistatic SAR) system. Background Art

[0002] Navigation Satellite System-Based Bistatic Interferometric Synthetic Aperture Radar (GNSS-InBSAR) is a bistatic SAR system that uses navigation satellites as external radiation sources and deploys receivers on or near the ground to receive echoes from target scenes. With the continuous improvement and expansion of navigation satellite systems, the number of navigation satellites continues to increase. At least 16 navigation satellites can be observed from different angles above any location on the Earth's surface. By selecting different satellites, three-dimensional deformation monitoring can be achieved. Furthermore, compared to traditional low-orbit interferometric SAR, navigation satellite-based interferometric SAR systems offer significant advantages, including short re-orbit times, wide coverage, continuous spatial and temporal monitoring, and low cost.

[0003] The GNSS-InBSAR system uses PS (Permanent Scatterer) technology. During the 3D deformation inversion process, the interferometric phases of PS points at different angles are first extracted, and then the 3D deformation results are obtained through correlation. However, the system suffers from a low signal-to-noise ratio and the presence of multiple sources of error that affect the accuracy of the interferometric phases. Furthermore, system parameters dictate that the resolution of navigation satellite images is lower than that of traditional interferometric SAR, blurring image texture information and resulting in discontinuous interferometric phases at the extracted PS points. Therefore, the spatial filtering problem of navigation satellite-based bistatic SAR systems needs to be addressed. Summary of the Invention

[0004] Aiming at the problems existing in the three-dimensional deformation inversion processing of the GNSS-InBSAR system, the present invention provides an adaptive spatial domain filtering method based on a navigation satellite bistatic interferometric SAR system.

[0005] The technical solutions for implementing the present invention are as follows:

[0006] An adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system, the method comprising the following steps:

[0007] The first step is to use a navigation satellite bistatic interferometric SAR system to obtain several SAR images, and obtain the PS points of each SAR image and the interferometric phase corresponding to each PS point;

[0008] The second step is to determine the actual search area of ​​each PS point based on all the SAR images obtained in the first step, where the actual search area of ​​each PS point is the area that interacts with the interferometric phase of the PS point, and the size of the area is determined by the resolution unit of the PS point;

[0009] In the third step, non-local mean filtering is performed on other PS points located in the actual search area of ​​each PS point determined in the second step to complete the adaptive spatial domain filtering based on the navigation satellite bistatic interferometric SAR system.

[0010] In the first step, the method for obtaining the interference phase corresponding to each PS point is: establishing a model of the interference phase of the PS point, and obtaining the interference phase corresponding to each PS point based on the established model. The specific method is:

[0011] Assume there are M satellites with K repeated orbits. The imaging result of the GNSS-InBSAR system contains Q pixels. The phase model of the qth pixel of the kth repeated orbit of the mth GNSS satellite is expressed as:

[0012]

[0013] in, represents the target scattering phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; represents the geometric configuration phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; represents the atmospheric phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; represents the noise phase of the qth pixel of the kth re-orbit of the mth GNSS satellite.

[0014] For each PS point, due to the stable scattering characteristics and high signal-to-noise ratio of the corresponding pixel, the following assumptions hold:

[0015]

[0016] The interferometric phase of each PS point in the GNSS-InBSAR system is expressed as:

[0017]

[0018] ΔLp (k) (q)=Lp (k+1) (q)-Lp (k) (q)

[0019] Where λ represents the system wavelength, represents the atmospheric phase from the receiving path of the kth re-orbit of the mth GNSS satellite, represents the synthetic aperture center vector of the kth re-orbit of the mth GNSS satellite, LR represents the position vector of the receiver, and Lp (k) (q) represents the position vector of the qth pixel in the kth heavy track, ΔLp (k) (q) represents the three-dimensional deformation of the qth PS point between two re-orbits of the satellite;

[0020] The final interferometric phase term of the GNSS-based InSAR re-orbit can be divided into the following four parts: atmospheric phase, noise phase, terrain phase generated by the spatial baseline, and deformation phase.

[0021] In the second step, the search area is defined as:

[0022] As for deformation, in geology, it is believed that within a certain small range, the deformation variables have the same trend, which is reflected in the interference phase of the PS point. It is manifested as the PS points interacting with each other in a small range, and their interference phases have the same change trend. The search area is determined according to the size of the resolution unit of the PS point.

[0023] In the second step, the method for calculating the resolution unit of each PS point is:

[0024] The theoretical resolution of the PS point can be calculated based on the dual-base structure. Assume that the ambiguity function between the position vector A of the PS point A and the position vector B of the adjacent PS point B is expressed as follows:

[0025]

[0026]

[0027]

[0028] Among them, Φ TA and Φ RA are the unit vectors from the transmitter and receiver to PS point A, β is the bistatic angle, Θ is the direction along the β bisector, ω E and Ξ are the equivalent angular velocity and equivalent motion direction, p is the range pulse compression result, m A is the azimuthal pulse compression result, λ is the wavelength, and c is the speed of light;

[0029] The theoretical 3dB resolution unit area is defined as:

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

[0031] In the second step, the method for determining the actual search area of ​​each PS point based on all SAR images is as follows:

[0032] The actual search area of ​​the PS point is a circle with a value of 1.5*ρ. r is the radius, with PS as the center, where ρ r The range resolution of the GNSS-InSAR system can be calculated as follows:

[0033]

[0034] Ω r =Ξ×Z

[0035] Where Ξ is the equivalent direction of motion, Z is the unit vector along the celestial direction, B is the bandwidth of the transmitter signal, c is the speed of light, β is the bistatic angle, Θ is the direction along the β bisector, and Ω r and Ω a is a unit vector in the direction of range resolution and azimuth resolution.

[0036] The search area for each PS point is finally determined as:

[0037] Ρ(A)={B||A,B| <r},r=1.5*ρ r

[0038] For different satellites at the same target location, the size of the resolution unit is inconsistent, so the area of ​​the search area is also inconsistent.

[0039] In the third step, the method of performing non-local mean filtering on each PS point is:

[0040] The interference phase of each PS point itself and the interference phase of all PS points in the actual search area of ​​the PS point are weighted summed;

[0041] When performing weighted summation, the weight value is determined by A and B, where A is the distance between the PS point and all PS points in the actual search area of ​​the PS point, and B is the coherence between all PS points in the actual search area of ​​the PS point and the theoretical resolution unit where the PS point is located;

[0042] For all PS points of any star image, the entire actual search area set can be obtained; for example, for satellite S1, the number of its total PS points is Q1, and the PS point selection result is [PS1, PS2, ..., PS Q1 ], the phase of PS point is [Ψ(PS1),Ψ(PS2),…,Ψ(PS Q1 )].

[0043] The non-local means filtering of the GNSS-InSAR system can be expressed as a weighted summation:

[0044]

[0045] Among them, P represents the search area of ​​PS0, [PS1,PS2,…,PS i ] represents the PS point in P, Ψ(PS i ) indicates PS i The original phase of represents the phase after PS0 filtering, w irepresents the weight coefficient;

[0046]

[0047]

[0048]

[0049]

[0050] Where Z(i) represents the normalization constant, σ d and σ h represents the smoothing parameter, and γ represents the coherence coefficient. PS points with closer distances have a greater impact on the interference phase, and strong coherence indicates similar deformations. Therefore, we adjust the weight to be inversely proportional to distance and proportional to coherence.

[0051] Based on the process of steps 1 to 3 above, the adaptive spatial domain filtering process is completed.

[0052] Beneficial effects

[0053] (1) The method of the present invention realizes adaptive spatial domain filtering by defining a search area for PS points and performing non-local mean filtering on each PS point based on the search area;

[0054] (2) In the method of the present invention, the search area size of each PS point is determined according to the resolution unit size of the GNSS-InBSAR system. Different satellites have different search areas, which can realize adaptive spatial filtering processing.

[0055] (3) The method of the present invention solves the problem of eliminating multi-source errors through spatial filtering due to the system's dual-base configuration changes, low signal-to-noise ratio and two-dimensional resolution in the GNSS-InBSAR system, which plays an important role in the practical application of the GNSS-InBSAR system.

[0056] (4) The method of the present invention defines the weight of the non-local mean filter based on the fact that the interference phase of the PS point with a closer distance has a greater influence, and strong coherence can indicate that the features are similar. The interference phase of the PS point has a stronger influence, so the weight is inversely proportional to the distance and proportional to the coherence.

[0057] (5) The method of the present invention solves the problem that the PS point space is discontinuous and cannot be directly filtered in the spatial domain in the GNSS-InBSAR system, and performs spatial filtering between PS points through non-local mean filtering.

[0058] (6) The present invention discloses an adaptive spatial filtering algorithm based on a navigation satellite bistatic interferometric synthetic aperture radar. The invention first defines a search area for PS points and performs non-local mean filtering on each PS point based on the search area, thereby realizing adaptive spatial filtering processing. This algorithm solves the problem of eliminating multi-source errors through spatial filtering due to the bistatic configuration changes, low signal-to-noise ratio, and two-dimensional resolution of the GNSS-InBSAR system, and plays an important role in the practical application of the GNSS-InBSAR system. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 The GNSS-InBSAR system configuration of the embodiment of the present invention is as follows;

[0060] Figure 2 Schematic diagram of the search area of ​​the PS point in the embodiment of the present invention;

[0061] Figure 3 FIG1 is a diagram of an experimental device for implementing the present invention;

[0062] Figure 4 This is a comparison chart of the one-dimensional deformation accuracy of four satellites before and after filtering in the embodiment of the present invention. DETAILED DESCRIPTION

[0063] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0064] It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments may be combined with each other; and, based on the embodiments in this disclosure, all other embodiments obtained by persons of ordinary skill in the art without creative work are within the scope of protection of this disclosure.

[0065] It should be noted that various aspects of the embodiments within the scope of the appended claims are described below. It should be apparent that the aspects described herein can be embodied in a wide variety of forms, and any specific structure and / or function described herein is merely illustrative. Based on this disclosure, it should be understood by those skilled in the art that an 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 described herein can be used to implement an apparatus and / or practice a method. In addition, other structures and / or functionalities other than one or more of the aspects described herein can be used to implement this apparatus and / or practice this method.

[0066] like Figure 1As shown in the figure, the GNSS-InBSAR system architecture uses an in-orbit GNSS satellite as a transmitter, obtains the interferometric phase through the satellite's repeated orbit, and uses PS points to obtain high-precision interferometric phase. The PS-InSAR method is based on InSAR phase analysis of SAR stacked data sets. Points with low noise levels are extracted from the SAR image set, called permanent scatterer (PS) points. These points can maintain stable scattering characteristics over multiple satellite cycles. This application performs adaptive spatial filtering based on the actual search area of ​​the PS points.

[0067] The present application embodiment proposes an adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system, the steps of which are as follows:

[0068] Calculate the set of actual search areas of all PS points on any satellite image, where the actual search area of ​​a PS point is the area that interacts with the interferometric phase of the PS point;

[0069] Perform non-local mean filtering on each PS point within the search area.

[0070] In another embodiment of the present application, the process of obtaining the actual search area of ​​each PS point is as follows:

[0071] Take the distance l*ρ r is a constraint condition, the distance l*ρ r The deformation variable of any point in the inner region affects the deformation variable of point A, and the search area of ​​the PS point is determined based on this constraint condition;

[0072] In another embodiment of the present application, the specific process of filtering each PS point based on the non-local mean of the search area is as follows:

[0073] The interference phase of each PS point itself and the interference phase of the PS points in the search area are weighted and summed, where the weight value is determined by the distance between the PS point in the search area and the central PS point and the coherence of the theoretical resolution unit where the PS point is located.

[0074] This embodiment takes M satellites as an example to describe in detail the specific method of using the present invention to perform the adaptive spatial filtering method:

[0075] Step 1: Assume there are M satellites and K repeated orbits. The phase model of each pixel of the GNSS-InBSAR system imaging result can be expressed as:

[0076]

[0077] Where q, m, and k represent the pixel, satellite, and reorbit number, respectively. σ, g, a, and n represent the target scattering phase, geometric configuration phase, atmospheric phase, and noise phase, respectively.

[0078] For PS points, due to the stable scattering characteristics and high signal-to-noise ratio of the corresponding pixels, the following assumptions hold:

[0079]

[0080] The interferometric phase of the PS point in the GNSS-InBSAR system can be expressed as:

[0081]

[0082] Where λ represents the system wavelength, represents the synthetic aperture center vector of the kth re-orbit of the mth GNSS satellite, LR represents the position vector of the receiver, and Lp (k) (q) represents the position vector of the qth pixel in the kth heavy track.

[0083] The final interferometric phase term from GNSS-based InSAR re-trajectories can be divided into four components: atmospheric phase, noise phase, terrain phase due to the spatial baseline, and deformation phase. The noise phase affects the accuracy of the final deformation inversion. Therefore, spatial filtering is necessary to eliminate the effects of the noise phase.

[0084] Step 2: Calculate the set of actual search areas of all PS points on each of the M satellites.

[0085] The following is a detailed analysis and explanation of the calculation process of the actual search area of ​​a certain PS point:

[0086] In geology, deformation is considered to have the same trend within a small range. The interferometric phases of PS points interact with each other within a small range and have the same trend. The search area is affected by the resolution unit of the PS point.

[0087] The theoretical resolution can be calculated based on the dual-base structure. Assume that the ambiguity function between the position vector A of PS point A and the position vector B of its neighboring PS point B is expressed as:

[0088]

[0089]

[0090]

[0091] Among them, Φ TA and Φ RA are the unit vectors from the transmitter and receiver to PS point A, β is the bistatic angle, Θ is the direction along the β bisector, ω Eand Ξ are the equivalent angular velocity and equivalent motion direction, p is the range pulse compression result, m A is the azimuthal pulse compression result, λ is the wavelength, and c is the speed of light.

[0092] The theoretical 3dB resolution unit area is defined as:

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

[0094] The range resolution and azimuth resolution of the GNSS-InSAR system can be calculated by the following formula:

[0095]

[0096]

[0097] Ω r =Ξ×Z

[0098] Ω a =Θ×Z

[0099] Where Z is the unit vector along the sky. B is the bandwidth of the transmitter signal, T s is the synthetic aperture time, λ is the wavelength, Ω r and Ω a is a unit vector in the direction of range resolution and azimuth resolution.

[0100] Take the distance l*ρ r is a constraint condition, the distance l*ρ r The deformation variable of any point in the inner part affects the deformation variable of point A. Based on this constraint, the search area of ​​the PS point is determined, such as Figure 2 shown.

[0101] Ρ(A)={B||A,B| <r},r=l*ρ r

[0102] For different satellites at the same target location, the resolution unit size is inconsistent, so the search area area is also inconsistent. Then, spatial filtering can be performed on each PS point based on the search area P(A).

[0103] Step 3: Perform non-local mean filtering on each PS point in the search area.

[0104] The specific process of this step is:

[0105] For all PS points of any star image, the entire actual search area set can be obtained; for example, for satellite S1, the number of its total PS points is Q1, and the PS point selection result is [PS1, PS2, ..., PSQ1 ], the phase of PS point is [Ψ(PS1),Ψ(PS2),…,Ψ(PS Q1 )].

[0106] The non-local means filtering of the GNSS-InSAR system can be expressed as a weighted summation:

[0107]

[0108] Among them, P represents the search area of ​​PS0, [PS1,PS2,…,PS i ] represents the PS point in P, Ψ(PS i ) indicates PS i The original phase of represents the phase after PS0 filtering, w i represents the weight coefficient;

[0109]

[0110]

[0111]

[0112]

[0113] Where Z(i) represents the normalization constant, σ d and σ h represents the smoothing parameter, and γ represents the coherence coefficient. PS points with closer distances have a greater impact on the interference phase, and strong coherence indicates similar deformations. Therefore, we adjust the weight to be inversely proportional to distance and proportional to coherence.

[0114] Based on the process of steps 1 to 3 above, the adaptive spatial domain filtering process is completed.

[0115] Example

[0116] In this embodiment, the experimental simulation is conducted on a scene located in Changshu City, Jiangsu Province, China (31.7582N, 120.9323E). SAR images are acquired using 8 MEO satellites. Taking the No. 1 satellite (Beidou II MEO 3) as an example, the scene and equipment are as follows: Figure 3 The scene consists of a lake with an artificial road and buildings. Most of the scene is low-resolution vegetation.

[0117] First, calculate the search area. Assuming l = 1.5, based on the dual-base configuration, we can calculate the resolution unit size of a PS point, and then calculate the size of the search area. The average one-dimensional deformation inversion accuracy generated by the four BeiDou-2 satellites before and after filtering is shown in the figure below. Figure 4The points of different shapes in the figure represent the deformation inversion results of the PS point at that location. Figure 4 (a)(c)(e)(g) show the deformation inversion results of each PS point in the scene before filtering. Figure 4 (b)(d)(f)(h) show the deformation inversion results of each PS point in the filtered scene.

[0118] Finally, the mean square error of the one-dimensional deformation of the monitoring area before and after filtering for each satellite is shown in Table 1, with an average reduction of 18.8%. The one-dimensional deformation of multiple satellites was converted into three-dimensional deformation using a multi-angle correlation algorithm. Before filtering, the mean square error of the three-dimensional deformation of the monitoring area was 49.3 mm, 29.5 mm, and 20.3 mm (E, N, U), respectively. After spatial filtering using this algorithm, the mean square error of the three-dimensional deformation of the monitoring area was 34.4 mm, 24.3 mm, and 16.8 mm (E, N, U), respectively. The results show that the overall accuracy has been improved, and the GNSS-InBSAR system has strong three-dimensional interferometric processing capabilities. At the same time, it is demonstrated that the spatial filtering problem of the navigation satellite bistatic interferometric SAR system has been effectively solved.

[0119] Table 1

[0120] Sat.1 Sat.2 Sat.3 Sat.4 Sat.5 Sat.6 Sat.7 Sat.8 One-dimensional deformation mean square error before filtering (mm) 24.4 23.4 22.5 22.9 24.2 20.5 23.6 22.4 Mean square error of one-dimensional deformation after filtering (mm) 19.8 19.3 18.0 20.9 19.6 17.9 18.1 15.7

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

Claims

1. An adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system, characterized in that The steps of the method include: The first step is to use a navigation satellite bistatic interferometric SAR system to obtain several SAR images, and obtain the PS points of each SAR image and the interferometric phase corresponding to each PS point; The second step is to determine the actual search area of ​​each PS point based on all the SAR images obtained in the first step, where the actual search area of ​​each PS point is the area that interacts with the interferometric phase of the PS point, and the size of the area is determined by the resolution unit of the PS point; In the third step, non-local mean filtering is performed on other PS points within the actual search area of ​​each PS point determined in the second step, thereby completing the adaptive spatial filtering based on the navigation satellite bistatic interferometric SAR system. In the first step, the method for obtaining the interference phase corresponding to each PS point is: establishing a model of the interference phase of the PS point, and obtaining the interference phase corresponding to each PS point based on the established model; The method for establishing the interference phase model of the PS point is: Assume there are M satellites with K repeated orbits. The imaging result of the GNSS-InBSAR system contains Q pixels. The phase model of the qth pixel of the kth repeated orbit of the mth GNSS satellite is expressed as: in, represents the target scattering phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; represents the geometric configuration phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; represents the atmospheric phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; represents the noise phase of the qth pixel of the kth re-orbit of the mth GNSS satellite; Obtaining the interferometric phase of each PS point refers to obtaining the interferometric phase of each PS point in the GNSS-InBSAR system, which is expressed as: ΔLp (k) (q)=Lp (k+1) (q)-Lp (k) (q) Where λ represents the system wavelength, represents the atmospheric phase from the receiving path of the kth re-orbit of the mth GNSS satellite, represents the synthetic aperture center vector of the kth re-orbit of the mth GNSS satellite, LR represents the position vector of the receiver, and Lp (k) (q) represents the position vector of the qth pixel in the kth heavy track, ΔLp (k) (q) represents the three-dimensional deformation of the qth PS point between two re-orbits of the satellite; In the second step, the method for calculating the resolution unit of each PS point is: Assume that the fuzzy function expression between the position vector A of PS point A and the position vector B of its neighboring PS point B is: Among them, Φ TA and Φ RA are the unit vectors from the transmitter and receiver to PS point A, β is the bistatic angle, Θ is the direction along the β bisector, ω E and Ξ are the equivalent angular velocity and equivalent motion direction, p is the range pulse compression result, m A is the azimuthal pulse compression result, λ is the wavelength, and c is the speed of light; Then the resolution unit area of ​​PS point is: D(A) P ={B”|χ(A,B)>-3dB} B" is the bandwidth of the transmitter signal.

2. The adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system according to claim 1, characterized in that: The interferometric phase terms of the GNSS-InSAR system include atmospheric phase, noise phase, terrain phase generated by the spatial baseline, and deformation phase.

3. The adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system according to claim 2, characterized in that: In the second step, the method for determining the actual search area of ​​each PS point based on all SAR images is as follows: The actual search area of ​​the PS point is a circle with a value of 1.5*ρ. r is the radius, with PS as the center, where ρ r The range resolution of the GNSS-InSAR system is calculated as follows: Oh r =Ξ×Z Where Ξ is the equivalent direction of motion, Z is the unit vector along the celestial direction, c is the speed of light, β is the bistatic angle, Θ is the direction along the β bisector, and Ω r and Ω a is a unit vector in the direction of range resolution and azimuth resolution.

4. The adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system according to claim 3, characterized in that: In the second step, the search area of ​​each PS point is finally determined to be: P(A)={B||A,B| <r},r=1.5*ρ r 。 5. The adaptive spatial filtering method based on a navigation satellite bistatic interferometric SAR system according to claim 4, characterized in that: In the third step, the method for performing non-local mean filtering on each PS point is: The interference phase of each PS point itself and the interference phase of all PS points in the actual search area of ​​the PS point are weighted summed; When performing weighted summation, the weight value is determined by A' and B', where A' is the distance between the PS point and all PS points in the actual search area of ​​the PS point, and B' is the coherence between all PS points in the actual search area of ​​the PS point and the theoretical resolution unit where the PS point is located; For all PS points of any star image, the entire actual search area set is obtained.

Images