CN216411556U - Synthetic aperture radar system based on time series interference deformation measurement - Google Patents
Synthetic aperture radar system based on time series interference deformation measurement Download PDFInfo
- Publication number
- CN216411556U CN216411556U CN202120708621.0U CN202120708621U CN216411556U CN 216411556 U CN216411556 U CN 216411556U CN 202120708621 U CN202120708621 U CN 202120708621U CN 216411556 U CN216411556 U CN 216411556U
- Authority
- CN
- China
- Prior art keywords
- phase
- image
- deformation
- interference
- offset
- 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.)
- Active
Links
- 238000005259 measurement Methods 0.000 title claims abstract description 23
- 238000006243 chemical reaction Methods 0.000 claims description 11
- 238000000034 method Methods 0.000 abstract description 61
- 238000004364 calculation method Methods 0.000 abstract description 26
- 238000007781 pre-processing Methods 0.000 abstract description 5
- 230000009466 transformation Effects 0.000 abstract description 5
- 238000012544 monitoring process Methods 0.000 description 16
- 238000001914 filtration Methods 0.000 description 15
- 238000001228 spectrum Methods 0.000 description 14
- 238000010586 diagram Methods 0.000 description 9
- 230000006870 function Effects 0.000 description 8
- 238000012545 processing Methods 0.000 description 8
- 239000008186 active pharmaceutical agent Substances 0.000 description 7
- 238000012952 Resampling Methods 0.000 description 6
- 238000012937 correction Methods 0.000 description 5
- 238000013461 design Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 238000005457 optimization Methods 0.000 description 4
- 230000000737 periodic effect Effects 0.000 description 4
- 238000001514 detection method Methods 0.000 description 3
- 238000000819 phase cycle Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000001360 synchronised effect Effects 0.000 description 3
- 238000009825 accumulation Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 230000003595 spectral effect Effects 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000009436 residential construction Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
A synthetic aperture radar system based on time series interferometric deformation measurements, the synthetic aperture radar system comprising: carrying out data preprocessing on the SAR image set; removing flat ground and terrain phases from the interference phases to generate differential interference phases, and calculating pixel by pixel to generate a differential interference pattern; carrying out linear deformation phase estimation of time and space domains on the differential interference phase to obtain a time sequence deformation phase of each point target; and calculating phase transformation deformation according to the radar wavelength parameters, thereby obtaining a deformation measurement value of the synthetic aperture radar image. The method has the advantages of high deformation measurement precision, concise and efficient calculation steps and stable calculation result.
Description
Technical Field
The present invention relates to a Synthetic Aperture Radar (SAR), and more particularly, to a synthetic aperture radar system based on time series interferometric deformation measurements.
Background
Synthetic Aperture radar (sar), an active earth observation system, can be installed on flight platforms such as airplanes, satellites, spacecraft, etc., and can perform earth observation all day long and all day long, and has a certain ground surface penetration capability. Therefore, the SAR system has unique advantages in disaster monitoring, environmental monitoring, marine monitoring, resource exploration, crop estimation, mapping, military and other applications, and can play a role that other remote sensing means are difficult to play, so that the SAR system is more and more paid attention by countries in the world.
The satellite radar time sequence interference technology (TS-InSAR) is a remote sensing technology for detecting ground millimeter deformation in space of 800 kilometers by utilizing a satellite radar multi-temporal big data interference coherence principle, is the only surface monitoring technical means with high precision, all-time, all-weather, all-coverage and full-automatic at present, has an average detection range of more than 1500 square kilometers, can provide services such as accurate deformation measurement, state monitoring, abnormal early warning, post-disaster damage assessment and the like for the industries of China, land, traffic, surveying and mapping, railways, electric power, water conservancy, municipal administration, disaster reduction, residential construction, energy and the like, and efficiently solves the problem of pain points in the geographic information monitoring market.
The deformation monitoring and measuring technology has the characteristics of high precision, all-time, all-weather, full coverage, full automation and the like, but the technical problem that how to make the measured deformation more accurate is urgently needed to be solved at present due to the motion of the satellite and the space downward shooting visual angle.
SUMMERY OF THE UTILITY MODEL
In view of the above, the primary objective of the present invention is to provide a synthetic aperture radar system based on time-series interferometric deformation measurement, so as to at least partially solve the above technical problems.
In order to achieve the above object, the present invention provides a synthetic aperture radar system based on time series interferometric deformation measurement, comprising:
the antenna assembly comprises an active phase antenna array and a transceiving module, and is used for configuring the phase and the gain of a transmitting wave according to a beam position parameter transmitted by the central control assembly, receiving a radio frequency echo signal and performing down-conversion on the radio frequency echo signal to the central control assembly;
the central control component comprises a processor and a memory and is used for managing the distribution of task parameters and the time sequence of system operation; the central control assembly comprises a phase interference unit, and the phase interference unit is used for removing the flat land phase and the terrain phase from the interference phase of the interference pattern of the SAR image acquired by the antenna assembly, generating a differential interference phase and calculating to generate the differential interference pattern.
Based on the technical scheme, compared with the prior art, the synthetic aperture radar system based on the time series interference deformation measurement has at least one of the following beneficial effects:
the SAR system has high deformation measurement precision and concise and efficient calculation steps, can overcome the problem of time decorrelation caused by the increase of time span and the increase of Doppler centroid frequency difference of two images in the tiny surface deformation measurement, and has stable calculation result;
the SAR system can eliminate the influence of factors such as noise, shadow and the like, eliminates the phenomenon of mismatching in the matching result, and ensures the accuracy of image registration;
the SAR system can greatly improve the accuracy of image registration and the accuracy of deformation calculation through rough estimation and fine estimation of image offset;
the SAR system not only carries out linear deformation phase estimation of time and space domains on differential interference phases of the interferogram, but also carries out nonlinear deformation phase estimation when needed, so that residual phases such as atmosphere and noise can be removed, and the time sequence deformation phase of each point target is obtained;
the SAR system of the utility model can use the existing high-precision control point data (deformation amount of synchronous observation) such as level, GPS and the like to correct the benchmark, thereby ensuring the accuracy of calculation.
Drawings
FIG. 1 is a block schematic diagram of a time series interferometric deformation measurement based synthetic aperture radar system of the present invention;
FIG. 2 is a process flow diagram of data processing of the TS-InSAR technique of embodiment 1 of the present invention;
FIG. 3 is a schematic diagram of initial offset estimation of an orbit image of a multi-temporal image fine registration module according to the present invention;
FIG. 4 is a schematic diagram illustrating anchor point distribution in the rough image offset estimation step according to the present invention;
FIG. 5 is a schematic diagram of the main and auxiliary image selection for rough image offset estimation according to the present invention;
FIG. 6 is a flowchart of the image offset rough estimation algorithm of the present invention;
FIG. 7 is a schematic flow chart of the present invention for eliminating gross error points by mean variance method.
Detailed Description
In describing particular embodiments, specific details of structures, properties, effects, or other features are set forth in order to provide a thorough understanding of the embodiments by one skilled in the art. However, it is not excluded that a person skilled in the art may implement the utility model in a specific case without the above-described structures, performances, effects or other features.
The flow chart in the drawings is only an exemplary flow demonstration, and does not represent that all the contents, operations and steps in the flow chart are necessarily included in the scheme of the utility model, nor does it represent that the execution is necessarily performed in the order shown in the drawings. For example, some operations/steps in the flowcharts may be divided, some operations/steps may be combined or partially combined, and the like, and the execution sequence shown in the flowcharts may be changed according to actual situations without departing from the gist of the present invention.
The TS-InSAR technology is an interferometric method for calculating the time domain deformation and the space domain deformation of a permanent scatterer in a long-time sequence SAR image set so as to extract high-precision time sequence deformation information. The data processing process flow comprises processing flows of differential interference calculation, time and space deformation quantity estimation and the like, and can accurately restore the satellite measurement deformation.
Specifically, the utility model discloses a synthetic aperture radar system based on time series interference deformation measurement, which comprises:
the antenna assembly comprises an active phase antenna array and a transceiving module, and is used for configuring the phase and the gain of a transmitting wave according to a beam position parameter transmitted by the central control assembly, receiving a radio frequency echo signal and performing down-conversion on the radio frequency echo signal to the central control assembly;
the central control component comprises a processor and a memory and is used for managing the distribution of task parameters and the time sequence of system operation; the central control assembly comprises a phase interference unit, and the phase interference unit is used for removing the flat land phase and the terrain phase from the interference phase of the interference pattern of the SAR image acquired by the antenna assembly, generating a differential interference phase and calculating to generate the differential interference pattern.
Wherein the task parameter distribution managed by the central control component comprises: transmit pulse parameters and/or antenna beam settings.
The memory of the central control component is a read-only memory, such as an independent ROM memory, and the following execution logic is solidified therein:
step 1: performing data preprocessing on all images in the SAR image set so as to select a main image and an auxiliary image, performing pre-filtering on all the main and auxiliary images, and calculating an interference phase to generate an interference pattern;
step 2: removing flat ground and terrain phases from the interference phases of the interference pattern obtained in the step 1 to generate differential interference phases, and calculating pixel by pixel to generate a differential interference pattern;
and step 3: performing linear deformation phase estimation of time and space domains on the differential interference phase of the differential interference image obtained in the step 2 to obtain a time sequence deformation phase of each point target;
and 4, step 4: and calculating phase transformation deformation according to the radar wavelength parameters, thereby obtaining a deformation measurement value of the SAR system.
Wherein the step of selecting the main image specifically comprises:
a. calculating time and space baselines among all image pairs to generate a time and space baseline distribution diagram;
b. one image centered on the temporal and spatial baselines is selected as the primary image.
The step 1 of performing data preprocessing on all images in the SAR image set comprises the following steps:
and registering and cutting all SAR images relative to the main image, and combining to generate a time sequence interferogram set.
The specific steps of registering and cutting all SAR images relative to the main image are as follows:
A. registering all images relative to a main image, requiring the combined image pairs, registering according to the main image, and cutting all images into areas with consistent ranges;
B. cutting all data into areas with consistent ranges;
C. and (3) performing image pair combination on all the registered interference image pairs and the main image according to a time sequence, calculating an interference phase pixel by pixel, and generating a time sequence interference image set.
Wherein, the concrete requirements of the step A meet the following regulations:
selecting a registration algorithm, setting registration parameters, and performing registration calculation on each image pair in all images;
the errors of azimuth direction and distance direction are required to be less than 0.25 pixel when the main and auxiliary images are registered, and homonymous points of the calculated registration polynomial are uniformly distributed on the whole scene image.
Wherein, the concrete requirements of the step B meet the following regulations:
all the public areas of the cut registered images are larger than or equal to the designed monitoring working range, and if the public areas are missing, data can be supplemented in time;
and selecting a common area in the registered image as an InSAR processing range, and cutting all images into areas in the same range.
Wherein, the main image in the step C is selected by, for example, the following steps:
selecting a proper SAR image as a main image and other images as auxiliary images according to a three-dimensional space distribution map formed by a time base line, a space base line and Doppler center frequency of the SAR images obtained from the same area to form N interference images; on the basis of fully considering the optimization of three factors of time base line, space base line and Doppler centroid frequency difference among N +1 SAR images, the following comprehensive correlation coefficient model is established:
wherein,
in the above formula, γmIn order to synthesize the correlation coefficient, K is a constant, N is the number of SAR interferograms after registration,
and B⊥，cRespectively representing the corresponding space vertical baseline and time baseline of the main image and the auxiliary image, tk，mAnd tcRespectively representing the time accumulation amount of the relative reference time corresponding to the main image and the slave image, fk，mAnd fcRespectively representing the corresponding numbers of the master image and the slave imageThe frequency of the center of mass of the plerian, x is the position of the azimuth, and c is the speed of light; when the comprehensive correlation coefficient reaches the maximum, the corresponding image is the main image of the selected public area, which is called the public main image for short.
Wherein, the step A, B of registering and cutting all the SAR images with respect to the main image specifically includes:
calculating the initial offset of the auxiliary image relative to the main image based on the track information in the SLC parameter file;
the roughly estimating the image offset specifically includes: roughly estimating distance and azimuth registration offset by using cross-correlation optimization of the main SLC data and the auxiliary SLC data, and roughly estimating a registration offset polynomial according to an offset file;
the method for precisely estimating the image offset specifically includes: optimizing and finely estimating distance and azimuth registration offset by using cross-correlation of the primary SLC data and the secondary SLC data, and finely estimating a registration offset polynomial according to an offset file;
based on the offset polynomials in the distance and azimuth directions in the offset parameter file, performing SLC resampling by using two-dimensional SINC interpolation;
and registering the DEM with the main image, and cutting the DEM range to be consistent with the main image.
Wherein the cross-correlation optimized coarse estimation of the primary and secondary SLC data is preferably pixel-level.
The step of roughly estimating the image offset is implemented by using an anchor point matching method, for example.
Wherein the step of roughly estimating the image offset specifically comprises:
adopting a window-based search method, namely intercepting single-vision complex image data with the size of NxN on a main image as a reference window, intercepting single-vision complex image data with the size of MxM on an auxiliary window as a search window, then moving the search window on the auxiliary image, and determining corresponding homonymous points corresponding to anchor points uniformly arranged on the main image at certain intervals according to a preset criterion; the image offset rough estimation unit matches the image homonymous point to 1 pixel.
The predetermined criterion is, for example, to perform rough offset estimation by using a correlation coefficient method and a spectrum maximum value method.
The spectrum maximum value method specifically comprises the following steps: and seeking the frequency spectrum maximum value of the two complex images as a registration measure, wherein when the two complex images are accurately registered, the quality of the formed interference fringes is highest. The specific measurement criteria of this step, for example, require that the following relation be satisfied:
f＝FFT(R·S*)；
where f is the relative maximum spectrum, FFT () is the fast fourier transform (complex conjugate), and R and S are the complex values of the corresponding pixels of the primary and secondary images, respectively.
In general, the basic idea of spectral maxima is to define the ratio of the maximum modulus of the complex interferogram to the sum of the moduli of the other frequency components, even though the relative maximum spectrum is the merit function used in the registration process. The specific operation is to take a point to be matched from the main image, intercept the sub-image in a larger window (search window, M × M) around the main image, as S1, take the corresponding point at the same position in the auxiliary image, intercept the sub-image in a smaller window (matching window, N × N) around the auxiliary image, as S2, where M >2 × N, to ensure the accuracy of the calculation. Sliding S2 point by point in two dimensions within S1, and calculating the interferogram spectrum at that position once per sliding, i.e., F ═ FFT (S1', S2). Here, S1' represents a sub-image of S1 and S2. Therefore, the maximum spectrum method actually uses the phase information of the InSAR images for registration, i.e. the closer the phases of the two images are, the larger the corresponding spectral value is, which is also called phase correlation.
Wherein, the fine estimation step is preferably at sub-pixel level when the image offset is fine estimated.
Wherein the step of performing a fine estimation of the image offset specifically comprises:
the image offset fine estimation is registered to a sub-pixel level on the basis of the offset of the auxiliary image relative to the main image, which is obtained by the rough estimation, and requires 1/10 pixels, and the unit adopts an oversampling method and a correlation coefficient FFT interpolation method;
and interpolating K times of data of the reference window of the main image and the auxiliary image search window, and solving the offset by using the interpolated data as the reference and adopting a correlation coefficient method, a frequency spectrum maximum value method or a minimum average fluctuation function method so as to match the interpolated data to the sub-pixel level.
Wherein, for the mismatching points, for example, a mean square error method or a polynomial fitting method is adopted to eliminate.
And in the image interpolation resampling step, the offset of each pixel of the full image is fitted through the offset of the anchor point, a new auxiliary image is obtained by utilizing the interpolated auxiliary image, and the auxiliary image is resampled according to the offset.
The image interpolation resampling step is realized by a bilinear interpolation method or a four-point cubic convolution interpolation method.
The specific steps in the steps of registering and cutting the DEM and the main image are as follows:
a. sampling the DEM to be consistent with the resolution of the main image;
b. the DEM and the main image are registered, and the registration precision is superior to 0.5 pixel;
c. calculating and generating a conversion lookup table from the DEM coordinate system to the SAR image coordinate system according to the registration relational expression;
d. and converting the DEM into the SAR image coordinate system by utilizing a polynomial fitting algorithm according to the conversion lookup table to generate the DEM under the image coordinate system.
In step C, the specific steps of pre-filtering all the main and auxiliary images, calculating the interference phase, and generating the interferogram are as follows:
a. in a frequency domain, intercepting common frequency bands of the main image and the auxiliary image to carry out pre-filtering to generate the filtered main image and the filtered auxiliary image;
b. and carrying out complex conjugate multiplication on the pixel pairs of the main and auxiliary images which are subjected to pre-filtering to generate interference phase values, and calculating pixel by pixel to generate an interference pattern.
Wherein, the step of pre-filtering all the main and auxiliary images in the step C, calculating the interference phase and generating the interferogram may further include:
the method comprises the following steps of screening CS point targets for pixels of a time series interference pattern set:
a. identifying the CS point target of the SAR data by adopting an amplitude dispersion index method or a signal-to-noise ratio method;
b. identifying DS points according to a Smirnov detection method in a KS theory;
c. and extracting the point target meeting the requirement of the condition from the interference image set to generate an interference phase sequence of the CS/DS point target.
The step of removing the flat ground and the terrain phase in the step 2 specifically includes, for example:
calculating a flat ground phase according to the space baseline parameters and the earth ellipsoid parameters;
calculating a terrain phase by utilizing the DEM after registration;
and removing the flat ground phase and the terrain phase from the interference phase to generate a differential interference phase, and calculating pixel by pixel to generate a differential interference pattern.
Wherein, the step of removing the flat ground and the terrain phase in the step 2 may further comprise the following steps:
visually checking each difference interference pattern, if the residual interference fringe is more than half wavelength, calculating the residual phase of the space base line, and removing.
The specific steps of calculating and removing the spatial baseline residual phase are as follows:
a. carrying out space baseline rough estimation on the differential interference pattern by utilizing the quadric surface model to obtain a rough estimation phase of the space baseline; subtracting the coarse estimation phase from the differential phase in the differential interference pattern to obtain a residual phase;
b. estimating the residual phase by utilizing fast Fourier transform to obtain a residual baseline phase;
c. adding the coarse estimated phase of the spatial baseline to the phase of the residual baseline to obtain an improved phase of the spatial baseline;
d. and removing residual flat ground phases from the flat ground phases by using the corrected spatial baseline phases, and calculating to obtain corrected flat ground phases and an interference pattern set.
Wherein, the step of removing the flat ground and the terrain phase in the step 2 may further comprise the following steps:
and performing linear deformation phase estimation of time and space domains on the differential interference phase of the interference pattern, performing nonlinear deformation phase estimation if required, and removing atmospheric or noise residual phases to obtain the time sequence deformation phase of each point target.
The specific calculation steps of the linear deformation phase estimation and the nonlinear deformation phase estimation are as follows:
firstly, CS point targets are connected to form a Delauney irregular triangular network (DTIN, or redundant network), and the difference phase difference of adjacent points is solved according to the connection relation between the points;
establishing a two-dimensional periodic chart of a CS point target according to the relation of a space baseline and a time baseline, maximizing a model correlation coefficient by taking the two-dimensional periodic chart as an objective function, and estimating a linear deformation rate and an elevation difference value between adjacent points; if the monitoring work design book only requires a linear deformation result, directly outputting the calculation of the vertical deformation amount of the result to generate a geological disaster body rate map;
removing two phase quantities in the step II from the differential interference phase to obtain a residual phase; carrying out spatial domain mean value filtering on the residual phase, and calculating to obtain a main image atmospheric phase; carrying out space domain low-pass filtering and time domain high-pass filtering on the residual phase without the main image atmospheric phase to obtain a slave image atmospheric phase, and further decomposing a nonlinear deformation phase;
and fourthly, adding the linear deformation phase obtained in the step III and the nonlinear deformation phase obtained in the step III, and combining time base line parameters to obtain the time sequence deformation phase of each CS point target.
And 4, converting the unwrapping phase into a deformation quantity in the sight line direction by combining external auxiliary data according to the radar wavelength parameter and after phase transformation deformation calculation, and converting the sight line deformation quantity into a vertical direction according to an included angle between the sight line and the vertical direction.
The method comprises a step of geocoding after phase transition shape change calculation in the step 4, wherein the geocoding step is to perform geocoding by using a Digital Elevation Model (DEM) product, and the method comprises the following specific steps:
a. utilizing a conversion lookup table from a DEM coordinate system to an SAR image coordinate system to complete inverse conversion of the monitoring result from the SAR image coordinate system to a geodetic coordinate system, namely geocoding is carried out on the vertical deformation of the monitoring result;
b. and (4) collecting all point targets after the geocoding, converting time units of the deformation into adults, generating annual deformation rate, and calculating pixel by pixel to generate a geological disaster body rate map.
And step 4, after phase transformation shape change calculation, a step of reference correction is further included, and the disaster body rate of the point target after geocoding is used for correcting the reference by utilizing the existing level and GPS high-precision control point data.
The reference correction comprises the following specific steps:
a. taking the synchronous leveling measurement result as a reference, and calculating the average value of the difference between the target deformation amount and the measured value of the point on the adjacent point, namely the integral deviation value between the target deformation amount and the measured value;
b. and adding the integral deviation value obtained in the last step into the deformation value of each point target, and correcting the integral deviation of InSAR result deformation caused by non-uniform reference points to finish the benchmark correction.
In order that the objects, technical solutions and advantages of the present invention will become more apparent, the present invention will be further described in detail with reference to the accompanying drawings in conjunction with the following specific embodiments.
As shown in fig. 2, the synthetic aperture radar system based on time-series interferometric deformation measurement of the present embodiment mainly includes an antenna module and a central control module, and is configured to perform processing procedures such as data preprocessing, differential interference calculation, time and space deformation estimation, and deformation calculation. As will be described in greater detail below.
Data preprocessing
A method for selecting a main image
The TS-InSAR method of this embodiment selects a single main image. On the premise of meeting the requirements of a space baseline and a time baseline, the working steps of SAR main image selection and image pair combination are as follows:
a. calculating time and space baselines among all image pairs to generate a time and space baseline distribution diagram;
b. a scene with centered temporal and spatial baselines is selected as the main image.
Image registration, cutting and combination
And all SAR images carry out registration and cutting on the main image, and are combined to generate a time series interferogram set. The method comprises the following specific steps:
a. all images are registered to the main image. And registering the combined image pairs according to the main image, and cutting all the images into areas with consistent ranges. The following specifications should be met in particular:
and selecting a registration algorithm, setting registration parameters, and performing registration calculation on each image pair.
The errors of azimuth direction and distance direction are required to be less than 0.25 pixel when the main and auxiliary images are registered, and homonymous points of the calculated registration polynomial are uniformly distributed on the whole scene image.
b. All data is cropped to areas of consistent extent. The cutting requirements should be satisfied:
all the clipped public areas of the registered images are larger than or equal to the designed monitoring working range, and if the public areas are missing, data can be supplemented in time.
And selecting a common area in the registered image as an InSAR processing range, and cutting all images into areas in the same range.
c. And respectively carrying out image pair combination on all the registered interference image pairs and the main image according to a time sequence. And calculating interference phases pixel by pixel to generate a time sequence interference pattern set. The theory of registration is as follows:
master image selection
For N +1 SAR images acquired in the same area, a proper SAR image is selected as a main image according to a three-dimensional space distribution diagram consisting of a time base line, a space base line and a Doppler center frequency of the SAR images, and other images are used as auxiliary images to form N interferograms. In order to obtain interference phase with higher quality, before interference processing, an image with optimal space-time distribution is selected as a common main image of interference registration. The selection of the common main image relates to the correlation of the interference pair formed subsequently, and the time base line, the effective space base line and the Doppler centroid frequency difference are important factors influencing the interference correlation. For the measurement of interference deformation, the smaller the length of the effective space baseline is, the better the effective space baseline is; for the measurement of the tiny surface deformation, the image must have a longer time span, and the increase of the time span easily causes the time decorrelation; in addition, the larger the difference of the Doppler centroid frequencies of the two images is, the larger the decorrelation degree is. Therefore, on the basis of fully considering the optimization of three factors of time base line, space base line and Doppler centroid frequency difference among the N +1 SAR images, the following comprehensive correlation coefficient model can be established:
wherein,
in the above formula, γmIn order to synthesize the correlation coefficient, K is a constant, N is the number of SAR interferograms after registration,
and B⊥，cRespectively representing the corresponding space vertical baseline and time baseline of the main image and the auxiliary image, tk，mAnd tcRespectively representing the time accumulation amount of the relative reference time corresponding to the main image and the slave image, fk，mAnd fcRespectively representing the Doppler centroid frequency corresponding to the main image and the slave image, wherein x is the azimuth position, and c is the light speed; when the comprehensive correlation coefficient reaches the maximum, the corresponding image is the main image of the selected public area, which is called the public main image for short.
Image offset initialization
Fig. 3 is a schematic diagram of the estimation of the initial offset of the orbit image of the multi-temporal image fine registration module according to the present invention, and as shown in fig. 3, the image offset initialization unit calculates the initial offset of the secondary image relative to the primary image based on the orbit information in the SLC parameter file.
Coarse estimation of image offset
Description of the function:
the cross-correlation optimization of the primary and secondary SLC data is used to approximate (pixel-level) the distance and orientation registration offsets.
And (2) roughly estimating a registration offset polynomial according to the offset file.
Detailed design:
the purpose of the image offset rough estimation unit is to match the image homonym points to 1 pixel, and the module is intended to adopt an anchor point matching method, i.e. anchor points are uniformly arranged on the main image at certain intervals, and then the homonym points are searched on the auxiliary image, as shown in fig. 4.
In this embodiment, a window-based search method is adopted, that is, single-view complex image data (reference window) with a size of N × N is captured on a main image, single-view complex image data (search window) with a size of M × M is captured on an auxiliary window, the main window is moved on the auxiliary image, and a corresponding homonymy point of the complex image is determined according to a certain criterion, for example, the correlation coefficient is maximum, as shown in fig. 5.
In this embodiment, coarse estimation of image offset is performed by using a correlation coefficient method and a spectrum maximum value method.
And the frequency spectrum maximum value measure seeks the frequency spectrum maximum value of the two complex images as a registration measure, and when the two complex images are accurately registered, the formed interference fringe has the highest quality. The specific calculation formula is as follows:
f＝FFT(R·S*)；
where f is the relative maximum spectrum, FFT () is the fast fourier transform (complex conjugate), and R and S are the complex values of the corresponding pixels of the primary and secondary images, respectively.
The flow of estimation by the principle of maximum correlation coefficient is shown in fig. 6.
Image offset fine estimation
Description of the function:
optimizing and finely estimating distance and azimuth registration offset by using cross-correlation of main SLC data and auxiliary SLC data (sub-pixel level);
and (2) finely estimating a registration offset polynomial according to the offset file.
Detailed design:
the image offset fine estimation is registered to a sub-pixel level on the basis of the offset of the auxiliary image relative to the main image, which is obtained by rough estimation, wherein 1/10 pixels are generally required, and the unit adopts an oversampling method and a correlation coefficient FFT interpolation method.
And interpolating K times of data of a reference window of the main image and a search window of the auxiliary image, solving the offset by using the interpolated data as a reference and adopting a correlation coefficient method, a frequency spectrum maximum value method or a minimum average fluctuation function method, matching the offset to a sub-pixel level, and ensuring that the size of the reference window is a power of 2 when the frequency spectrum maximum value method is adopted. Since the coarse registration is already one pixel, it can be considered that the search window size is 4 more than the reference window size, when the search range is [ -2,2 ].
Assuming the offset of the auxiliary image relative to the main image during coarse registration, the position of the anchor point of the main image, and the offset of the image point calculated by the oversampling method, the coordinates of the anchor point of the main image on the auxiliary image after fine registration are:
the flow of the oversampling registration is basically similar to the coarse registration flow, and only the offset is solved by three coarse registration methods after the original window data is interpolated.
Due to the influence of factors such as noise, shadow and the like, some points in the matching result have a mismatching phenomenon, and the mismatching points must be detected to ensure the accuracy of image registration. The gross error rejection method comprises a mean variance method and a polynomial fitting method.
According to the formula of InSAR,
the pixel offset is
From theoretical analysis, it can be seen that the variance of the offset of the anchor point (the distance to the pixel coordinate) in the same row on the primary image should be better than 1 pixel on the secondary image. Therefore, the degree of deviation from the mean value can be used as a basis for eliminating the rough difference points. Meanwhile, the whole-column anchor control is based on the fact that the variance Δ is smaller than a certain value L.
The flowchart of the mean variance method for eliminating the coarse difference points is shown in fig. 7.
Image interpolation resampling
Description of the function:
the image interpolation module performs SLC resampling by using two-dimensional SINC interpolation based on the offset polynomial of the distance and the azimuth in the offset parameter file.
Detailed design:
the basic idea of the image interpolation resampling unit is to fit the offset of each pixel of the full image through the offset of the anchor point, and then obtain a new auxiliary image by using an interpolated auxiliary image. Because the offsets of the auxiliary images corresponding to all the pixels in the main image cannot be completely the same and the offset of each pixel cannot be found through registration, a second-order or multi-order polynomial is used for fitting the relationship between the image coordinates and the offset of the anchor pixel of the main image to obtain a polynomial coefficient in the pixels of which the offset is obtained through limited registration in the main image. By using the polynomial coefficients, the offset of each pixel in the main image can be calculated, and then the auxiliary image is resampled according to the offset.
After obtaining the exact offsets for the N anchor points, a second order (or higher order, here second order is taken as an example) polynomial fitting parameter is used:
and calculating an offset for the coordinates of each pixel in the main image by using the fitting parameters obtained by the formula, and finding the matched position in the auxiliary image according to the offset.
When the accurate offset between the two complex images is determined, the secondary image needs to be resampled (interpolated) so that the pixels in the secondary image are aligned with the pixels in the main image one by one. Since the nearest neighbor interpolation method cannot meet the requirement of interpolation precision for neutron pixel registration in InSAR interference processing, the bilinear interpolation method and the four-point cubic convolution interpolation method are adopted in the embodiment.
Registering and cropping a Digital Elevation Model (DEM) with the main image
And registering the DEM with the main image, and cutting the DEM range to be consistent with the main image. The method comprises the following specific steps:
a. DEM should be sampled to the resolution consistent with the main image;
b. the DEM and the main image are registered, and the registration precision is superior to 0.5 pixel;
c. calculating and generating a conversion lookup table from the DEM coordinate system to the SAR image coordinate system according to the registration relational expression;
d. and converting the DEM into the SAR image coordinate system by utilizing a polynomial fitting algorithm according to the conversion lookup table to generate the DEM under the image coordinate system.
Fourth interferogram phase calculation
And pre-filtering all the main and auxiliary images, and calculating interference phases to generate an interference pattern. The method comprises the following specific steps:
a. and (4) pre-filtering. In a frequency domain, intercepting common frequency bands of the main image and the auxiliary image to carry out pre-filtering to generate the filtered main image and the filtered auxiliary image;
b. and (4) interference phase calculation. And carrying out complex conjugate multiplication on the pixel pairs of the main and auxiliary images which are subjected to pre-filtering to generate interference phase values, and calculating pixel by pixel to generate an interference pattern.
Fifthly, selecting targets of CS point and DS point
And (3) carrying out CS point target screening on the pixels of the time series interference pattern set. The method comprises the following specific steps:
CS point target identification. The SAR data CS point target is preferably identified by methods such as an amplitude dispersion exponential method, a signal-to-noise ratio method and the like. One or more methods are preferably selected by combining the types of the ground objects in the monitoring area so as to improve the accuracy of CS point target identification;
and b, identifying DS point targets. Identifying DS points according to a Smirnov detection method in a KS theory;
and c, generating a CS/DS point target interference phase sequence. And extracting the point target meeting the requirement of the condition from the interference image set to generate an interference phase sequence of the CS/DS point target.
(II) differential interference calculation
First level land and space phase removal
Calculating a flat ground phase according to the space baseline parameters and the earth ellipsoid parameters; and calculating the terrain phase by utilizing the DEM after registration. And removing the flat ground phase and the terrain phase from the interference phase to generate a differential interference phase, and calculating pixel by pixel to generate a differential interference pattern.
A novel space baseline refinement
Visually checking each difference interference pattern, if the residual interference fringe is more than half wavelength, calculating the residual phase of the space base line, and removing. The method comprises the following specific steps:
a. carrying out space baseline rough estimation on the differential interference pattern by utilizing the quadric surface model to obtain a rough estimation phase of the space baseline; subtracting the coarse estimation phase from the differential phase in the differential interference pattern to obtain a residual phase;
b. estimating the residual phase by utilizing fast Fourier transform to obtain a residual baseline phase;
c. adding the coarse estimated phase of the spatial baseline to the phase of the residual baseline to obtain an improved phase of the spatial baseline;
d. and removing residual flat ground phases from the flat ground phases by using the corrected spatial baseline phases, and calculating to obtain corrected flat ground phases and an interference pattern set.
Estimation of the deformation in the time/space domain
And linear deformation phase estimation of time and space domains is carried out on the differential interference phase of the interference pattern, if required, nonlinear deformation phase estimation is carried out, residual phases such as atmosphere and noise are removed, and the time sequence deformation phase of each point target is obtained. The calculation steps of the TS-InSAR are as follows:
estimation of parameters between adjacent points
And connecting the CS point targets to form a Delauney irregular triangular network (DTIN, or redundant network), and solving the difference phase difference of adjacent points according to the connection relation between the points.
Linear deformation phase and residual elevation phase calculation
And establishing a two-dimensional periodic chart of the CS point target according to the relation between the space baseline and the time baseline, maximizing the correlation coefficient of the model by taking the two-dimensional periodic chart as an objective function, and estimating the linear deformation rate and the elevation difference between adjacent points. If the monitoring work design book only requires a linear deformation result, the calculation of the vertical deformation amount of the result can be directly output, and a geological disaster body velocity map is generated.
Computing nonlinear deformation phase and atmospheric phase
Two phase quantities in the second step are removed from the differential interference phase to obtain a residual phase. And carrying out spatial domain mean value filtering on the residual phase, and calculating to obtain the atmospheric phase of the main image. And performing space domain low-pass filtering and time domain high-pass filtering on the residual phase without the main image atmospheric phase to obtain a slave image atmospheric phase, and further decomposing a nonlinear deformation phase.
Fourthly, calculating the time series deformation phase
And (4) adding the linear deformation phase in the step (II) and the nonlinear deformation phase in the step (III), and combining time base line parameters to obtain the time sequence deformation phase of each CS point target.
Fourth deformation variable calculation
Firstly, calculating deformation quantity of visual line and converting vertical direction
According to the radar wavelength parameters, after phase transformation deformation calculation, according to the requirement, combining with external auxiliary data, converting the unwrapping phase into the deformation of the sight line direction, and then converting the sight line direction deformation into the vertical direction according to the included angle between the sight line and the vertical direction.
② geocoding
The geocoding method can utilize DEM products for geocoding. The method comprises the following specific steps:
a. and (4) utilizing a conversion lookup table from the DEM coordinate system to the SAR image coordinate system to complete the inverse conversion of the monitoring result from the SAR image coordinate system to the geodetic coordinate system, namely geocoding the vertical deformation of the monitoring result.
b. And (4) collecting all point targets after the geocoding, converting time units of the deformation into adults, generating annual deformation rate, and calculating pixel by pixel to generate a geological disaster body rate map.
Third, benchmark correction
The disaster body rate of the point target after geocoding can utilize the existing high-precision control point data (deformation amount of simultaneous observation) such as level, GPS and the like to correct the benchmark, and the specific steps are as follows:
a. taking the synchronous leveling measurement result as a reference, and calculating the average value of the difference between the target deformation amount and the measured value of the point on the adjacent point, namely the integral deviation value between the target deformation amount and the measured value;
b. and adding the integral deviation value obtained in the last step into the deformation value of each point target, and correcting the integral deviation of InSAR result deformation caused by non-uniform reference points to finish the benchmark correction.
The above-mentioned embodiments are intended to illustrate the objects, technical solutions and advantages of the present invention in further detail, and it should be understood that the above-mentioned embodiments are only exemplary embodiments of the present invention and are not intended to limit the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (1)
1. A synthetic aperture radar system based on time series interferometric deformation measurements, comprising:
the antenna assembly comprises an active phase antenna array and a transceiving module, and is used for configuring the phase and the gain of a transmitting wave according to a beam position parameter transmitted by the central control assembly, receiving a radio frequency echo signal and performing down-conversion on the radio frequency echo signal to the central control assembly;
the central control component comprises a processor and a memory and is used for managing the distribution of task parameters and the time sequence of system operation; the central control assembly comprises a phase interference unit which removes flat ground and terrain phases from interference phases of an interference pattern of the SAR image acquired by the antenna assembly, generates a differential interference phase and calculates and generates the differential interference pattern; the memory is a read-only memory, and execution logic is solidified in the memory; wherein the task parameter distribution managed by the central control component comprises: transmit pulse parameters and/or antenna beam settings.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202120708621.0U CN216411556U (en) | 2021-04-08 | 2021-04-08 | Synthetic aperture radar system based on time series interference deformation measurement |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202120708621.0U CN216411556U (en) | 2021-04-08 | 2021-04-08 | Synthetic aperture radar system based on time series interference deformation measurement |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN216411556U true CN216411556U (en) | 2022-04-29 |
Family
ID=81280323
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202120708621.0U Active CN216411556U (en) | 2021-04-08 | 2021-04-08 | Synthetic aperture radar system based on time series interference deformation measurement |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN216411556U (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115494501A (en) * | 2022-09-29 | 2022-12-20 | 南宁市勘测设计院集团有限公司 | Urban infrastructure deformation monitoring method based on high-resolution PS-InSAR |
| CN116148855A (en) * | 2023-04-04 | 2023-05-23 | 之江实验室 | Method and system for removing atmospheric phase and resolving deformation of time sequence InSAR |
-
2021
- 2021-04-08 CN CN202120708621.0U patent/CN216411556U/en active Active
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115494501A (en) * | 2022-09-29 | 2022-12-20 | 南宁市勘测设计院集团有限公司 | Urban infrastructure deformation monitoring method based on high-resolution PS-InSAR |
| CN116148855A (en) * | 2023-04-04 | 2023-05-23 | 之江实验室 | Method and system for removing atmospheric phase and resolving deformation of time sequence InSAR |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113340191B (en) | Time series interference SAR deformation quantity measuring method and SAR system | |
| CN112986993B (en) | InSAR deformation monitoring method based on space constraint | |
| Lanari et al. | Generation of digital elevation models by using SIR-C/X-SAR multifrequency two-pass interferometry: The Etna case study | |
| US6046695A (en) | Phase gradient auto-focus for SAR images | |
| US5923278A (en) | Global phase unwrapping of interferograms | |
| CN101770027B (en) | Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion | |
| US6011505A (en) | Terrain elevation measurement by interferometric synthetic aperture radar (IFSAR) | |
| CN109375222B (en) | Synthetic aperture radar interferometric ionosphere phase estimation and compensation method | |
| CN107102333B (en) | Satellite-borne InSAR long and short baseline fusion unwrapping method | |
| CN104316920A (en) | High-precision sea surface height extracting method of radar altimeter through small incidence angle interference | |
| CN216411556U (en) | Synthetic aperture radar system based on time series interference deformation measurement | |
| CN112051571A (en) | LOS (line of sight) deformation variable estimation method of novel differential InSAR (interferometric synthetic Aperture Radar) | |
| CN107918127A (en) | A kind of road slope deformation detecting system and method based on vehicle-mounted InSAR | |
| CN112711021B (en) | Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method | |
| CN105180852B (en) | GB SAR deformation monitoring methods based on triple steppings | |
| CN109239710B (en) | Method and device for acquiring radar elevation information and computer-readable storage medium | |
| CN115201825B (en) | Atmospheric delay correction method in InSAR (interferometric synthetic aperture radar) inter-seismic deformation monitoring | |
| CN110441770A (en) | Three-dimensional deformation measurement method based on multi-section MIMO-SAR joint observation | |
| CN109116321A (en) | A kind of phase filtering method and height measurement method of spaceborne interference imaging altimeter | |
| CN111721241A (en) | GNSS-InBSAR and GB-InSAR cross-system fusion three-dimensional deformation measurement method | |
| CN109085556B (en) | High-frequency ground wave radar wave field forming method based on first-order and second-order peak ratios | |
| Rossi et al. | High-resolution InSAR building layovers detection and exploitation | |
| Mao et al. | Estimation and compensation of ionospheric phase delay for multi-aperture InSAR: An azimuth split-spectrum interferometry approach | |
| WO2005010556A1 (en) | Radar position and movement measurement for geophysical monitoring | |
| CN117073621A (en) | Urban settlement monitoring method integrating InSAR and Beidou |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| GR01 | Patent grant | ||
| GR01 | Patent grant |