CN116699610A - Beidou InSAR Atmospheric Error Compensation Method Based on Region Division - Google Patents

Abstract

The invention discloses a Beidou InSAR atmospheric error compensation method based on region division. The atmospheric error phase of the propagation region is considered a constant; and for the estimation area, the atmospheric refractive index of the estimation area is estimated by combining imaging and phase results of different satellites and PS point selection results, each PS point in the scene is traversed, the atmospheric refractive index of each point is estimated respectively, the connectivity of the atmospheric refractive index of the scene is judged by adopting a median absolute deviation method, and the atmospheric refractive index correction is carried out on points which do not meet the connectivity, namely outliers. And finally, adding phases in the propagation area and the estimation area to obtain the atmospheric phase of each target, so as to realize atmospheric phase compensation and further improve the accuracy of deformation inversion and the effective early warning of disasters.

Description

Beidou InSAR atmospheric error compensation method based on regional division

Technical Field

The invention belongs to the technical field of bistatic synthetic aperture radars, and particularly relates to a method for compensating an atmospheric error phase by a Beidou satellite bistatic InSAR system.

Background

The InSAR system of Beidou satellite (BeiDou-InSAR, beiDou based Interferometric Synthetic Aperture Radar System) can be used to monitor scene deformation. According to the system, an on-orbit Beidou satellite is used as a transmitter, a static receiving mechanism is arranged on the ground to form a double-base SAR system, and then deformation monitoring is realized by using a heavy-orbit SAR image. The system inherits the advantages of the Beidou positioning system and the radar system, and can realize three-dimensional deformation measurement of an opposite scene through a single device. Compared with the traditional deformation detection method, the system has the advantages of low cost, low monitoring period and the like.

The implementation of three-dimensional deformation inversion requires interference processing of echoes, however, at different data acquisition moments, due to different scene parameters such as temperature, humidity, atmospheric pressure, wind speed, rainwater and the like, the atmospheric refractive index can be changed, so that signal propagation is affected, and echo delay is generated. For two echoes of the same target at different moments, a phase error, namely an atmospheric error, is introduced to influence deformation inversion accuracy, so that in order to realize high-accuracy deformation inversion, an algorithm is required to be provided for compensating the atmospheric error phase of the Beidou InSAR system.

Disclosure of Invention

In view of the above, the invention provides a Beidou InSAR atmospheric error compensation method based on regional division. The propagation path of the receiver to the PS point target is divided into a propagation region and an estimation region. The atmospheric error phase of the propagation region is considered a constant; and for the estimation area, the atmospheric refractive index of the estimation area is estimated by combining imaging and phase results of different satellites and PS point selection results, each PS point in the scene is traversed, the atmospheric refractive index of each point is estimated respectively, the connectivity of the atmospheric refractive index of the scene is judged by adopting a median absolute deviation method, and the atmospheric refractive index correction is carried out on points which do not meet the connectivity, namely outliers. And finally, adding phases in the propagation area and the estimation area to obtain the atmospheric phase of each target, so as to realize atmospheric phase compensation and further improve the accuracy of deformation inversion and the effective early warning of disasters.

The specific scheme of the invention is as follows:

the Beidou InSAR atmospheric error compensation method based on regional division comprises the following steps:

step one, acquiring a PS point set by combining multi-star data;

step two, dividing a propagation area;

step three, estimating area division;

estimating the atmospheric refractive index in the estimation area;

and fifthly, judging the atmosphere refractive index connectivity.

Further, in the fifth step, the method for determining the atmospheric refractive index connectivity is as follows: a median absolute deviation method. And calculating the median absolute deviation of the atmospheric refractive index in the target point supporting area, and considering the point as an outlier when the difference between the atmospheric refractive index and the median of the point is more than three times of the median absolute deviation. The method for compensating the outliers comprises the following steps: and performing linear fitting by using two non-outliers nearest to the outlier, wherein the atmospheric refractive index obtained by fitting is used as a new atmospheric refractive index of the outlier.

The invention has the following beneficial effects:

the invention solves the problem of low three-dimensional deformation inversion precision caused by the atmospheric phase error of the Beidou satellite double-base InSAR system during interference processing, compensates the atmospheric phase error and improves the deformation inversion precision.

Drawings

FIG. 1 is a schematic diagram of an atmospheric phase according to an embodiment of the present invention.

Fig. 2 is a schematic diagram of estimation block division according to an exemplary embodiment of the present invention.

FIG. 3 is a schematic diagram of two exemplary point locations according to an exemplary embodiment of the present invention.

Fig. 4 shows the accuracy of an exemplary point 1 before uncompensated according to an exemplary embodiment of the invention.

Fig. 5 shows the accuracy of an exemplary point 2 before uncompensated, in accordance with an exemplary embodiment of the invention.

FIG. 6 is an exemplary point 1 of the present invention with accuracy after compensating for atmospheric phase errors.

FIG. 7 is an exemplary point 2 of the present invention with accuracy after compensating for atmospheric phase errors.

Detailed Description

The invention will now be described in detail by way of example with reference to the accompanying drawings.

The Beidou InSAR atmospheric error compensation method based on regional division specifically comprises the following steps:

step one: multi-star data combined acquisition PS point set

The method comprises the following steps:

for K satellites S _a1 L S _aK ]The corresponding PS point selection result is thatThe set of points for the atmospheric phase compensation is:

calculating the bistatic path difference of each point on each day according to the respective system configuration:

step two: propagation region division

The method comprises the following steps:

in the Beidou InSAR system, the atmospheric phase of the target is:

wherein ,n_atm The atmospheric refractive index of the current position at the current time is represented. R represents a bistatic range difference.

The introduction of atmospheric errors in the interferometric phase is due to changes in the atmospheric refractive index, which can be expressed as:

in the actual data processing process, it can be assumed that an is within a small range _atm Is a constant value. Thus, equation (4) can be expressed as an accumulated version:

in the formula (5), the atmospheric refractive index within the upward change Δr of the distance is considered to be constant. Dividing the total transmission path into multiple segments, converting the integral of the path into discrete accumulation, so that only the delta n of each segment needs to be estimated _atm The atmospheric error can be compensated.

However, in practical situations, there will be a region between the receiver and the monitored scene where there is no target, such as monitoring deformation of the shore on one side of the river, where no target above the river can be used to estimate the atmospheric phase, and Δn corresponding to each Δr cannot be obtained _atm . The propagation of the receiver to the target is thus divided into two segments: the propagation region and the estimation region are the paths from the receiver to the left boundary of the estimation region, and the schematic diagram is shown in fig. 1, and the atmospheric error in the propagation region is considered to be a constant.

Step three: estimating a region division

The method comprises the following steps:

a schematic diagram of the estimation block division is shown in fig. 2: for each object in the scene, its estimated region is calculated independently.

Taking any PS point Q therein as an example, as shown in fig. 2: the distance of ΔR/2 is first divided to the left and right in the distance direction with the target point as the center, as shown by the red line in FIG. 2. Then, a circle is drawn by taking the center of the target point position and R_ { atm } as the radius, namely, the absolute distance is limited, and the corresponding area is shown as a light blue circle in fig. 2. The area formed by the circle and the two arcs is the estimated area, as shown by the yellow part in fig. 2. By limiting the system distance to the range and the absolute distance between the targets, an estimated zone range for the target point can be obtained.

In the process of performing atmospheric compensation on the Q point, the influence of Q on an estimation result needs to be removed. Because the imaging results of different satellites are combined in the phase calculation process, the focusing positions of the same target point under different angles have deviation, and all the points need to be removed in the atmospheric refractive index estimation process. The concept of a guard region is presented here, considering a circular region with radius r_ { pro } as a guard region around the center of the Q point, as shown by the green circle in fig. 2. In estimating the atmospheric refractive index, it is necessary to remove the dots in the green region and their corresponding phases.

Step four: atmospheric refractive index estimation within an estimation zone

The method comprises the following steps:

in the Beidou InSAR system, the atmospheric error in the propagation area is considered to be a constant, the atmospheric refractive index in the estimation area is constant, and any point Q is arranged in the estimation area of the target point Q _k The atmospheric phase of (2) can be expressed as:

wherein ,represents a constant atmospheric phase within an estimated region of Q, R ₀ (Q; n) represents the shortest bistatic range difference in the estimated region of Q, which can be represented as R ₀ (Q；n)＝R(Q；n)-ΔR/2。Δn _atm (Q; n) represents the atmospheric refractive index in the estimated region of Q. Utilize the point set->Data of (a) can be implemented to +.> and Δn_atm (Q; n).

wherein ,representing all Q _k Average of the point bistatic range differences, +.>Representing all Q _k And (3) averaging the point atmospheric phases. Since an accurate atmospheric phase cannot be obtained in practical situations, the phase extracted by the PS point is used for replacing the phase:

where the number of the elements in the process is,comprises deformation phase->However, the deformation phase occurs only locally, and the occurrence area is smaller than the estimation area, so that the parameter estimation result in the estimation area is not affected. Even if the deformation area is large, the atmospheric parameter eta can be corrected through subsequent connectivity judgment, and the influence of the deformation phase is eliminated. Thus, the first and second substrates are bonded together,can be used as an estimated input value for atmospheric error.

Step five: atmospheric refractive index connectivity determination

The method comprises the following steps:

in the process of atmospheric refractive index estimation, deformation phase can influence the estimation result, cause outliers and influence the compensation of atmospheric errors, so that the outliers in the deformation phase need to be compensated. And judging the connectivity of the atmospheric refractive index by using a median absolute deviation method, and correcting the point which does not meet the connectivity. The median absolute deviation method has the advantages of strong outlier detection capability, simple calculation, high calculation efficiency and the like, is a simple and exquisite outlier detection method, and comprises the following main steps:

for a target value n ₀ The support area is defined as:

n ₀ the median of (2) is defined as:

med(n ₀ )＝median(n∈S _MAD (n ₀ )) (10)

the median absolute deviation is defined as:

MAD(n ₀ )＝c _MAD ×median _i＝1,L ,(|x _i -median _j＝1,L ,(n ₀ )|) (11)

where i, j=1, 2,.. _MAD = 1.4826. When the difference between the point and the median is greater than three times the absolute deviation of the median, the point is considered to be an outlier. Based on this, the atmospheric refractive index can be corrected as:

wherein ,representing a linear fit of two non-outliers closest to the outlier.

So far, the atmospheric error estimation result of the Q point on the nth day is as follows:

all PS points, i.e. point sets, for all satellites throughout the dayThe estimation areas of the three points are calculated respectively, then the atmospheric refractive index estimation and connectivity judgment are carried out, and the atmospheric error estimation results of all the points can be obtained, and the final differential phase is the phase caused by deformation:

in this embodiment, a deformation scene of 1200 m×1200 m is selected, the transmitter selects Beidou2IGSO3, and the two PS points are [157,218], [304,266], and the positions of the two PS points are shown as red points in fig. 3. The interferometric phase accuracy before two exemplary point offsets is shown in fig. 4 and 5, and the average accuracy of two PS points is: 17.4168mm and 9.4958mm. The results after the atmospheric error compensation are shown in fig. 6 and 7, and the accuracy of the two PS points is 7.5528mm and 3.0722mm.

Table 1 precision comparison

According to the results of table 1, it can be seen that the deformation inversion accuracy after being processed by the Beidou InSAR atmospheric phase error compensation method based on region division provided by the invention is improved compared with that before being processed, and the effectiveness of the invention is proved.

In summary, the above embodiments are only preferred embodiments of the present invention, and are not intended to limit the scope of the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (6)

1. The Beidou InSAR atmospheric phase error compensation method based on regional division is characterized by comprising the following steps of:

step one, acquiring a PS point set by combining multi-star data;

step two, dividing a propagation area;

step three, estimating area division;

estimating the atmospheric refractive index in the estimation area;

and fifthly, judging the atmosphere refractive index connectivity.

2. The method of claim 1, wherein in step one, for K satellites [ S ] _a1 …S _aK ]The corresponding PS point selection result is thatThe set of points for the atmospheric phase compensation is:

calculating the bistatic path difference of each point on each day according to the respective system configuration:

3. the method of claim 1, wherein in the second step, in the beidou InSAR system, the atmospheric phase of the target is:

wherein ,n_atm The atmospheric refractive index of the current position at the current moment is represented; r represents a bistatic range difference;

the introduction of atmospheric errors in the interferometric phase is due to changes in the atmospheric refractive index, which can be expressed as:

in the actual data processing process, it can be assumed that an is within a small range _atm Is a constant value; thus, it can be expressed in the form of an accumulation:

4. the method of claim 1, wherein in the third step, the estimated region is calculated for each object in the scene independently. Taking any PS point Q as an example, taking a target point as a center, firstly dividing the distance of delta R/2 to the left and right along the distance, and then drawing a circle by taking the center of a target point and R_ { atm } as a radius, namely limiting the absolute distance and estimating the range of the point;

in the process of performing atmospheric compensation on the Q point, the influence of Q on an estimation result needs to be removed. Because the imaging results of different satellites are combined in the phase calculation process, the focusing positions of the same target point under different angles have deviation, and all the points need to be removed in the atmospheric refractive index estimation process. The concept of a guard region is presented here, considering a circular region with radius r_ { pro } as a guard region, centered on the Q-point.

5. The method as claimed in claim 1, wherein in the fourth step, in the Beidou InSAR system, the atmospheric error in the propagation region is considered as a constant, the atmospheric refractive index in the estimation region is constant, and the atmospheric refractive index in the estimation region of the target point Q is any point Q _k The atmospheric phase of (2) can be expressed as:

wherein ,represents a constant atmospheric phase within an estimated region of Q, R ₀ (Q; n) represents the shortest bistatic range difference in the estimated region of Q, which can be represented as R ₀ (Q；n)＝R(Q；n)-ΔR/2。Δn _atm (Q; n) represents the atmospheric refractive index in the estimated region of Q; utilize the point set->Data of (a) can be implemented to +.> and Δn_atm An estimate of (Q; n);

wherein ,representing all Q _k Average of the point bistatic range differences, +.>Representing all Q _k Averaging the spot atmospheric phase; since an accurate atmospheric phase cannot be obtained in practical situations, the phase extracted by the PS point is used for replacing the phase:

where the number of the elements in the process is,comprises deformation phase->However, the shapeThe variable phase is only locally generated, and the generation area is smaller than the estimation area, so that the parameter estimation result in the estimation area is not influenced; even if the deformation area is large, the atmospheric parameter eta can be corrected through subsequent connectivity judgment, and the influence of the deformation phase is eliminated; thus (S)>Can be used as an estimated input value for atmospheric error.

6. The method of claim 1, wherein in the fifth step, during the estimation of the atmospheric refractive index, the deformation phase affects the estimation result, causes an outlier, affects the compensation of the atmospheric error, and therefore requires the compensation of the outlier; the connectivity of the atmospheric refractive index is judged by a median absolute deviation method, and the point which does not meet the connectivity is corrected; the median absolute deviation method has the advantages of strong outlier detection capability, simple calculation, high calculation efficiency and the like, is a simple and exquisite outlier detection method, and comprises the following main steps:

for a target value n ₀ The support area is defined as:

n ₀ the median of (2) is defined as:

med(n ₀ )＝median(n∈S _MAD (n ₀ ))

the median absolute deviation is defined as:

MAD(n ₀ )＝c _MAD ×median _{i＝1，…，} (|x _i -median _{j＝1，…，} (n ₀ )|)

where i, j=1, 2,.. _MAD = 1.4826. When the difference between the point and the median is more than three times of the absolute deviation of the median, the point is considered as an outlier; based on this, the atmospheric refractive index can be corrected as:

wherein ,representing a linear fit of two non-outliers closest to the outlier;

so far, the atmospheric error estimation result of the Q point on the nth day is as follows:

all PS points, i.e. point sets, for all satellites throughout the dayAnd (3) respectively calculating an estimation area of each point, and then carrying out atmospheric refractive index estimation and connectivity judgment to realize atmospheric error compensation of all points.