CN112882030B - InSAR imaging interference integrated processing method - Google Patents
InSAR imaging interference integrated processing method Download PDFInfo
- Publication number
- CN112882030B CN112882030B CN202110039286.4A CN202110039286A CN112882030B CN 112882030 B CN112882030 B CN 112882030B CN 202110039286 A CN202110039286 A CN 202110039286A CN 112882030 B CN112882030 B CN 112882030B
- Authority
- CN
- China
- Prior art keywords
- imaging
- interference
- phase
- antenna
- insar
- 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
- 238000003384 imaging method Methods 0.000 title claims abstract description 144
- 238000003672 processing method Methods 0.000 title claims abstract description 21
- 238000012545 processing Methods 0.000 claims abstract description 28
- 238000000034 method Methods 0.000 claims abstract description 23
- 238000009825 accumulation Methods 0.000 claims abstract description 13
- 230000001427 coherent effect Effects 0.000 claims abstract description 13
- 238000001914 filtration Methods 0.000 claims abstract description 13
- 238000005259 measurement Methods 0.000 claims abstract description 9
- 238000010587 phase diagram Methods 0.000 claims abstract description 9
- 238000004422 calculation algorithm Methods 0.000 description 17
- 230000008569 process Effects 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 238000013507 mapping Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000012876 topography Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 238000004804 winding Methods 0.000 description 1
Images
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
- 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
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
-
- 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/9094—Theoretical aspects
-
- 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/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
- G01S7/418—Theoretical aspects
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 provides an InSAR imaging interference integrated processing method, which comprises the following steps: positioning the geographic space of the image in an absolute geographic space coordinate system, and constructing a back projection imaging model for integrating InSAR imaging and interferometric height measurement; carrying out coherent accumulation imaging by using a back projection imaging model, outputting a digital orthophoto map, and carrying out conjugate multiplication on the digital orthophoto map to obtain a residual interference phase; carrying out phase filtering on the residual interference phase and outputting an interference phase diagram; performing ultra-sparse control point area network combined interference calibration according to the interference phase diagram and DEM auxiliary conditions, and calculating interference calibration result parameters; and establishing an elevation inversion model based on the back projection imaging model, and inverting the ground elevation in the elevation inversion model by using the interference calibration result parameters. The method has the remarkable advantages that InSAR interference processing is simpler, self-registration is completed during imaging, DEM and SAR data matching is avoided, multi-scene image splicing is avoided, and an output product is more user-friendly.
Description
Technical Field
The invention relates to the technical field of synthetic aperture radars, in particular to an InSAR imaging interference integrated processing method.
Background
Synthetic Aperture Radars (SAR) synthesize a large antenna array by platform motion, thereby greatly improving the transverse resolution and realizing two-dimensional high-resolution imaging to the ground. The interferometric Synthetic Aperture Radar (InSAR) combines the SAR technology and the interferometric technology, and utilizes the difference of target observation visual angles to invert the ground elevation information through interferometric phase difference. Due to the advantages of all-time and all-weather, the SAR/InSAR has wide application potential in the fields of surveying and mapping, homeland resource investigation, ocean monitoring, forestry and the like.
With the development of SAR/InSAR towards multi-mode, high resolution, wide swath, miniaturization, intelligent imaging and detection, the traditional SAR/InSAR imaging and processing algorithm can not meet the requirements of new system, new mode SAR imaging resolution and high precision. A ground Digital Elevation Model (DEM) represents the topographic relief degree, and has non-negligible influence on the high-resolution multimode SAR/InSAR motion compensation precision, the imaging resolution and the interference processing precision.
First, the influence of terrain fluctuation on SAR/InSAR imaging accuracy is seen. The traditional SAR imaging is based on the flat ground assumption, residual motion errors caused by areas with low and medium resolution or little topographic relief can be ignored, but for multi-mode SAR with high resolution, large synthetic aperture, large track offset, short acting distance and the like, the flat ground assumption causes large motion errors, thereby affecting SAR imaging accuracy.
In terms of DEM-assisted SAR imaging algorithm, traditional SATA, PTA, FD and the like perform terrain-dependent compensation by a sub-aperture dividing method, and terrain elevation information in the whole sub-aperture irradiation range is replaced by reference terrain in the sub-aperture range or in the beam center, so that terrain-assisted SAR imaging is realized to a certain extent.
Then, the influence of the topography on the accuracy of the interference processing is observed. The basic principle of InSAR height measurement is that ground elevation information is inverted through high-precision interference phases, single-vision complex images obtained by double antennas or two times of observation are subjected to the steps of registration, filtering, land removal, phase expansion and the like, unwrapped interference phase values are output, interference parameter calibration is solved by using ground control point information, and finally the ground elevation value is calculated according to an InSAR height inversion model. Therefore, the processing efficiency and accuracy of each step of the interference processing are influenced by the topographic relief, and the interference processing performance is obviously improved under the initial DEM auxiliary condition.
Aiming at the aspect of an InSAR processing algorithm assisted by DEM, the initial DEM information is used for estimating parameters such as an interference baseline, interference phase offset and the like, the density degree of phase winding is reduced in the phase unwrapping processing step, the DEM information is separately used for the InSAR processing process, and the efficiency and the precision are improved to a certain extent.
The traditional SAR imaging algorithm assisted by DEM is generally based on a two-dimensional frequency domain imaging algorithm, such as CS, RD and omega-k imaging algorithms, the approximation of a classical orientation space-variant motion compensation frequency domain algorithm signal model assisted by terrain is considered, motion errors caused by terrain compensation in frequency domain must be assumed, terrain elevation in a sub-aperture range or a beam center is replaced by reference terrain in the whole sub-aperture irradiation range, and the accuracy is difficult to meet the requirement in application scenes of high resolution, long synthetic aperture, large track offset and the like.
A traditional DEM auxiliary InSAR processing algorithm is separately used for interference baseline estimation, interference phase bias parameter estimation and phase unwrapping and terrain phase unwrapping. However, the auxiliary DEM information does not cover interference SAR processing steps such as registration, filtering and interference calibration, systematic auxiliary registration, filtering, phase expansion, interference calibration, elevation inversion and other full-flow interference processing of external DEM data are not seen, and the auxiliary DEM processing potential is not fully exerted to a certain extent by the aid of the independent auxiliary InSAR processing method.
Disclosure of Invention
Technical problem to be solved
In view of the above, the present invention provides an interferometric synthetic aperture radar (InSAR) imaging integration processing method to solve at least one of the above technical problems.
(II) technical scheme
The invention provides an InSAR imaging interference integrated processing method, which comprises the following steps: positioning the geographic space of the image in an absolute geographic space coordinate system, and constructing a back projection imaging model for integrating InSAR imaging and interferometric height measurement; carrying out coherent accumulation imaging by using the back projection imaging model, outputting a digital orthophoto map, and carrying out conjugate multiplication on the digital orthophoto map to obtain a residual interference phase; carrying out phase filtering on the residual interference phase and outputting an interference phase diagram; performing ultra-sparse control point area network combined interference calibration according to the interference phase diagram and DEM auxiliary conditions, and calculating interference calibration result parameters; and establishing an elevation inversion model based on the backward projection imaging model, and inverting the ground elevation in the elevation inversion model by using the interference calibration result parameters.
In some embodiments, the absolute geospatial coordinate system is a WGS84 spatial rectangular coordinate system, the locating the geospatial of the image comprising: and converting the external DEM elevation value into the external DEM elevation under the WGS84 space rectangular coordinate system through geocoding.
In some embodiments, the back projection imaging model comprises a primary antenna and a secondary antenna, and the performing coherent accumulation imaging with the back projection imaging model comprises: in InSAR imaging, the phase item of the crossing moment of the center of the azimuth beam is reserved, and the imaging signal models corresponding to the main antenna and the auxiliary antenna are respectively
Wherein,
respectively calculating the slant distances of the main antenna and the auxiliary antenna at any time according to the external DEM;
the main antenna and the auxiliary antenna which are respectively at the time of passing through the center of the azimuth beam are calculated to obtain slant distances according to the external DEM;
and
compressing the contribution of the image to the current pixel point for the distance of the main antenna and the auxiliary antenna at the crossing moment of the azimuth beam center; λ is the radar carrier frequency wavelength;
and
doppler phase history compensation items of the main antenna and the auxiliary antenna are respectively used for completing airborne SAR motion compensation while azimuth imaging is carried out; f. of 1 (x i ，y j ) And f 2 (x i ，y j ) Imaging signal models corresponding to the main antenna and the auxiliary antenna respectively;
and (4) removing the external terrain phase while imaging for the terrain phase item.
In some embodiments, said conjugate multiplying said digital orthophoto map to obtain a residual interference phase comprises: and carrying out conjugate multiplication on points corresponding to the two registered complex images of the main antenna and the auxiliary antenna to obtain a residual interference phase, wherein the residual interference phase is as follows:
wherein, delta phi is a residual interference phase;
the real slant distances of the main antenna and the auxiliary antenna at the crossing moment of the center of the azimuth beam are respectively.
In some embodiments, the interference calibration result parameters include a baseline length, a baseline tilt angle, and an initial phase offset, the ultra-sparse control point area network combined interference calibration is performed according to the interference phase map and DEM auxiliary conditions, and the interference calibration result parameters are calculated, where the interference calibration result parameters include: performing inversion on baseline parameters through interference phase gradients, combining a plurality of control point information and external DEM information in a multi-scene of the area network, and calculating the length of the baseline and the inclination angle of the baseline; and calculating initial phase offset based on a least square method by using the control point, the homologous point and the external DEM parameter.
In some embodiments, the back projection imaging model includes a primary antenna and a secondary antenna, the building of an elevation inversion model based on the back projection imaging model, inverting a ground elevation in the elevation inversion model using the interferometric calibration result parameters, the steps including: based on the absolute geospatial coordinate system, lifting the reference height upwards to construct an imaging reference surface; on the imaging reference surface, a target point is respectively reconstructed into a first imaging point and a second imaging point by a main antenna and an auxiliary antenna, so that interference image registration is realized and interference processing of image splicing is avoided; calculating an interference phase difference caused by the height difference according to the height difference between the real height of the target point and the reference plane to obtain the height difference; and compensating the initial DEM elevation by using the height difference, and inverting the ground elevation.
In some embodiments, the height difference-induced interference phase difference is:
wherein,
is an interference phase difference;
is the phase of the target point;
is the phase of the first imaging point;
is the initial phase offset; λ is the radar carrier frequency wavelength;
is a main dayThe slope distance from the line phase center to the target point;
the slant distance from the phase center of the auxiliary antenna to the target point;
the slant distance from the phase center of the main antenna to the first imaging point;
the slant distance from the phase center of the auxiliary antenna to the second imaging point; b is the length of a base line, specifically the distance between the main antenna and the auxiliary antenna; theta is a radar incident angle; alpha is a baseline inclination angle; Δ h is the height difference.
(III) advantageous effects
Compared with the prior art, the invention has the beneficial effects that:
in consideration of the fact that in recent years, SAR/InSAR imaging modes are more and more complex (beam bunching, circular trace, three-dimensional, MIMO, array and the like), imaging resolution is gradually developed from decimeter to centimeter, and the imaging resolution is more prone to being applied to light and small platforms such as unmanned planes, and higher requirements are put forward on airborne SAR motion compensation algorithms; the InSAR height measurement processing from images to DOM and DEM data products which can be applied by users is too specialized and complicated, and particularly, interference fringes are dense in an area with large topographic relief, so that difficulty is brought to phase unwrapping. In addition, the important factor for limiting the application of the InSAR in the surveying and mapping field is that the number of control points is large during large-area surveying and mapping, and the large workload is brought to the field of topographic survey and mapping due to the fact that the InSAR control points need to be manually distributed and synchronously acquired by radar data acquisition.
The invention has the idea that under the geographic information space coordinate system of a final data product, an imaging frame of backward projection point-by-point pulse-by-pulse coherent accumulation is adopted, so that the motion error caused by terrain can be compensated pixel by pixel, and the full aperture compensates the space-variant error quantity of each pixel position, and theoretically, the algorithm has no motion compensation error at all; meanwhile, with the assistance of an external DEM, the method can compensate mismatch caused by topographic relief, not only can remove flat ground, but also can remove most of topographic interference phases, so that the interference phases are extremely sparse and even not fuzzy, and the step of unwrapping the interference phases is avoided; in the stage of large-area interference measurement of the area network, the DEM auxiliary function can be played, and the number of control points is reduced. Compared with the traditional processing method, the method has the advantages of simpler processing process, higher processing precision, more friendly user application, more sparse dependent control points and the like.
The technical scheme of the invention has the remarkable advantages of strong mode adaptability, high precision of motion compensation and phase maintenance, simpler interference processing, self-registration of imaging completed simultaneously, DEM and SAR data matching avoided, multi-scene image splicing avoided, and relatively more user-friendly output products.
Drawings
Fig. 1 is a flowchart of an InSAR imaging interference integrated processing method according to an embodiment of the present invention.
Fig. 2 is an operation flowchart of an InSAR imaging interference integrated processing method according to an embodiment of the present invention.
Fig. 3 is a geometric schematic diagram of an InSAR back projection imaging model according to an embodiment of the present invention.
FIG. 4 is a flowchart illustrating operations of joint interferometric calibration of an ultra-sparse control point area network according to an embodiment of the present invention.
FIG. 5 is a geometric schematic of an elevation inversion model according to an embodiment of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to specific embodiments and the accompanying drawings.
In the present disclosure, the terms "include" and "comprise," as well as derivatives thereof, mean inclusion without limitation; the term "or" is inclusive, meaning and/or.
In this specification, the various embodiments described below which are used to describe the principles of the present disclosure are by way of illustration only and should not be construed in any way to limit the scope of the disclosure. The following description with reference to the accompanying drawings is provided to assist in a comprehensive understanding of exemplary embodiments of the present disclosure as defined by the claims and their equivalents. The following description includes various specific details to aid understanding, but such details are to be regarded as illustrative only. Accordingly, those of ordinary skill in the art will recognize that various changes and modifications of the embodiments described herein can be made without departing from the scope and spirit of the present disclosure. Moreover, descriptions of well-known functions and constructions are omitted for clarity and conciseness. Moreover, throughout the drawings, the same reference numerals are used for similar functions and operations.
Fig. 1 is a flowchart of an InSAR imaging interference integrated processing method according to an embodiment of the present invention. Fig. 2 is an operation flowchart of an InSAR imaging interference integrated processing method according to an embodiment of the present invention.
With reference to fig. 2, to explain fig. 1 in detail, the invention provides an InSAR imaging and interference integrated processing method, which includes steps S1 to S5.
S1, positioning the geographic space of the image in an absolute geographic space coordinate system, and constructing a back projection imaging model for integration of InSAR imaging and interferometric height measurement.
The absolute geospatial coordinate system is a WGS84 rectangular spatial coordinate system, and in the embodiment of the invention, the geographic space of the image is located in the absolute geospatial coordinate system.
It should be appreciated that in the area network imaging context, where conventional area network imaging methods require binning, the motion compensation of the imaging method includes local orientation and distance coordinate systems. The WGS84 spatial rectangular coordinate system is an international universal geocentric coordinate system, and in the embodiment of the present invention, the WGS84 spatial rectangular coordinate system is used as an imaging coordinate system. On the basis, each pixel of the image adopts gridding division of rasterization on the reference ellipsoid of the WGS84 space rectangular coordinate system, so that no shot exists. In fact, the method is limited by computer resources, and data scenes can be divided physically, but because each scene is completed based on the WGS84 space rectangular coordinate system, the method does not need multi-scene image splicing.
Fig. 3 is a geometric schematic diagram of an InSAR back projection imaging model according to an embodiment of the present invention.
The invention establishes an InSAR back projection imaging model, and a WGS84 space rectangular coordinate system is adopted as a coordinate system to position the geographic space of an image. Specifically, through geocoding, the external DEM elevation value is converted into a WGS84 space rectangular coordinate system and is used as an imaging coordinate system, and the external DEM elevation in the imaging coordinate system is obtained.
Due to the fact that elevation errors exist between the real elevation and the external DEM elevation under the imaging coordinate system, the real elevation and the external DEM elevation are different. As shown in FIG. 3, a target point P is arbitrarily selected from the rasterized geocoded image, and the actual elevation and the external DEM elevation of the target point are respectively h i And h d To illustrate.
It should be understood that the synthetic aperture time is the time a point object passes through a beam. In fig. 2, there are two antennas consisting of a main antenna and a sub-antenna, and the distance between the two antennas in the synthetic aperture radar system is the base length. In addition, the reference slant range of the interferometric synthetic aperture radar system refers to the distance from the main antenna to the reference point of the imaging space at the middle position of the synthetic aperture length in the interferometric synthetic aperture radar system.
As shown in FIG. 3, the true slant ranges of the primary antenna and the secondary antenna at the time of crossing the center of the azimuth beam based on the true elevation of the target point in the azimuth synthetic aperture time of the target point are respectively shown as
Correspondingly, in the azimuth synthetic aperture time of the target point, the calculated slant distances of the external DEM of the main antenna and the auxiliary antenna at the azimuth beam center crossing moment are respectively expressed as
The external DEM calculation slant distances of the main antenna and the auxiliary antenna set at any time in the azimuth synthetic aperture time of the target point are respectively expressed as
It will be appreciated that the azimuth beam center crossing time is just one particular time of the arbitrary times.
And S2, carrying out coherent accumulation imaging by using the back projection imaging model, outputting a digital orthophoto map, and carrying out conjugate multiplication on the digital orthophoto map to obtain a residual interference phase.
The basic idea of the back projection imaging algorithm is to find out a corresponding accumulation curve for coherent accumulation by calculating the two-way time delay from each pixel in an imaging area to an SAR antenna platform in an aperture length, thereby recovering the target function of each pixel. From the perspective of signal processing, the back projection imaging algorithm is a process of signal coherent accumulation, because the concept of phase is replaced by time delay, the back projection imaging algorithm is independent of frequency, thereby avoiding Fresnel approximation or other geometric approximation, and the back projection imaging algorithm is suitable for SAR imaging of any motion of a carrier.
In the embodiment of the invention, the original radar echo signal is set to be represented as s (r, x), wherein r and x respectively represent pixel coordinates of the distance direction and the azimuth direction of a target point. Performing range matching filtering p on the radar original echo signal * (-r) obtaining a range-compressed image of
s M (r，x)＝s(r，x)*p * (-r)＝∫ τ s(τ，x)p * (τ-r)dτ
Wherein s is M (r, x) is a distance compressed image; p is a radical of * (-r) represents a distance direction matching filter function, specifically, a conjugate complex number is taken for the distance direction r; s (r, x) p * (-r) represents the convolution of the radar original echo signal with the range-wise matched filter function; the third term of the equation represents the expansion formula of the convolution.
And (3) integrating (namely coherent accumulation) any pixel point in the rasterized geocoded image within the azimuth synthetic aperture time of the pixel point, so that the imaging of the pixel point along the azimuth direction can be realized.
And synthesizing the distance compressed image and the imaging of the pixel point along the azimuth direction, and finally obtaining the image which is a complex image. Without considering the local azimuth motion compensation error and the elevation error, the phase of the complex image only remains the phase of the backscattering coefficient, and does not contain the phase at the reference slope used for inverting the ground elevation information in the conventional interferometric SAR.
In InSAR imaging, the phase item of the crossing moment of the center of the azimuth beam is reserved, and the imaging signal models respectively corresponding to the main antenna and the auxiliary antenna are
Wherein,
respectively calculating the slant distances of the main antenna and the auxiliary antenna at any time according to the external DEM;
the main antenna and the auxiliary antenna which are respectively at the time of passing through the center of the azimuth beam are calculated to obtain slant distances according to the external DEM;
and
compressing the contribution of the image to the current pixel point for the distance of the main antenna and the auxiliary antenna at the crossing moment of the azimuth beam center; λ is the radar carrier frequency wavelength;
and
doppler phase history compensation items of the main antenna and the auxiliary antenna are respectively used, and high-precision airborne SAR motion compensation is completed while azimuth imaging is carried out; f. of 1 (x i ，y j ) And f 2 (x i ，y j ) Imaging signal models corresponding to the main antenna and the auxiliary antenna respectively;
and (4) removing the external terrain phase while imaging for the terrain phase item.
In the formula of the imaging signal model, coherent accumulation in synthetic aperture time can be realized after compensation of Doppler phase history compensation terms. The contribution of the distance compression image of the main antenna and the auxiliary antenna at the crossing moment of the azimuth beam center to the current pixel point can be accurately calculated through interpolation.
Based on f 2 (x i ，y j ) The imaging process of the secondary antenna eliminates the terrain interference phase. In addition, because the dual antennas are imaged under a uniform coordinate system, a uniform spatial baseline exists in the imaging process. That is, the imaging process compensates for both topographic relief and the amount of phase mismatch caused by the spatial baseline.
Through the steps, in an absolute geographic space coordinate system, residual motion errors caused by terrain fluctuation and orientation space variation are compensated by using an imaging overall frame of backward projection point-by-point pulse-by-pulse coherent accumulation, and imaging results are directly output to a digital orthophoto map in a geographic reference coordinate system.
Further, the point conjugate multiplication corresponding to the two complex images after the registration of the main antenna and the auxiliary antenna is carried out to obtain an interference phase, wherein the interference phase is
Wherein, delta phi is an interference phase;
the real slant distances of the main antenna and the auxiliary antenna at the crossing moment of the center of the azimuth beam are respectively.
The interference phase delta phi is the residual interference phase that removes the initial topography.
It should be noted that when the initial DEM error is smaller than the InSAR ambiguity height, the residual interference phase will not be wrapped, and the phase unwrapping processing step in the traditional DEM-assisted InSAR processing algorithm is avoided.
And S3, performing phase filtering on the residual interference phase, and outputting an interference phase diagram.
And (4) performing interference phase filtering on the residual interference phase by adopting a conventional filtering method, and outputting an interference phase image. The conventional filtering method may be, for example, median filtering or Goldstein frequency domain filtering, and the invention is not limited thereto.
And S4, performing ultra-sparse control point area network combined interference calibration according to the interference phase diagram and the DEM auxiliary condition, and calculating interference calibration result parameters.
The interference calibration result parameters comprise a base line length, a base line inclination angle and initial phase offset.
FIG. 4 is a flowchart illustrating operations of joint interferometric calibration of an ultra-sparse control point area network according to an embodiment of the present invention. As shown in FIG. 4, external DEM assisted sparse control point joint interferometric scaling includes sub-steps S41-S42.
S41, a baseline parameter is inverted through interference phase gradient, information of a plurality of control points and external DEM information are combined in multiple scenes of the area network, and the length of the baseline and the inclination angle of the baseline are calculated.
Referring to fig. 4, the baseline parameters are inverted by the interferometric phase gradients, that is, the phase change rate calculation of the differential interferometric phase domain is performed based on the interference data of the area network. Furthermore, a control point baseline estimation equation is established by combining a plurality of control point information in a multi-scene of the area network, and an external DEM baseline estimation equation is established by combining external DEM information in the multi-scene of the area network.
In the control point baseline estimation equation, the functional relation between InSAR interference phase and baseline parameters is
Wherein, delta r is InSAR interference phase; b is the base length; theta is a radar incident angle; alpha is a baseline inclination angle; λ is the radar carrier frequency wavelength; delta phi is a residual interference phase;
is the initial phase offset.
The InSAR interference phase is subjected to distance direction derivation to obtain
It can be seen that in the differential interference phase domain, the initial phase is biased
And the elimination is smooth.
d theta/dr represents the degree of the change of the incidence angle along with the radar slant range, is related to the elevation value of the target and is obtained by calculating the elevation value of the external DEM
Where dh/dr represents the degree of relief.
In an external DEM baseline estimation equation, an exogenous DEM error model h = h epsilon + hc + xi is established
Wherein h is c The constant error of the DEM is adopted, and h epsilon is an external terrain elevation value with an error; xi is random error, h is exogenous DEM error.
And substituting the formula of the external source DEM error model into the formula of d theta/dr to calculate dh/dr. Due to constant error h c The degree of relief is not affected and therefore the base length B and the base inclination angle alpha can be calculated.
It should be noted that theoretically, an equation can be established for the elevation value of each pixel point in the image, and the influence of a large random error in the external DEM error can be eliminated by using a massive constraint equation and a least square method, so that the high-precision estimation of the baseline parameters is realized.
And S42, calculating initial phase offset based on a least square method by using the control point, the homologous point and the external DEM parameter.
After the above step S41, the base line length and the base line inclination are known quantities, the base line parameters are eliminated in the interferometric scaling equation, and only the unknown number of the interferometric phase offset remains.
Specifically, the exogenous DEM error model is substituted into a calibration equation to obtain
Wherein, the second term of the equation is a scaling equation,
the above equations may be solved, for example, using a least squares or optimization method, which is not limiting in particular to the invention. In addition, in order to reduce the influence of random errors in elevation errors, a large amount of elevation information constraint equations can be used, so that initial phase offset parameters with high precision can be calculated.
In some embodiments, in order to further improve the calibration precision of the interference parameters and improve the robustness of the method, when the area network is subjected to large-area interference measurement, the homonymy point information of the overlapping area between adjacent images is extracted, and a homonymy point constraint equation and an external DEM constraint equation are combined to solve baseline parameters (including baseline length and baseline inclination angle) and initial phase offset parameters.
The external DEM-assisted ultra-sparse control point area network combined interference scaling equation can be written in the form of
The above equations may be solved, for example, using least squares or optimization, to obtain the baseline length B, baseline tilt angle α, and initial phase offset
And S5, establishing an elevation inversion model based on the back projection imaging model, and inverting the ground elevation in the elevation inversion model by using the interference calibration result parameters.
Although the backward projection imaging algorithm can realize integrated signal processing of imaging and interference height measurement and introduces external DEM assistance, the Doppler phase history compensation item is compensated for azimuth focusing in the imaging process by the backward projection imaging algorithm, so that the phase at the reference slant distance used for inverting the ground elevation information in the traditional interference SAR is not contained in the complex image obtained by imaging.
FIG. 5 is a geometric schematic of an elevation inversion model according to an embodiment of the invention. This step includes substeps S51-S54.
And S51, lifting the reference height to construct an imaging reference surface based on the absolute geospatial coordinate system.
In the embodiment of the invention, the absolute geographic space coordinate system is a WGS84 space rectangular coordinate system, and on the basis, the invention constructs an imaging reference surface.
Referring to fig. 5, the coordinate system is a reference coordinate system newly established based on the WGS84 space rectangular coordinate system, specifically, the z-axis of the WGS84 space rectangular coordinate system is raised from 0 point to a reference height p 'where p' is located, and the reference height is raised to the reference coordinate system, wherein the z-axis is the reference height direction, the y-axis points to the flight direction of the radar platform, and the x-axis is determined by the right-hand rule. The imaging reference plane is the xOy plane in fig. 5.
And S52, reconstructing a target point into two imaging points by the main antenna and the auxiliary antenna on the imaging reference surface respectively, and realizing interference image registration and avoiding an interference processing step of image splicing.
Referring to FIG. 5, in the reference coordinate system, the main antenna is the main antenna, the sub-antenna is the sub-antenna, A 1 And A 2 Respectively, the phase centers of the main and sub antennas.
Through step S51, the absolute elevation of the xOy imaging reference plane relative to the WGS84 spatial rectangular coordinate system is z p′ Thus the phase center A of the main antenna 1 Relative elevation with respect to a reference coordinate system is H-z p′ 。
The target point is P point, the height difference between the P point and the xOy imaging reference surface is delta h, and the height difference is also initial DElevation error of EM. Because the height difference delta h exists between the real height of the target point P and the reference plane, the target point P is respectively reconstructed as an imaging point Q by the main antenna and the auxiliary antenna on the imaging reference plane 1 ，Q 2 。
And S53, calculating an interference phase difference caused by the height difference according to the height difference between the actual height of the target point and the reference plane to obtain the height difference.
Referring to fig. 5, in the reference coordinate system, θ is the antenna phase center a 1 The side view angle when imaging the target point P, i.e. the angle of incidence of the radar beam. Antenna phase center A 1 And an antenna phase center A 2 The distance between the two is the base length B, the included angle between the base length B and the horizontal direction is alpha, and R is the antenna phase center A 1 The slope to the target point P.
The interference phase difference caused by the height difference Δ g can be converted into the target point P and the image point Q 1 The interference phase difference between them, in combination with the geometrical relationship of FIG. 5, is
Wherein,
is an interference phase difference;
is the phase of the target point P;
is an imaging point Q 1 The phase of (d);
is the initial phase offset; lambda is the radar carrier frequency wavelength;
is the antenna phase center A 1 To the target pointThe pitch of P, i.e., R in FIG. 4;
the slant distance from the antenna phase center A2 to a target point P;
is the antenna phase center A 1 To the imaging point Q 1 The skew distance of (2);
from the antenna phase center A2 to the imaging point Q 2 The skew distance of (a); b is the base length; theta is a radar incident angle; alpha is a baseline inclination angle;
is the initial phase offset; Δ h is the height difference.
The interference phase difference
The formula establishes a constraint relation between the residual interference phase and the height difference delta h.
Further, the above formula is converted to obtain a height difference Δ h of
And S54, compensating the initial DEM elevation by using the height difference, and inverting the ground elevation.
The initial DEM elevation value is compensated by using the height difference delta g obtained by the formula to obtain an accurate DEM elevation value, namely the real elevation of the ground, so that the ground elevation is inverted.
In summary, the invention provides an interferometric synthetic aperture radar (InSAR) imaging and interference integrated processing method, which is characterized in that under a geographic information space rectangular coordinate system, based on DEM assistance, a back projection imaging model is adopted to perform InSAR high-precision coherent imaging and interference integrated processing, interference processing steps such as registration and phase unwrapping are completed during imaging, and by using DEM assistance information, inSAR area network combined interference calibration is assisted, the number of control points is greatly reduced, and InSAR refined high-precision imaging and interference processing are realized under the DEM assistance condition.
The flowchart and block diagrams in the figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present disclosure. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams or flowchart illustration, and combinations of blocks in the block diagrams or flowchart illustration, can be implemented by special purpose hardware-based systems which perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.
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 (7)
1. An InSAR imaging interference integrated processing method is characterized by comprising the following steps:
positioning the geographic space of the image in an absolute geographic space coordinate system, and constructing a back projection imaging model for integrating InSAR imaging and interferometric height measurement;
performing coherent accumulation imaging by using the backward projection imaging model, outputting a digital orthophoto map, and performing conjugate multiplication on the digital orthophoto map to obtain a residual interference phase;
carrying out phase filtering on the residual interference phase and outputting an interference phase diagram;
performing ultra-sparse control point area network combined interference calibration according to the interference phase diagram and DEM auxiliary conditions, and calculating interference calibration result parameters;
and establishing an elevation inversion model based on the back projection imaging model, and inverting the ground elevation in the elevation inversion model by using the interference calibration result parameters.
2. The InSAR imaging and interference integrated processing method according to claim 1, wherein the absolute geospatial coordinate system is a WGS84 spatial rectangular coordinate system, and the locating the geospatial of the image comprises:
and converting the external DEM elevation value into the external DEM elevation under a WGS84 space rectangular coordinate system through geocoding.
3. The InSAR imaging and interference integrated processing method according to claim 1, wherein the back projection imaging model comprises a main antenna and a secondary antenna, and the performing coherent accumulation imaging by using the back projection imaging model comprises:
in InSAR imaging, the phase item of the crossing moment of the center of the azimuth beam is reserved, and the imaging signal models respectively corresponding to the main antenna and the auxiliary antenna are
Wherein,
respectively calculating the slant distances of the main antenna and the auxiliary antenna at any time according to the external DEM;
respectively calculating slant distances of the main antenna and the auxiliary antenna at the time of passing through the center of the azimuth beam according to an external DEM;
and
respectively compressing the contribution of the image to the current pixel point for the distance of the main antenna and the auxiliary antenna at the crossing moment of the azimuth beam center, wherein x represents the coordinate of the pixel of the azimuth direction of the target point; λ is the radar carrier frequency wavelength;
and
doppler phase history compensation items of the main antenna and the auxiliary antenna are respectively used for completing airborne SAR motion compensation while azimuth imaging is carried out; f. of 1 (x i ，y j ) And f 2 (x i ，y j ) Imaging signal models corresponding to the main antenna and the auxiliary antenna respectively;
and (4) removing the external terrain phase while imaging for the terrain phase item.
4. The InSAR imaging and interference integrated processing method according to claim 3, wherein the obtaining of the residual interference phase by conjugate multiplication of the digital orthophoto map comprises:
and carrying out conjugate multiplication on points corresponding to the two registered complex images of the main antenna and the auxiliary antenna to obtain a residual interference phase, wherein the residual interference phase is as follows:
wherein, deltaPhi is a residual interference phase;
the true slant distances of the main antenna and the auxiliary antenna at the crossing moment of the azimuth beam center are respectively.
5. The InSAR imaging interference integrated processing method according to claim 1, wherein the interferometric calibration result parameters include baseline length, baseline dip and initial phase offset, the ultra-sparse control point area network combined interferometric calibration is performed according to the interferometric phase diagram and DEM auxiliary conditions, and interferometric calibration result parameters are calculated, the steps include:
performing inversion on baseline parameters through interference phase gradients, combining a plurality of control point information and external DEM information in a multi-scene of the area network, and calculating the length of the baseline and the inclination angle of the baseline;
and calculating initial phase offset based on a least square method by using the control point, the homologous point and the external DEM parameter.
6. The InSAR imaging and interference integrated processing method according to claim 1, wherein the back projection imaging model comprises a main antenna and an auxiliary antenna, the elevation inversion model is established based on the back projection imaging model, and the ground elevation is inverted in the elevation inversion model by using the interference calibration result parameters, and the method comprises the following steps:
based on the absolute geospatial coordinate system, lifting the reference height upwards to construct an imaging reference surface;
on the imaging reference surface, a target point is respectively reconstructed into a first imaging point and a second imaging point by a main antenna and an auxiliary antenna, so that interference image registration is realized and interference processing of image splicing is avoided;
calculating an interference phase difference caused by the height difference according to the height difference between the real height of the target point and the reference plane to obtain the height difference;
and compensating the initial DEM elevation by using the height difference, and inverting the ground elevation.
7. The InSAR imaging and interference integrated processing method according to claim 6, wherein the interference phase difference caused by the height difference is as follows:
wherein,
is an interference phase difference;
is the phase of the target point;
is the phase of the first imaging point;
is the initial phase offset; λ is the radar carrier frequency wavelength;
the slant distance from the phase center of the main antenna to the target point;
the slant distance from the phase center of the auxiliary antenna to the target point;
the slant distance from the phase center of the main antenna to the first imaging point;
the slant distance from the phase center of the auxiliary antenna to the second imaging point; b is the length of a base line, specifically the distance between the main antenna and the auxiliary antenna; theta is a radar incident angle; alpha is a baseline inclination angle; Δ h is the height difference.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110039286.4A CN112882030B (en) | 2021-01-12 | 2021-01-12 | InSAR imaging interference integrated processing method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110039286.4A CN112882030B (en) | 2021-01-12 | 2021-01-12 | InSAR imaging interference integrated processing method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112882030A CN112882030A (en) | 2021-06-01 |
| CN112882030B true CN112882030B (en) | 2023-01-03 |
Family
ID=76044922
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110039286.4A Active CN112882030B (en) | 2021-01-12 | 2021-01-12 | InSAR imaging interference integrated processing method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112882030B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113589282B (en) * | 2021-07-12 | 2023-08-08 | 中国科学院国家空间科学中心 | Method for removing ground effect by satellite-borne interference imaging altimeter based on image domain transformation |
| CN113640797B (en) * | 2021-08-09 | 2022-04-12 | 北京航空航天大学 | Front squint height measurement method for reference stripe mode InSAR |
| CN114609633B (en) * | 2022-03-17 | 2023-09-01 | 电子科技大学 | SAR height measurement method by circumferential beam focusing mode interference |
| CN118011344B (en) * | 2024-04-02 | 2024-06-11 | 中国科学院空天信息创新研究院 | SAR interference calibration method assisted by external DEM |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5463397A (en) * | 1993-10-25 | 1995-10-31 | Hughes Aircraft Company | Hyper-precision SAR interferometry using a dual-antenna multi-pass SAR system |
| CN103018741A (en) * | 2012-12-11 | 2013-04-03 | 电子科技大学 | Interferometric synthetic aperture radar (InSAR) imaging and flat ground removing integral method based on back projection |
| CN103033813A (en) * | 2012-12-31 | 2013-04-10 | 长沙理工大学 | Initial phase offset real-time estimation method based on external coarse precision DEM (Dynamic Effect Model) |
| CN103630898A (en) * | 2013-03-27 | 2014-03-12 | 中国科学院电子学研究所 | Method for estimating multi-baseline interferometry SAR phase bias |
| CN105929399A (en) * | 2016-04-25 | 2016-09-07 | 电子科技大学 | Interference SAR data imaging and elevation estimation method |
-
2021
- 2021-01-12 CN CN202110039286.4A patent/CN112882030B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5463397A (en) * | 1993-10-25 | 1995-10-31 | Hughes Aircraft Company | Hyper-precision SAR interferometry using a dual-antenna multi-pass SAR system |
| CN103018741A (en) * | 2012-12-11 | 2013-04-03 | 电子科技大学 | Interferometric synthetic aperture radar (InSAR) imaging and flat ground removing integral method based on back projection |
| CN103033813A (en) * | 2012-12-31 | 2013-04-10 | 长沙理工大学 | Initial phase offset real-time estimation method based on external coarse precision DEM (Dynamic Effect Model) |
| CN103630898A (en) * | 2013-03-27 | 2014-03-12 | 中国科学院电子学研究所 | Method for estimating multi-baseline interferometry SAR phase bias |
| CN105929399A (en) * | 2016-04-25 | 2016-09-07 | 电子科技大学 | Interference SAR data imaging and elevation estimation method |
Non-Patent Citations (5)
| Title |
|---|
| DEM辅助的山区InSAR相位解缠;刘辉等;《测绘科学技术学报》;20171231;第34卷(第02期);全文 * |
| DEM辅助的零中频多基线InSAR干涉处理;刘辉等;《河南师范大学学报(自然科学版)》;20180710;第46卷(第05期);全文 * |
| InSAR DEM Reconstruction Based on Backprojection Algorithm in Two Converse Flights;Xiaoning HU等;《2019 6th Asia-Pacific Conference on Synthetic Aperture Radar (APSAR)》;20191129;全文 * |
| 一种机载双天线InSAR干涉参数定标新方法;靳国旺等;《测绘学报》;20100215;第39卷(第01期);全文 * |
| 基于曲面投影的毫米波InSAR数据成像方法;韦顺军等;《雷达学报》;20150215;第04卷(第01期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112882030A (en) | 2021-06-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112882030B (en) | InSAR imaging interference integrated processing method | |
| CN109633648B (en) | Multi-baseline phase estimation device and method based on likelihood estimation | |
| CN107102333B (en) | Satellite-borne InSAR long and short baseline fusion unwrapping method | |
| CN106526593B (en) | Sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR | |
| CN109239710B (en) | Method and device for acquiring radar elevation information and computer-readable storage medium | |
| CN107861920A (en) | cloud data registration method | |
| CN112711021B (en) | Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method | |
| JP2590689B2 (en) | Interferometric synthetic aperture radar system and terrain change observation method | |
| Wang et al. | Large-scale orthorectification of GF-3 SAR images without ground control points for China’s land area | |
| Crosetto et al. | Radargrammetry and SAR interferometry for DEM generation: validation and data fusion. | |
| Dong et al. | Radargrammetric DSM generation in mountainous areas through adaptive-window least squares matching constrained by enhanced epipolar geometry | |
| CN110703252B (en) | Digital elevation model correction method for interferometric synthetic aperture radar shadow area | |
| CN109946682B (en) | GF3 data baseline estimation method based on ICESat/GLAS | |
| KR100441590B1 (en) | Method of generating DEM for Topography Measurement using InSAR | |
| CN115712095A (en) | SAR satellite three-dimensional positioning error correction method and system based on single angular reflection | |
| Schmitt et al. | Towards airborne single pass decimeter resolution SAR interferometry over urban areas | |
| CN111696207B (en) | Multi-baseline DEM fusion method based on guided filtering | |
| CN111738135A (en) | SAR image feature extraction method considering slant range and geographic coordinate mapping deviation | |
| Li et al. | Accuracy of airborne IFSAR mapping | |
| Nascetti et al. | Radargrammetric digital surface models generation from high resolution satellite SAR imagery: Methodology and case studies | |
| CN117761695B (en) | Multi-angle SAR three-dimensional imaging method based on self-adaptive partition SIFT | |
| Benyi et al. | The principles of positioning with space-borne SAR images | |
| Biswas et al. | Investigation of ground deformation with PSInSAR approach in an unstable urban area Naples, Italy using X-band SAR images | |
| Sosnovsky et al. | Modification of inversed vortex phase field unwrapping algorithm for the InSAR height measurements | |
| CN117741698A (en) | Multi-scene information joint baseline calibration method and system |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |