CN117075056A - Ground-based synthetic aperture slope radar atmosphere correction method - Google Patents
Ground-based synthetic aperture slope radar atmosphere correction method Download PDFInfo
- Publication number
- CN117075056A CN117075056A CN202310819068.1A CN202310819068A CN117075056A CN 117075056 A CN117075056 A CN 117075056A CN 202310819068 A CN202310819068 A CN 202310819068A CN 117075056 A CN117075056 A CN 117075056A
- Authority
- CN
- China
- Prior art keywords
- phase
- atmosphere
- point
- radar
- atmospheric
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 238000012937 correction Methods 0.000 title claims description 21
- 238000012544 monitoring process Methods 0.000 claims abstract description 28
- 238000012216 screening Methods 0.000 claims abstract description 11
- 238000010587 phase diagram Methods 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 15
- 238000012545 processing Methods 0.000 claims description 9
- 239000006185 dispersion Substances 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 2
- 238000005259 measurement Methods 0.000 description 10
- 238000005516 engineering process Methods 0.000 description 9
- 239000011159 matrix material Substances 0.000 description 6
- 230000002159 abnormal effect Effects 0.000 description 5
- 230000000694 effects Effects 0.000 description 3
- 230000010365 information processing Effects 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 238000005305 interferometry Methods 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009828 non-uniform distribution Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 101150050759 outI gene Proteins 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/40—Means for monitoring or calibrating
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
The invention relates to the technical field of radar surveying, in particular to a method for correcting the atmosphere of a ground-based synthetic aperture slope radar. The method comprises the steps of screening stable atmosphere estimated PS points from pixel points of a two-dimensional complex scattering radar image, and further obtaining interference phases of the stable atmosphere estimated PS points; and setting a phase residual decision threshold based on a binary atmosphere model and the interference phase, rejecting the PS points which do not meet the phase residual decision threshold, resetting the phase residual decision threshold based on the PS points which are not rejected, repeating the rejecting if the relative variation of the phase residual decision threshold meets an iteration condition, otherwise, ending the iteration, completing optimal atmosphere compensation based on the PS points which are not rejected after the iteration is ended, and finally completing deformation inversion. The PS point finally selected by the method is higher in quality and more reliable, the atmospheric estimation accuracy is improved, and the deformation inversion accuracy of the monitoring area is further improved.
Description
Technical Field
The invention relates to the technical field of radar surveying, in particular to a method for correcting the atmosphere of a ground-based synthetic aperture slope radar.
Background
Landslide and debris flow disasters are particularly serious, and in addition, mines, dams and the like bring economic benefits to us, and meanwhile, various potential safety hazards such as landslide, dam deformation, settlement displacement and the like are accompanied. Before the safety accident occurs, the occurrence of the disaster can be predicted in advance by a deformation monitoring means, and the technical and economic losses are effectively reduced. The deformation monitoring technology is used for realizing early warning and alarming of hidden danger disasters in a monitoring field by slowly monitoring the surface displacement of a deformation monitoring field and combining with priori data.
Deformation monitoring radar is a microwave measurement method developed in the last decade. It has unique measurement advantages compared with the traditional means. The non-contact remote monitoring method can effectively solve the problem that contact measurement equipment such as a differential GPS, a level gauge and the like cannot be applied to dangerous situations. Compared with the optical measurement method, the method has the advantages of all weather and weather, can be used under the condition of extremely low visibility, and is not easily influenced by factors such as rain, fog, illumination and the like.
Deformation monitoring radars can be classified into large-scale long-range measurements and local short-range measurements. The remote measurement is mainly based on satellite-borne radar interferometry technology, GPS technology and the like. The defects are that the distance is long, the precision is low, the method is not suitable for high-precision on-site monitoring, and the measurement period is long. It is only suitable for wide area census and preliminary investigation, and cannot be adapted to emergency measurements. The local short-distance measurement is mainly Ground-based synthetic aperture radar (GB-based Synthetic Aperture Radar, SAR), total station and the like, and is suitable for small-scale monitoring, short in measurement period and high in precision.
GB-SAR is widely applied to deformation monitoring fields of mining areas, dams, glaciers, buildings and the like. The working wave bands of the microwave antenna mainly comprise microwave frequency bands of X, ka, ku and the like. With the continuous development of GB-SAR technology, the accuracy of deformation monitoring is also improved continuously, and the disturbance of the atmospheric phase is not allowed to be ignored. Atmospheric changes have become a major disturbance factor in GB-SAR radar deformation monitoring. Because the topography of the monitoring environment is complex, the water vapor, the humidity and the temperature show non-uniform distribution characteristics in space, the atmospheric phase error correction method based on uniform atmosphere assumption can cause deviation of atmospheric correction, and the deformation monitoring precision of the GB-SAR is reduced. Currently, GB-SAR atmospheric correction methods are mainly divided into two categories: one is a method based on meteorological data, which corrects the atmospheric phase according to the relation between the refractive index of electromagnetic wave propagation and the temperature, air pressure and humidity data in a scene; another class of function model-based methods generally assumes that the atmospheric effects have a strong spatial correlation, i.e., that the atmosphere is spatially uniform, where the atmospheric phase compensation is accomplished by parameterized modeling with stable PS-point interference phases.
In summary, the GB-SAR needs to compensate for measurement errors caused by atmospheric disturbance when performing deformation monitoring, and compared with the meteorological data method, the PS-point-based atmospheric compensation algorithm has advantages in that no cooperative target or additional sensor data are needed, the method has stronger independence and is easier to implement, and is widely applied to the GB-SAR deformation monitoring field. Although the parameterized model atmosphere phase compensation method based on the PS point has many advantages, in describing the space atmosphere phase of GB-SAR with a wide range of azimuth observation angles, the method needs to accurately describe the atmosphere phase with complex space variation by using the PS point as accurately as possible, and the key is a reliable PS point selection strategy and method.
Disclosure of Invention
In view of the defects and shortcomings of the existing PS point-based atmospheric compensation technology, the patent provides a ground-based synthetic aperture slope radar atmospheric correction method, which is based on a modeling parameter atmospheric model and a PS point technology, stable atmospheric estimated PS points are screened and extracted through a multi-threshold, unstable atmospheric estimated PS points are removed by an iterative algorithm, finally, the quality of the selected PS points is higher and more reliable, the atmospheric estimation accuracy is greatly improved, and further, the deformation inversion precision of a monitoring area is improved.
In order to achieve the above object, the present invention provides the following technical solutions:
the ground synthetic aperture slope radar atmosphere correction method is characterized by comprising the following steps of:
s1, acquiring a two-dimensional complex scattering radar image, performing interference processing on the radar image to obtain an interference phase diagram of the radar image, performing two-dimensional space unwrapping on the interference phase diagram to obtain an interference phase of a pixel point of the radar image, screening out a stable atmosphere estimated PS point from the pixel point of the radar image, and acquiring the interference phase of the stable atmosphere estimated PS point;
s2, setting a phase residual error judgment threshold based on an interference phase of a binary atmosphere model and the stable atmosphere estimation PS point, rejecting the stable atmosphere estimation PS point which does not meet the phase residual error judgment threshold, resetting the phase residual error judgment threshold based on the stable atmosphere estimation PS point which is not rejected, repeating the rejecting if the relative variation of the phase residual error judgment threshold meets an iteration condition, otherwise, ending the iteration, and finally taking the stable atmosphere estimation PS point which is not rejected as the PS point for optimal atmosphere compensation;
and S3, performing optimal atmosphere compensation based on the optimal atmosphere parameters of the PS points for optimal atmosphere compensation and the interference phase diagram to finish deformation inversion.
As a preferable scheme of the invention, the method for correcting the atmosphere of the ground-based synthetic aperture slope radar specifically comprises the following steps of:
acquiring a two-dimensional complex scattering radar image, and marking the two-dimensional complex scattering radar image as Img k And generates a radar complex scattering image sequence { Img } k From the sequence { Img } k Extracting two adjacent frame images as main and auxiliary images of interference processing, taking complex conjugate of the auxiliary images and multiplying the auxiliary images with the main image, and obtaining an interference phase map Img after taking the main value phase inpha For the interference phase diagram Img inpha Two-dimensional space unwrapping is carried out to obtain the interference phase diagram Img inpha The interference phase phi of the pixel points of (2) can be expressed as:
φ＝φ def +φ atm +φ noise ；
wherein phi is def Is the deformation phase phi of the pixel point atm Is the atmospheric error phase of the pixel point, phi noise Noise phase for pixel points;
as a preferable scheme of the invention, the method for correcting the atmospheric pressure of the ground-based synthetic aperture slope radar specifically comprises the following steps of:
first, the sequence { Img } is calculated k Average correlation coefficient between adjacent frame imagesAnd then according to the radar radiation image Img k Is defined and the image Img k Amplitude A of the middle pixel point (i, j) k (i, j) and phase->Calculating the sequence { Img } k Average relative amplitude +.>Then, the sequence { Img is calculated k Standard deviation delta of amplitude time series of pixel point (i, j) in } A Mean value m from the amplitude time series A Is the amplitude dispersion D of the time sequence complex scattering image A The method comprises the steps of carrying out a first treatment on the surface of the Finally, respectively setting threshold value threshold T 1 、T 2 And T 3 Screening out pixels of the radar image while satisfying +.> And D is A ≤T 3 As the stable atmosphere estimated PS point, and simultaneously obtaining the interference phase of the stable atmosphere estimated PS point.
As a preferable scheme of the invention, the method for correcting the atmosphere of the ground synthetic aperture slope radar comprises the step of threshold value T 1 、T 2 And T 3 The settings of (2) are as follows: the threshold value threshold T 1 To set 0.9, T 2 Set to 5 and T 3 Setting the PS point to be 0.1, wherein the PS point finally screened out at the moment is used as a stable atmosphere estimated PS point; the threshold value threshold T 1 Set to 0.8, T 2 Set to 0 and T 3 And setting the PS point to be 0.2, wherein the PS point screened at the moment is used as a monitoring point for deformation inversion.
As a preferable scheme of the invention, the interference of the stable atmosphere estimated PS point is eliminated by a spatial filtering algorithm in the ground synthetic aperture slope radar atmosphere correction methodNoise phase phi of phases phi noise Thereafter, the deformation phase Φ of the PS point is estimated due to the stable atmosphere def Approximately 0 such that the interference phase phi of the stable atmospheric estimate PS point is equal to its atmospheric error phase phi atm Equal.
As a preferable scheme of the invention, the method for correcting the ground-based synthetic aperture slope radar atmosphere specifically comprises the following steps of:
s21, constructing a binary atmosphere model based on the distance between the radar and the target and the included angle between the target and the radar aperture center, and calculating the atmosphere error phase phi of the stable atmosphere estimated PS point based on the model atm Thereby obtaining an interference phase phi of the stable atmosphere estimated PS point;
s22, establishing a simultaneous equation system of the interference phase, the distance and the azimuth angle of the stable atmospheric estimated PS point, and solving an atmospheric model parameter vector beta by adopting a least square method to obtain an optimized estimated value of the vector betaAnd calculates a vector X composed of the distance and azimuth angle of the stable atmospheric estimated PS point and the optimal estimated value +.>A vector product which is the atmospheric phase error estimated value phi atm ；
S23, based on the atmospheric phase error estimation value phi atm Performing compensation processing on the interference phase of the stable atmosphere estimated PS point to obtain an atmosphere compensation residual phase phi su b, and calculating the atmosphere compensation residual phase phi sub Is the mean of the samples of (2)And variance->Based on the sample mean->Variance->Phase residual decision threshold coefficient setting phase residual decision threshold +.>The calculation formula is as follows:
wherein T is 0 Is a threshold coefficient;
s24, according to the thresholdEliminating the stable atmospheric estimated PS point, if the absolute value of the interference phase of the stable atmospheric estimated PS point is +.>Removing;
s25, repeatedly executing the stable atmospheric estimation PS points which are not removed in the step S24 in the steps S21-S24 to obtain a phase residual error judgment threshold corresponding to the current iteration number mm is more than or equal to 2, and calculating the decision threshold +.>Relative amount of change of (2)If->Iteration threshold N is less than or equal to 0 When 0 therein<N 0 And (3) continuing iteration, otherwise, ending the iteration, and taking the stable atmospheric estimated PS points which are not rejected as PS points for optimal atmospheric compensation when the iteration is ended.
As a preferable scheme of the invention, the method for correcting the ground-based synthetic aperture slope radar atmosphere comprises the following steps of:
wherein phi is atm Represents the atmospheric error phase, lambda is the wavelength of a radar transmitting signal, r is the distance between the radar and a target, theta is the line-of-sight included angle between the target and the center of the radar aperture, and beta 0 Is a constant term, beta 1 Representing coefficients related to the first power of distance, beta 2 Representing coefficients related to the second power of distance, beta 3 Representing coefficients related to the azimuthal first power, beta 4 Representing coefficients related to the azimuthal second power, beta 5 Representing coefficients that are jointly related to the distance first power and the azimuth first power.
As a preferable scheme of the invention, the method for correcting the atmospheric pressure of the ground-based synthetic aperture slope radar comprises the following steps:
Φ＝Xβ+ε
wherein,β＝[β 0 ,β 1 ,…,β 5 ] T ，ε＝[ε 0 ,ε 1 ,…,ε n ] T phi is an n X1-dimensional vector composed of interference phases of n PS points, X is an n X6-dimensional vector composed of distances and azimuth angles of constant 1 and n PS points, beta is a 6X 1-dimensional vector composed of atmospheric model parameters beta, and epsilon is an n X1-dimensional vector composed of phase errors.
As a preferable scheme of the invention, the method for correcting the atmosphere of the ground synthetic aperture slope radar comprises the step of calculating the threshold coefficient T 0 Set to 2.
As a preferable scheme of the invention, the method for correcting the atmosphere of the ground-based synthetic aperture slope radar comprises the step of iterating the threshold value N 0 Set to 0.95.
Compared with the prior art, the invention has the beneficial effects that:
according to the invention, stable atmosphere estimated PS points are screened through multiple threshold values, and then the unstable atmosphere estimated PS points are further removed through an iterative algorithm, so that accurate selection of higher quality PS points is realized, on the basis, GB-SAR atmospheric correction is realized through an atmospheric correction method based on the PS points and a parameterized model, the correction precision is greatly improved, and the urgent requirement of the existing GB-SAR atmospheric correction technology for accurately correcting complex atmosphere in deformation monitoring is met.
Description of the drawings:
FIG. 1 is a schematic flow chart of a method for correcting the atmospheric pressure of a ground-based synthetic aperture slope radar according to embodiment 1;
FIG. 2 is a flow chart of the information processing of the method for correcting the atmospheric pressure of the ground-based SAR slope radar in the embodiment 2;
FIG. 3 is an atmospheric estimated phase curve of a certain distance band before the elimination of the abnormal PS points in the embodiment 1 and the embodiment 2;
FIG. 4 is an atmospheric estimated phase curve of a certain distance band after the elimination of the abnormal PS points in the embodiment 1 and the embodiment 2;
FIG. 5 shows the different atmosphere phase compensation deformation curves of the monitoring points of the embodiment 1 and the embodiment 2.
Detailed Description
The present invention will be described in further detail with reference to test examples and specific embodiments. It should not be construed that the scope of the above subject matter of the present invention is limited to the following embodiments, and all techniques realized based on the present invention are within the scope of the present invention.
Example 1
As shown in FIG. 1, a flow chart of a method for correcting the atmospheric pressure of a ground-based synthetic aperture slope radar comprises the following specific steps:
s1, acquiring a two-dimensional complex scattering radar image, performing interference processing on the radar image to obtain an interference phase diagram of the radar image, performing two-dimensional space unwrapping on the interference phase diagram to obtain an interference phase of a pixel point of the radar image, screening out a stable atmosphere estimated PS point from the pixel point of the radar image, and acquiring the interference phase of the stable atmosphere estimated PS point;
specifically, a two-dimensional complex scattering radar image is acquired and recorded as Img k And generates a radar complex scattering image sequence { Img } k From the sequence { Img } k In }Extracting two adjacent frame images as main and auxiliary images of interference processing, taking complex conjugate of the auxiliary images and multiplying the auxiliary images with the main image, and obtaining an interference phase map Img after taking the main value phase inpha The calculation formula is as follows:
wherein, the interference phase diagram Img is obtained by a complex conjugate operator and a phase principal value inpha Two-dimensional space unwrapping is carried out to obtain the interference phase diagram Img inpha The interference phase phi of the pixel points of (2) can be expressed as:
φ＝φ def +φ atm +φ noise ；
wherein phi is def Is the deformation phase phi of the pixel point atm Is the atmospheric error phase of the pixel point, phi noise Noise phase for pixel points;
further, acquiring the interference phase of the pixel point of the radar image includes the following steps:
first, the sequence { Img } is calculated k Average correlation coefficient between adjacent frame imagesBased on the image sequence { Img } k Calculating each radar image Img from adjacent frame images in the sequence } k Correlation coefficient gamma of middle pixel point k The calculation formula is as follows:
wherein, the size of w1 and w2 is half of the sliding window length of two coordinate dimensions of the radar image pixel point matrix respectively, if the image sequence { Img k The length of the sequence is L, and a correlation coefficient sequence matrix { gamma } with the length of L-1 can be obtained k -wherein k = 1, …, L-1, and calculating the { γ k Average value of the average value to obtain average correlation coefficientThe calculation formula is as follows:
and then according to the radar radiation image Img k Is defined and the image Img k Amplitude A of the middle pixel point (i, j) k (i, j) and phaseCalculating the sequence { Img } k Average relative amplitude +.>Set radar radiation image Img k If the dimension of (i) is M×N, the complex scattering image value Img corresponding to the pixel (i, j) can be calculated k (i, j) and the calculation formula is:
wherein A is k (i, j) andrespectively representing the amplitude and the phase of the pixel point (i, j) in the kth frame image, and obtaining the image sequence { Img by calculation k Average relative amplitude of pixel points (i, j) in ∈>The calculation formula is as follows:
wherein,for the average value of the amplitude of the pixel points (i, j) in the L frame sequence images, the calculation formula is as follows:
for the average value of the amplitude mean value of all the pixel points, the calculation formula is as follows:
then, the sequence { Img is calculated k Standard deviation delta of amplitude time series of pixel point (i, j) in } A Mean value m from the amplitude time series A Is the amplitude dispersion D of the time sequence complex scattering image A Calculating amplitude dispersion D of time sequence complex scattering image A The formula is as follows:
wherein delta A Standard deviation, m, representing pixel amplitude time sequence A The average value of the pixel amplitude time sequence;
finally, respectively setting threshold value threshold T 1 、T 2 And T 3 Screening out pixel points of the radar image and meeting the requirements simultaneously And D is A ≤T 3 As the stable atmosphere estimated PS point, and simultaneously obtaining the interference phase of the stable atmosphere estimated PS point.
Further, in the present embodiment, the threshold value T 1 、T 2 And T 3 The settings of (2) are as follows: the threshold value threshold T 1 To set 0.9, T 2 Set to 5 and T 3 Setting the PS point to be 0.1, wherein the PS point finally screened out at the moment is used as a stable atmosphere estimated PS point; the threshold value threshold T 1 Set to 0.8, T 2 Set to 0 and T 3 And setting the PS point to be 0.2, wherein the PS point screened at the moment is used as a monitoring point for deformation inversion.
Further, the noise phase phi in the interference phase phi of the stable atmosphere estimated PS point is eliminated by a spatial filtering algorithm noise Thereafter, the deformation phase φ of the PS point is estimated due to the stable atmosphere def Approximately 0 such that the interference phase phi of the stable atmospheric estimate PS point is equal to its atmospheric error phase phi atm Equal.
S2, setting a phase residual error judgment threshold based on an interference phase of a binary atmosphere model and the stable atmosphere estimation PS point, rejecting the stable atmosphere estimation PS point which does not meet the phase residual error judgment threshold, resetting the phase residual error judgment threshold based on the stable atmosphere estimation PS point which is not rejected, repeating the rejecting if the relative variation of the phase residual error judgment threshold meets an iteration condition, otherwise, ending the iteration, and finally taking the stable atmosphere estimation PS point which is not rejected as the PS point for optimal atmosphere compensation;
specifically, S2 includes the following steps:
s21, constructing a binary atmosphere model based on the distance between the radar and the target and the sight angle between the target and the radar aperture center, calculating and acquiring the atmosphere error phase based on the model, and simultaneously acquiring the interference phase, wherein the model construction formula is as follows:
wherein phi is atm Represents the atmospheric error phase, lambda is the wavelength of a radar transmitting signal, r is the distance between the radar and a target, theta is the line-of-sight included angle between the target and the center of the radar aperture, and beta 0 Is a constant term, beta 1 Representing coefficients related to the first power of distance, beta 2 Representing coefficients related to the second power of distance, beta 3 Representing coefficients related to the azimuthal first power, beta 4 Representing coefficients related to the azimuthal second power, beta 5 Representing the first power of distance and the first azimuthCoefficients of power joint correlation; because of the interference phase phi of the stable atmosphere estimated PS point and the atmosphere phase phi thereof at this time atm Equal, so the interference phase phi of the stable atmosphere estimated PS point is obtained at the same time;
s22, establishing a simultaneous equation set of the interference phase of the stable atmosphere estimated PS point, the distance and the azimuth angle of the interference phase, wherein the equation set is as follows:
Φ＝Xβ+ε
wherein,β＝[β 0 ,β 1 ,…,β 5 ] T ，ε＝[ε 0 ,ε 1 ,…,ε n ] T phi is an n multiplied by 1 dimensional vector formed by n PS point interference phases, X is an n multiplied by 6 dimensional vector formed by a constant 1 and n PS distances and azimuth angles, and a 6 multiplied by 1 dimensional vector beta, namely an atmosphere model parameter, epsilon is an n multiplied by 1 dimensional vector formed by phase errors;
based on the equation set, solving an atmospheric model parameter vector beta by adopting a least square method to obtain an optimized estimation value of the vector betaThe calculation formula is as follows:
and calculates a vector X consisting of the distance and azimuth angle of the stable atmospheric estimated PS point and the optimized estimated valueA vector product which is the atmospheric phase error estimated value phi atm The calculation formula is as follows:
s23, based on the atmospheric phase error estimation value phi atm Estimating an interference phase of PS points for the stable atmosphereThe bit is compensated to obtain the atmosphere compensation residual error phase phi sub The calculation formula is as follows:
Φ sub ＝Φ-Φ atm ；
and calculates the atmospheric compensation residual phase phi sub Is the mean of the samples of (2)And variance->Based on the sample mean->Variance->Phase residual decision threshold coefficient setting phase residual decision threshold +.>The calculation formula is as follows:
wherein T is 0 Is a threshold coefficient;
further, T 0 The larger the culled point is, the fewer the point is, according to the normal too distribution theory, when T 0 At 3, 0.27% of PS spots are eliminated, and T is taken as 0 At 2, 4.55% of the PS spots are eliminated, and when T 0 When the number is 1, 31.73% of PS points are removed, the number of abnormal PS points are removed has influence on the iteration speed and the accuracy of atmospheric phase error estimation, but too many PS points are removed, part of stable PS points are also removed together, and the accuracy of atmospheric phase error estimation is also influenced, so according to engineering experience, T 0 The method is set to be about 2 for removing, and the removing effect is ideal.
S24, according to the thresholdEliminating the stable atmospheric estimated PS point, if the absolute value of the interference phase of the stable atmospheric estimated PS point is +.>Removing;
s25, repeatedly executing the stable atmospheric estimation PS points which are not removed in the step S24 in the steps S21-S24 to obtain a phase residual error judgment threshold corresponding to the current iteration number mm is more than or equal to 2, and calculating the decision threshold +.>Relative amount of change of (2)If->Iteration threshold N is less than or equal to 0 When 0 therein<N 0 And (3) continuing iteration, otherwise, ending the iteration, and taking the stable atmospheric estimated PS points which are not rejected as PS points for optimal atmospheric compensation when the iteration is ended.
Further, the iteration threshold N 0 The value of (2) is set to 0.95 because of N 0 The more the PS points are stabilized, but the number of remaining stable PS points is reduced, and the more the sample points involved in the phase error estimation are, the smaller the estimated error is, so according to engineering experience and actual monitoring scenario, embodiment N 0 A value of 0.95 can achieve a good estimation effect.
And S3, performing optimal atmosphere compensation based on the optimal atmosphere parameters of the PS points for optimal atmosphere compensation and the interference phase diagram to finish deformation inversion.
Comparing the atmospheric estimated phase curves of the strips in fig. 3 and fig. 4 at a certain distance, after the abnormal PS points are removed, the estimated atmospheric phase errors are more practical, i.e. the phase compensation errors of all PS participating in the atmospheric phase error estimation are the smallest in a statistical sense. In addition, by using the method for correcting the atmospheric pressure of the ground-based synthetic aperture slope radar, as shown in fig. 5, deformation errors of monitoring points (actually non-deformation points, namely, accumulated deformation amount should be 0 mm) caused by atmospheric pressure phase errors are well corrected.
Example 2
As shown in FIG. 2, the information processing flow chart of the method for correcting the atmospheric pressure of the ground-based synthetic aperture slope radar comprises the following specific steps:
s1, acquiring an original echo of a detection area by a ground-based synthetic aperture radar through a mechanical or electronic scanning working mode, and acquiring a two-dimensional complex scattering radar image by using a radar imaging algorithm, wherein the two-dimensional complex scattering radar image is recorded as Img k Wherein k is a radar scan number, also known as a radar image number;
s2, processing the main Lei Dafu scattering image and the auxiliary Lei Dafu scattering image by utilizing an interferometry technology to obtain an interference phase map, and performing two-dimensional space unwrapping on the interference phase map by a phase unwrapping algorithm to obtain an interference phase of a pixel point of the Lei Dafu scattering image, wherein the interference phase of the pixel point can be expressed as:
φ＝φ def +φ atm +φ noise
wherein phi is the interference phase of the target point def Is the deformation phase of the target point phi atm Is the atmospheric phase phi noise Representing the noise phase. In the actual information processing process, the noise phase phi noise Effective cancellation can typically be performed using spatial filtering algorithms, and thus the interference phase of the target can be approximated as follows:
φ＝φ def +φ atm ；
s3, constructing a binary atmospheric parameter model containing a distance and an azimuth angle, and giving out a model parameter estimation method;
specifically, S3 includes the following steps:
s31, constructing a binary atmospheric parameter model containing a distance and an azimuth angle, wherein the binary atmospheric parameter model is shown in the following formula:
wherein phi is atm Represents the atmospheric error phase, lambda is the wavelength of a radar transmitting signal, r is the distance between the radar and a target, theta is the line-of-sight included angle between the target and the center of the radar aperture, and beta 0 Is a constant term, beta 1 Representing coefficients related to the first power of distance, beta 2 Representing coefficients related to the second power of distance, beta 3 Representing coefficients related to the azimuthal first power, beta 4 Representing coefficients related to the azimuthal second power, beta 5 Representing coefficients that are jointly related to the distance first power and the azimuth first power.
S32, for the stable PS point, the deformation phase can be approximately 0, so that the interference phase of the stable PS point is approximately the atmospheric phase, and if n PS points are shared in the detection area, a simultaneous equation system of the atmospheric phase, the distance and the azimuth angle can be established:
Φ＝Xβ+ε
wherein,β＝[β 0 ,β 1 ,…,β 5 ] T ，ε＝[ε 0 ,ε 1 ,…,ε n ] T . Phi is an n X1-dimensional vector formed by n PS-point interference phases, X is an n X6-dimensional vector formed by a constant 1 and n PS distances and azimuth angles, beta is a 6X 1-dimensional vector to be estimated, and epsilon is an n X1-dimensional vector formed by phase errors. Solving beta by least square:
the estimated value of the atmospheric phase is:
s4, based on two-dimensional complex scattering radar image sequence { Img } k And (3) performing multi-threshold screening to finish stable PS point extraction through amplitude dispersion, average correlation coefficient and relative amplitude mean value by using the PS technology.
Specifically, S4 includes the following steps:
s41, calculating a correlation coefficient gamma according to the main and auxiliary complex scattering images, wherein the formula is as follows:
wherein, γ is the coherence coefficient of the pixel point (i, j), x represents complex conjugate, and w1 and w2 are half of the sliding window length of two dimensions of the image matrix respectively.
Image sequence { Img } of length L k The correlation coefficient sequence matrix { gamma } with the length of L-1 can be calculated k (k=1, …, L-1), and thus an average correlation coefficient can be obtainedThe formula is as follows:
setting threshold T 1 Screening outIs used as a candidate point for the PS point.
S42, assume a radar radiation image Img k The dimension of (i) is M x N, and the complex scattering image value corresponding to the pixel point (i, j) can be expressed as
Wherein A is k (i, j) andrepresenting the amplitude and phase of the pixel (i, j) in the kth frame image, respectively.
Image sequence { Img } of length L k Average relative amplitude of each pixel point (i, j)Can be calculated by the following formula:
wherein,the average value of the amplitude of the pixel points (i, j) in the L frame sequence images is calculated by the following formula:
wherein,the average value of the amplitude mean value of all the pixels is calculated by the following formula:
setting threshold T 2 Screening out the candidate point set obtained in the step S42PS point of (c).
S43, calculating amplitude dispersion D of time sequence complex scattering image A The formula is as follows:
wherein delta A Standard deviation, m, representing pixel amplitude time sequence A The mean value of the pixel amplitude time series.
Setting threshold T 3 D is screened out from the candidate point set obtained in the step S42 A ≤T 3 PS point of (c).
S44, step S41 is performed to obtain threshold T 1 Set to 0.9, step S42 is performed to obtain the threshold T 2 Set to 5, step S43 is performed 3 Set to 0.1, the PS point selected at this time is taken asSample points of the atmospheric estimate;
s45, step S41, the threshold value T 1 Set to 0.8, step S42 is performed to obtain a threshold T 2 Set to 0, step S43 is performed 3 And setting the PS point to be 0.2, wherein the PS point screened at the moment is used as a monitoring point for deformation inversion.
S5, carrying out optimal estimation on parameters of a two-dimensional atmosphere model based on interference phases of stable PS points, calculating residual phases of the PS point interference phases and the atmosphere error phases, constructing a PS point rejection threshold based on statistical mean and variance of the residual phases, rejecting PS points with residual phases exceeding the threshold, repeating the atmosphere parameter estimation and PS rejection processes by using the residual PS points until iteration is terminated when the relative variation before and after two iterations of the residual phase threshold is smaller than a preset threshold, and obtaining the optimal atmosphere error phases under the optimal parameter estimation condition.
Specifically, S5 specifically includes the following steps:
s51, based on the atmospheric correction PS points screened by the threshold in the step S44, obtaining the atmospheric correction PS point interference phases from the interference phase diagram obtained in the step 2, and then performing two-dimensional space unwrapping on the interference phases by using an InSAR phase unwrapping algorithm (such as a branch-cut method, a quality-guided graph method, a minimum network cost flow method and the like);
s52, adopting the parameter model in the step S32 and the optimal solution method in the step S33 for the atmospheric correction PS point interference phase after unwrapping in the step S51, and obtaining the following results:
wherein m is the optimal calculation iteration sequence number of the atmospheric parameter, and corresponds toFor the mth iteration, calculate the atmospheric parameters, Φ m Correcting a subset of PS-point interference phases Φ, X for the atmosphere m Is phi m The corresponding coefficient matrix vector is also a subset of matrix X. When m=1, Φ 1 ＝Φ，/>X 1 ＝X。
Atmospheric parameters obtained by the mth iteration estimationThe atmospheric phase error can be calculated as:
s53, performing compensation processing on the interference phase of the PS point of the atmosphere estimation sample by using the atmosphere phase obtained in the step S52 to obtain an atmosphere compensation residual phaseThe method comprises the following steps:
sample PS point is estimated for atmosphere used for mth iteration estimationThe corresponding value in the (B) is subjected to statistics of sample mean and variance to obtain
S54, setting a phase residual error judgment threshold
Wherein T is 0 The threshold coefficient may be set to 2.
According to the thresholdRemoving the interference phase of the PS point of the atmospheric estimation sample, i.e. if +.>The PS point is needed to be removed from the atmospheric compensation sample points, the removing process is traversed through all the atmospheric estimation sample PS points, and the threshold relative change amount is calculated at the same time>/>
S55, whenIf so, repeating steps S52 to S54, otherwise, terminating the iteration, usingPhase diagram Img with interference inpha The deformation inversion can be completed, wherein the iteration threshold value is 0<N 0 Less than or equal to 1 and can be set to 0.95.
Comparing the atmospheric estimated phase curves of the strips in fig. 3 and fig. 4 at a certain distance, after the abnormal PS points are removed, the estimated atmospheric phase errors are more practical, i.e. the phase compensation errors of all PS participating in the atmospheric phase error estimation are the smallest in a statistical sense. In addition, by using the method for correcting the atmospheric pressure of the ground-based synthetic aperture slope radar, as shown in fig. 5, deformation errors of monitoring points (actually non-deformation points, namely, accumulated deformation amount should be 0 mm) caused by atmospheric pressure phase errors are well corrected.
The foregoing description of the preferred embodiments of the invention is not intended to be limiting, but rather is intended to cover all modifications, equivalents, and alternatives falling within the spirit and principles of the invention.
Claims (10)
1. The ground synthetic aperture slope radar atmosphere correction method is characterized by comprising the following steps of:
s1, acquiring a two-dimensional complex scattering radar image, performing interference processing on the radar image to obtain an interference phase diagram of the radar image, performing two-dimensional space unwrapping on the interference phase diagram to obtain an interference phase of a pixel point of the radar image, screening out a stable atmosphere estimated PS point from the pixel point of the radar image, and acquiring the interference phase of the stable atmosphere estimated PS point;
s2, setting a phase residual error judgment threshold based on an interference phase of a binary atmosphere model and the stable atmosphere estimation PS point, rejecting the stable atmosphere estimation PS point which does not meet the phase residual error judgment threshold, resetting the phase residual error judgment threshold based on the stable atmosphere estimation PS point which is not rejected, repeating the rejecting if the relative variation of the phase residual error judgment threshold meets an iteration condition, otherwise, ending the iteration, and finally taking the stable atmosphere estimation PS point which is not rejected as the PS point for optimal atmosphere compensation;
and S3, performing optimal atmosphere compensation based on the optimal atmosphere parameters of the PS points of the optimal atmosphere compensation and the interference phase diagram to finish deformation inversion.
2. The ground-based synthetic aperture slope radar atmosphere correction method according to claim 1, wherein acquiring the interference phase of the pixel point of the radar image comprises the steps of:
acquiring a two-dimensional complex scattering radar image, and marking the two-dimensional complex scattering radar image as Img k And generates a radar complex scattering image sequence { Img } k From the sequence { Img } k Extracting two adjacent frame images as main and auxiliary images of interference processing, taking complex conjugate of the auxiliary images and multiplying the auxiliary images with the main image, and obtaining an interference phase map Img after taking the main value phase inpha For the interference phase diagram Img inpha Two-dimensional space unwrapping is carried out to obtain the interference phase diagram Img inpha The interference phase phi of the pixel points of (2) is expressed as:
φ＝φ def +φ atm +φ noise ；
wherein phi is def Is the deformation phase phi of the pixel point atm Is the atmospheric error phase of the pixel point, phi noise Is the noise phase of the pixel.
3. The method for correcting the atmosphere of the ground-based synthetic aperture slope radar according to claim 2, wherein the step of screening out the stable atmosphere estimated PS points from the pixel points of the radar image comprises the steps of:
first, the sequence { Img } is calculated k Average correlation coefficient between adjacent frame imagesAnd then according to the radar radiation image Img k Is defined and the image Img k Amplitude A of the middle pixel point (i, j) k (i, j) and phase->Calculating the sequence { Img } k Average relative amplitude +.>Then, the sequence { Img is calculated k Standard deviation delta of amplitude time series of pixel point (i, j) in } A Mean value m from the amplitude time series A Is the amplitude dispersion D of the time sequence complex scattering image A The method comprises the steps of carrying out a first treatment on the surface of the Finally, respectively setting threshold value threshold T 1 、T 2 And T 3 Screening out pixels of the radar image while satisfying +.> And D is A ≤T 3 As the stable atmosphere estimated PS point, and simultaneously obtaining the interference phase of the stable atmosphere estimated PS point.
4. A method of ground-based synthetic aperture slope radar atmosphere correction as claimed in claim 3, wherein said threshold T 1 、T 2 And T 3 The settings of (2) are as follows: the threshold value threshold T 1 To set 0.9, T 2 Set to 5 and T 3 Setting the PS point to be 0.1, wherein the PS point finally screened out at the moment is used as a stable atmosphere estimated PS point; the threshold value threshold T 1 Set to 0.8, T 2 Set to 0 and T 3 And setting the PS point to be 0.2, wherein the PS point screened at the moment is used as a monitoring point for deformation inversion.
5. A method of ground-based synthetic aperture slope radar atmosphere correction according to claim 2, wherein the noise phase Φ in the interference phase Φ of the stable atmospheric estimate PS point is eliminated by a spatial filtering algorithm noise So that the interference phase phi of the stable atmosphere estimated PS point and the atmosphere error phase phi thereof atm Equal.
6. The method for correcting the atmosphere of the ground-based synthetic aperture slope radar according to claim 2, wherein S2 specifically comprises the steps of:
s21, constructing a binary atmosphere model based on the distance between the radar and the target and the included angle between the target and the radar aperture center, and calculating the atmosphere error phase phi of the stable atmosphere estimated PS point based on the model atm Thereby obtaining an interference phase phi of the stable atmosphere estimated PS point;
s22, establishing a simultaneous equation system of the interference phase, the distance and the azimuth angle of the stable atmospheric estimated PS point, and solving an atmospheric model parameter vector beta by adopting a least square method to obtain an optimized estimated value of the vector betaAnd calculates a vector X composed of the distance and azimuth angle of the stable atmospheric estimated PS point and the optimal estimated value +.>A vector product which is the atmospheric phase error estimated value phi atm ；
S23, based on the atmospheric phase error estimation value phi atm Compensating the interference phase of the stable atmosphere estimated PS point,obtaining the atmospheric compensation residual error phase phi sub And calculates the atmospheric compensation residual phase phi sub Is the mean of the samples of (2)And variance->Based on the sample mean->Variance->Phase residual decision threshold coefficient setting phase residual decision threshold +.>The calculation formula is as follows:
wherein T is 0 Is a threshold coefficient;
s24, according to the thresholdEliminating the stable atmospheric estimated PS point, if the absolute value of the interference phase of the stable atmospheric estimated PS point is +.>Removing;
s25, repeatedly executing the stable atmospheric estimation PS points which are not removed in the step S24 in the steps S21-S24 to obtain a phase residual error judgment threshold corresponding to the current iteration number mm is more than or equal to 2, and calculating the decision threshold +.>Relative amount of change of (2)If->When 0 < N 0 And (3) continuing iteration, otherwise, ending the iteration, and taking the stable atmospheric estimated PS points which are not rejected as PS points for optimal atmospheric compensation when the iteration is ended.
7. The method for correcting the atmosphere of a ground-based synthetic aperture slope radar according to claim 6, wherein the binary atmosphere model construction formula is:
wherein phi is atm Represents the atmospheric error phase, lambda is the wavelength of a radar transmitting signal, r is the distance between the radar and a target, theta is the line-of-sight included angle between the target and the center of the radar aperture, and beta 0 Is a constant term, beta 1 Representing coefficients related to the first power of distance, beta 2 Representing coefficients related to the second power of distance, beta 3 Representing coefficients related to the azimuthal first power, beta 4 Representing coefficients related to the azimuthal second power, beta 5 Representing coefficients that are jointly related to the distance first power and the azimuth first power.
8. The method of claim 6, wherein the system of equations in step S22 is:
Φ＝Xβ+ε
wherein Φ= [ phi ] 1 ，φ 2 ，…，φ n ] T ，β＝[β 0 ，β 1 ，…，β 5 ] T ，ε＝[ε 0 ，ε 1 ，…，ε n ] T Phi is an n X1-dimensional vector composed of interference phases of n PS points, X is an n X6-dimensional vector composed of distances and azimuth angles of constant 1 and n PS points, beta is a 6X 1-dimensional vector composed of atmospheric model parameters beta, and epsilon is an n X1-dimensional vector composed of phase errors.
9. The method for atmospheric correction of a ground-based synthetic aperture slope radar as defined in claim 6, wherein said threshold coefficient T in step S23 0 Set to 2.
10. The method for atmospheric correction of a ground-based synthetic aperture slope radar of claim 6, wherein said iteration threshold N 0 Set to 0.95.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310819068.1A CN117075056A (en) | 2023-07-05 | 2023-07-05 | Ground-based synthetic aperture slope radar atmosphere correction method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310819068.1A CN117075056A (en) | 2023-07-05 | 2023-07-05 | Ground-based synthetic aperture slope radar atmosphere correction method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117075056A true CN117075056A (en) | 2023-11-17 |
Family
ID=88714166
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310819068.1A Pending CN117075056A (en) | 2023-07-05 | 2023-07-05 | Ground-based synthetic aperture slope radar atmosphere correction method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117075056A (en) |
-
2023
- 2023-07-05 CN CN202310819068.1A patent/CN117075056A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107102333B (en) | Satellite-borne InSAR long and short baseline fusion unwrapping method | |
| CN112965146B (en) | Quantitative precipitation estimation method combining meteorological radar and rainfall barrel observation data | |
| CN113960595A (en) | Surface deformation monitoring method and system | |
| CN110780327B (en) | Marine target cooperative positioning method based on satellite-borne AIS and infrared camera | |
| CN112904337B (en) | Side slope deformation time sequence monitoring method based on Offset Tracking technology | |
| CN112284332B (en) | High-rise building settlement monitoring result three-dimensional positioning method based on high-resolution INSAR | |
| KR101804522B1 (en) | Apparatus and Method for SAR Offset Tracking using Multiple-Displacement estimated Kernel | |
| CN111856459B (en) | Improved DEM maximum likelihood constraint multi-baseline InSAR phase unwrapping method | |
| CN112444188B (en) | Multi-view InSAR sea wall high-precision three-dimensional deformation measurement method | |
| CN111220980B (en) | Ground-based SAR nonlinear atmospheric phase compensation method | |
| Izumi et al. | Iterative atmospheric phase screen compensation for near-real-time ground-based InSAR measurements over a mountainous slope | |
| CN114397629A (en) | Foundation arc interference synthetic aperture radar atmospheric interference phase correction method | |
| CN113009490A (en) | Radar three-dimensional wind field inversion method based on high-resolution mode dynamic constraint | |
| CN115079172A (en) | MTInSAR landslide monitoring method, equipment and storage medium | |
| Feng et al. | A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica | |
| CN108646244B (en) | Analysis method and system for measuring five-dimensional deformation of building | |
| CN117075056A (en) | Ground-based synthetic aperture slope radar atmosphere correction method | |
| JP4937791B2 (en) | Radar image processing device | |
| CN109946682B (en) | GF3 data baseline estimation method based on ICESat/GLAS | |
| CN115457022B (en) | Three-dimensional deformation detection method based on live-action three-dimensional model front-view image | |
| CN115343709B (en) | Terrain inversion method suitable for distributed spaceborne SAR system | |
| CN116485857A (en) | High-time-resolution glacier thickness inversion method based on multi-source remote sensing data | |
| CN110940984B (en) | Dual-polarization radar ratio differential phase shift rapid estimation method based on variational analysis | |
| CN114004870A (en) | Unmanned aerial vehicle SAR image geocoding method using optical orthographic image | |
| CN113568008A (en) | Sunflower-8 satellite-based minute-level precipitation real-time inversion estimation method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |