CN110031841B - ECMWF-based InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system - Google Patents
ECMWF-based InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system Download PDFInfo
- Publication number
- CN110031841B CN110031841B CN201910257685.0A CN201910257685A CN110031841B CN 110031841 B CN110031841 B CN 110031841B CN 201910257685 A CN201910257685 A CN 201910257685A CN 110031841 B CN110031841 B CN 110031841B
- Authority
- CN
- China
- Prior art keywords
- delay
- zenith
- atmospheric
- research area
- value
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 52
- 238000012937 correction Methods 0.000 title claims abstract description 38
- 238000011160 research Methods 0.000 claims abstract description 100
- 238000004364 calculation method Methods 0.000 claims abstract description 51
- 239000005436 troposphere Substances 0.000 claims abstract description 23
- 238000001035 drying Methods 0.000 claims description 16
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 14
- 238000012821 model calculation Methods 0.000 claims description 12
- 239000007788 liquid Substances 0.000 claims description 6
- 239000000126 substance Substances 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 3
- 239000000835 fiber Substances 0.000 abstract description 2
- 230000000694 effects Effects 0.000 description 11
- 238000005516 engineering process Methods 0.000 description 5
- 238000003384 imaging method Methods 0.000 description 4
- 238000013517 stratification Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 230000008520 organization Effects 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 230000001427 coherent effect Effects 0.000 description 1
- 230000014509 gene expression Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000012950 reanalysis Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- CCEKAJIANROZEO-UHFFFAOYSA-N sulfluramid Chemical group CCNS(=O)(=O)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)F CCEKAJIANROZEO-UHFFFAOYSA-N 0.000 description 1
Images
Landscapes
- Radar Systems Or Details Thereof (AREA)
- Image Processing (AREA)
Abstract
The invention discloses an InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system based on ECMWF (equal cost multi-aperture fiber). The method comprises the following steps: calculating the troposphere zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area according to the ECMWF weather forecast data and the DEM data of the research area; obtaining the total zenith delay of the research area corresponding to the coarse scale of each SAR image according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area; constructing a layered interpolation model according to the total zenith delay of each SAR image; obtaining an atmospheric delay phase value of each point on the interference SAR image; and removing the atmospheric delay phase value from the interference SAR unwrapping image to obtain an atmospheric delay corrected interference image. The invention improves the resolution ratio of atmospheric delay correction and has accurate calculation result.
Description
Technical Field
The invention relates to the technical field of InSAR, in particular to an InSAR atmospheric delay correction method and system based on ECMWF.
Background
InSAR (Interferometric Synthetic Aperture Radar) is a Radar technology applied to mapping and remote sensing. The method is a technology for acquiring the elevation information of the earth surface by utilizing a synthetic aperture radar to carry out coherent processing on two complex value images (images with both amplitude and phase) observed in the same area. InSAR technology has been widely studied and applied in many fields, such as topographic survey, earthquake, landslide, volcano, large area surface subsidence, marine mapping, disaster monitoring, etc.
In application, the InSAR technology returns a signal through electromagnetic waves, and the electromagnetic waves are influenced by tropospheric atmosphere, particularly water vapor scattering, in the atmospheric propagation process, so that the propagation path of the electromagnetic waves is bent, and the radar echo signal is subjected to phase delay, namely atmospheric delay. For InSAR, the imaging time and the viewing angle of radar are different, the atmospheric conditions (relative humidity, temperature, air pressure and the like) are changed in time and space, the atmospheric conditions experienced by electromagnetic wave propagation are different in the imaging of the radar and the atmospheric phase delay can affect the interference and differential interference phases, and the accuracy of the InSAR technology is affected.
At present, the atmospheric delay correction method in the InSAR includes a method based on MERIS (a medium resolution imaging spectrometer) and MODIS (a medium resolution imaging spectrometer) water vapor products, a method based on GPS data, a method based on a medium-scale atmospheric prediction model wrf (weather Research and Forecast model), and the like. The method is limited by data time resolution, the receiving time of the data is fixed, and the data time and the time of an actual interference diagram have deviation to influence a calculation result; the atmospheric delay correction by using GPS data has higher accuracy, but the data density of the GPS station in China is lower, so the resolution of the obtained atmospheric correction image is also lower. While the data resolution of the WRF model is also relatively low.
Therefore, it is an urgent technical problem to be solved in the art to provide a method and a system for InSAR atmosphere delay correction based on ECMWF.
Disclosure of Invention
In view of this, the invention provides an method and a system for InSAR atmospheric delay correction based on ECMWF, which improve the atmospheric delay correction resolution and have accurate calculation results.
In order to solve the technical problem, the invention provides an InSAR atmospheric delay correction method based on ECMWF, which is characterized by comprising the following steps:
calculating the troposphere zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area according to the ECMWF weather forecast data and the DEM data of the research area;
obtaining the total zenith delay of the research area corresponding to the coarse scale of each SAR image according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area;
constructing a layered interpolation model according to the total zenith delay of each SAR image, wherein the layered interpolation model is as follows: Δ L ═ S (h) + T (x) + w, where Δ L is the total zenith delay difference of two views of SAR images, S is the stratified delay component of atmospheric delay, T is the turbulent delay component of atmospheric delay, w is the residual delay value, x is a certain point on the SAR image, and h is the elevation value at the corresponding point x;
calculating the layering delay value S of each point on the interference SAR image in the research area by adopting an interpolation mode, wherein the method comprises the following steps:
iterative function s (h) qexp { -b (h-h)min)/(hmax-hmin) Finding the optimal parameter values of q and b, wherein q and b are index parameters of the area to be interpolated, h is the elevation value of the point to be interpolated, and h is the elevation value of the point to be interpolatedmaxIs the maximum elevation value, h, of the SAR image areaminThe minimum elevation value of the SAR image area is obtained;
obtaining the layering delay value S of each point in the research area by using an interpolation method according to the optimal parameter values of q and b;
separating the layered delay values S of all points in the research area according to a layered interpolation model, and obtaining turbulence delay values T and residual delay values w of all points in the research area by adopting an inverse distance weight interpolation method;
obtaining an atmospheric delay phase value of each point on the interference SAR image according to the layering delay value S of each point in the research area and the turbulence delay value T of each point;
and removing the atmospheric delay phase value from the interference SAR unwrapping image to obtain an atmospheric delay corrected interference image.
Optionally, calculating troposphere zenith atmosphere dry delay and zenith atmosphere wet delay of the research area according to ECMWF weather forecast data and DEM data of the research area, including:
using a formula
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
using a formula
Calculating zenith atmosphere wet retardation, wherein ZWD is zenith atmosphere wet retardation, k' and k are refractive constants, ρwIs liquid water density, TMIs the weighted average temperature, R, of the tropospherevIs a specific atmospheric constant, PWV is the reducible moisture content.
Optionally, calculating troposphere zenith atmosphere dry delay and zenith atmosphere wet delay of the research area according to ECMWF weather forecast data and DEM data of the research area, including:
using a formula
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
calculating the zenith atmospheric humidity delay by adopting a formula ZWD pi PWV, wherein ZWD is the zenith atmospheric humidity delay, PWV is the reducible water vapor content, and pi is a scale factor.
Optionally, pi is more than or equal to 6.0 and less than or equal to 6.5.
Optionally, obtaining the total zenith delay of each SAR image of the research area according to the zenith atmospheric dry delay and the zenith atmospheric wet delay of the research area, including: ZTD is ZHD + ZWD, ZTD being the total zenith delay.
The invention also provides an system for InSAR (interferometric synthetic aperture radar) atmospheric delay correction based on ECMWF (equal cost multi-aperture fiber), which is characterized by comprising the following steps: the system comprises a data acquisition module, a total zenith delay calculation module, an interpolation model calculation module, an atmospheric delay phase value generation module and an atmospheric delay correction module; wherein the content of the first and second substances,
the data acquisition module is connected with the total zenith delay calculation module and is used for acquiring ECMWF weather forecast data and DEM data of a research area and sending the ECMWF weather forecast data and the DEM data of the research area to the total zenith delay calculation module;
the total zenith delay calculation module is connected with the difference model calculation module and used for calculating the dry delay of the zenith atmosphere and the wet delay of the zenith atmosphere of the troposphere of the research area according to the ECMWF weather forecast data and the DEM data of the research area, obtaining the total zenith delay of each SAR image in the research area corresponding to the coarse scale of each SAR image according to the dry delay of the zenith atmosphere and the wet delay of the zenith atmosphere of the research area, and sending the total zenith delay of each SAR image in the research area to the interpolation model calculation module;
the interpolation model calculation module is connected with the atmospheric delay phase value generation module and used for constructing a layered interpolation model according to the total zenith delay of each SAR image and calculating the layered delay value S of each point in the research area and the turbulence delay value T of each point according to the layered interpolation model, and the method comprises the following steps: a model construction unit, a hierarchical delay value S calculation unit and a turbulence delay value T calculation unit, wherein,
the model building unit is used for building a layered interpolation model, and the layered interpolation model is as follows: Δ L ═ S (h) + T (x) + w, Δ L is the total zenith delay difference of two views of SAR images, S is the stratified delay portion of atmospheric delay, T is the turbulent delay portion of atmospheric delay, w is the residual delay value, x is a certain point on the SAR image, and h is the elevation value at the corresponding point x;
the hierarchical delay value S calculation unit is used for calculating the hierarchical delay value S of each point on the interference SAR image in the research area in an interpolation mode, and firstly, according to an iteration function S (h), qexp, b (h-h)min)/(hmax-hmin) Finding the optimal parameter values of q and b, wherein q and b are index parameters of the area to be interpolated, h is the elevation value of the point to be interpolated, and h is the elevation value of the point to be interpolatedmaxIs the maximum elevation value, h, of the SAR image areaminThe minimum elevation value of the SAR image area is obtained; then obtaining the layered delay values S of all points in the research area by adopting an interpolation method through the optimal parameter values of q and b, and sending the layered delay values S of all points in the research area to an atmospheric delay phase value generation module;
the turbulence delay value T calculation unit is used for separating the layered delay value S of each point in the research area according to the layered interpolation model, obtaining the turbulence delay value T and the residual delay value w of each point in the research area by adopting an inverse distance weight interpolation method, and sending the turbulence delay value T of each point in the research area to the atmosphere delay phase value generation module;
the atmosphere delay phase value generation module is connected with the atmosphere delay correction module and used for obtaining an atmosphere delay phase value of each point on the interference SAR image according to the layering delay value S of each point in the research area and the turbulence delay value T of each point and sending the atmosphere delay phase value to the atmosphere delay correction module;
and the atmospheric delay correction module is used for removing the atmospheric delay phase value from the interference SAR unwrapping image to obtain an interference image after atmospheric delay correction.
Optionally, the total zenith delay calculating module includes an atmosphere dry delay calculating unit and an atmosphere wet delay calculating unit, wherein,
the atmosphere dry delay calculation unit is used for adopting a formula according to ECMWF weather forecast data and DEM data of a research area
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
the atmosphere moisture delay calculation unit is used for adopting a formula according to ECMWF weather forecast data and DEM data of a research area
Calculating zenith atmosphere wet retardation, wherein ZWD is zenith atmosphere wet retardation, k' and k are refractive constants, ρwIs liquid water density, TMIs the weighted average temperature, R, of the tropospherevIs a specific atmospheric constant, PWV is the reducible moisture content.
Optionally, the total zenith delay calculating module includes an atmosphere dry delay calculating unit and an atmosphere wet delay calculating unit, wherein,
the atmosphere dry delay calculation unit is used for adopting a formula according to ECMWF weather forecast data and DEM data of a research area
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
and the atmosphere moisture delay calculation unit is used for calculating the zenith atmosphere moisture delay by adopting a formula ZWD & Π & PWV according to ECMWF weather forecast data and DEM data of a research area, wherein ZWD is the zenith atmosphere moisture delay, PWV is the reducible water vapor content, and Π is a scale factor.
Optionally, the total zenith delay calculating module further includes a total delay calculating module; and the total delay calculation module is used for obtaining the total zenith delay of each SAR image of the research area according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area, wherein ZTD is ZHD + ZWD, and ZTD is the total zenith delay.
Compared with the prior art, the method and the system for InSAR atmosphere delay correction based on ECMWF provided by the invention at least realize the following beneficial effects:
according to the method and the system for correcting the atmospheric delay of the troposphere, InSAR atmospheric delay calculation is carried out by utilizing ECMWF data, the used ECMWF data has the forecast data time interval of 3 hours and is updated once every 6 hours, the timeliness is strong, the spatial resolution can reach 0.125 degrees at most, and the method and the system are meteorological forecast data with the highest resolution at present. The invention utilizes real-time effective meteorological parameters to calculate accurate and high-precision atmospheric delay. In addition, the invention simultaneously considers the distribution rules of the atmospheric turbulence effect and the layering effect, and utilizes the layering interpolation model to carry out interpolation calculation on the delay on the spatial height distribution, thereby improving the resolution ratio and obtaining the calculation result of the atmospheric delay with higher accuracy.
Other features of the present invention and advantages thereof will become apparent from the following detailed description of exemplary embodiments thereof, which proceeds with reference to the accompanying drawings.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description, serve to explain the principles of the invention.
Fig. 1 is a flow chart of a method for correcting atmospheric delay of an InSAR troposphere based on ECMWF provided by the invention;
fig. 2 is a system block diagram of the ECMWF-based InSAR tropospheric atmosphere delay correction according to an embodiment of the present invention.
Detailed Description
Various exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be noted that: the relative arrangement of the components and steps, the numerical expressions and numerical values set forth in these embodiments do not limit the scope of the present invention unless specifically stated otherwise.
The following description of at least one exemplary embodiment is merely illustrative in nature and is in no way intended to limit the invention, its application, or uses.
Techniques, methods, and apparatus known to those of ordinary skill in the relevant art may not be discussed in detail, but are intended to be part of the specification where appropriate.
In all examples shown and discussed herein, any particular value should be construed as merely illustrative, and not limiting. Thus, other examples of the exemplary embodiments may have different values.
It should be noted that: like reference numbers and letters refer to like items in the following figures, and thus, once an item is defined in one figure, further discussion thereof is not required in subsequent figures.
The European middle-term Weather forecasting center (ECMWF) is an international organization comprising 24 European Union member countries, and is a unique international Weather forecasting research and business organization in the world today. The third generation ERA-Interim of ECMWF reanalysis data is currently updated in real time, using four-dimensional variational analysis, the model contains 60 layers in the vertical direction, the highest layer reaches 0.1hPa, providing assimilation analysis data four times a day, including various atmospheric parameters such as temperature, pressure, humidity, etc. The inventor considers that the data based on the ECMWF corrects the atmospheric delay of the troposphere by InSAR by using the higher spatial resolution and the time resolution thereof, and improves the resolution of atmospheric delay correction and the accuracy of a calculation result.
Fig. 1 is a flow chart of a method for correcting atmospheric delay of an InSAR troposphere based on ECMWF provided by the invention. As shown in fig. 1, the method provided by the present invention comprises:
step S101: calculating the troposphere zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area according to the ECMWF weather forecast data and the DEM data of the research area; the ECMWF weather forecast data comprises parameters such as air pressure, temperature, humidity and longitude and latitude. Optionally, the method of the present invention selects weather forecast data with a resolution of 0.125 ° from ECMWF weather forecast data when applied.
Zenith atmospheric dry Delay (Zenith static Delay, ZHD) can be calculated from data such as ground point latitude, elevation (height) and surface air pressure, optionally, formula can be adopted
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; and H represents the elevation of a ground point, and the unit is km.
The Zenith atmospheric moisture Delay (ZWD) can be obtained by calculating the reducible water vapor content PWV, wherein the reducible water vapor content PWV is the physical quantity of the water vapor condition in the atmosphere and can be directly obtained by the existing space-based monitoring system.
Alternatively, formulas may be employed
Calculating zenith atmospheric moisture delay, wherein ZWD is zenith atmospheric moisture delay, k is a refractive constant, is liquid water density, is a weighted average temperature of the troposphere, is a specific atmospheric constant, and PWV is reducible moisture content.
Optionally, the zenith atmospheric moisture delay may be calculated using a formula ZWD ═ Π · PWV, where ZWD is the zenith atmospheric moisture delay, PWV is the reducible water vapor content, and Π is the scaling factor. Furthermore, pi is more than or equal to 6.0 and less than or equal to 6.5. For a coarse transition between ZWD and PWV, pi may be chosen to have an average transition factor of 6.2.
Step S102: obtaining the total zenith delay corresponding to the coarse size of each SAR image in the research area according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area; tropospheric atmospheric delay is mainly composed of zenith dry delay and zenith wet delay. In this method, ZTD is defined as total zenith delay, where ZTD is ZHD + ZWD.
Step S103: constructing a hierarchical interpolation model (SDI model) according to the total zenith delay of each SAR image, wherein the hierarchical interpolation model comprises the following steps: Δ L ═ S (h) + T (x) + w, where Δ L is the difference in total zenith delay of two views of SAR images, i.e. the atmospheric delay value on the interferometric SAR image, S is the stratified delay part of the atmospheric delay, T is the turbulent delay part of the atmospheric delay, w is the residual delay value, x is a certain point on the SAR image, and h is the elevation value at the corresponding point x.
The inventor considers that the resolution of 0.125-degree weather forecast data in the ECMWF is higher and is closer to that of the SAR image, but the resolution of the SAR image is generally several meters or dozens of meters, and the resolution calculated by the ECMWF data alone is far insufficient, so that a layered interpolation model is further constructed. Considering that the turbulence fluctuation of temperature and humidity is caused by weather conditions and strong steam change, the vertical wind shear force or the strong convection effect generated by a thin turbulence layer in the cloud is a random variable quantity, so that the refractive index of the atmosphere shows spatial heterogeneity, thereby causing local phase gradient and causing turbulence effect. In addition, the earth's atmospheric profile exhibits varying degrees of stratification, i.e., a stratification effect, over a wide range of vertical scales, which stratification may be severely disturbed by tropospheric turbulence when local scales (e.g., interferograms) are considered. Therefore, the inventor takes the layering effect and the turbulence effect into consideration when the delay result is calculated, and constructs the layering interpolation model to improve the accuracy of atmospheric delay correction.
Step S104: calculating the layering delay value S of each point on the interference SAR image in the research area by adopting an interpolation mode, wherein the method comprises the following steps:
iterative function s (h) qexp { -b (h-h)min)/(hmax-hmin) Finding the optimal parameter values of q and b, wherein q and b are index parameters of the area to be interpolated, h is the elevation value of the point to be interpolated, and h is the elevation value of the point to be interpolatedmaxIs the maximum elevation value, h, of the SAR image areaminThe minimum elevation value of the SAR image area is obtained;
step S105: obtaining the layering delay value S of each point in the research area by using an interpolation method according to the optimal parameter values of q and b; in this step, the optimal parameters q and b are found through an iterative function, so that the layered delay components at different heights in the research area can be separated, that is, the layered delay values S of each point in the research area are obtained by adopting an interpolation method.
Step S106: separating the layered delay values S of all points in the research area according to a layered interpolation model, and obtaining turbulence delay values T and residual delay values w of all points in the research area by adopting an inverse distance weight interpolation method; when the turbulence delay value T is calculated, the inventor considers that a spline or bilinear interpolation method is suitable for the condition of gradual change in a large area range, but is not suitable for large change in a short horizontal distance, so that the two methods are not suitable for extreme weather conditions, and therefore the turbulence delay value T of each point in a research area is calculated by adopting an IDW (inverse distance weight) interpolation mode so as to reasonably calculate the turbulence delay value T.
Step S107: obtaining an atmospheric delay phase value of each point on the interference SAR image according to the layering delay value S of each point in the research area and the turbulence delay value T of each point;
step S108: and removing the atmospheric delay phase value from the interference SAR unwrapping image to obtain an atmospheric delay corrected interference image.
The method for correcting the atmospheric delay of the troposphere utilizes ECMWF data to calculate InSAR atmospheric delay, the used ECMWF data has the forecast data time interval of 3 hours and is updated once every 6 hours, the timeliness is strong, the spatial resolution can reach 0.125 degrees at most, and the method is the meteorological forecast data with the highest resolution at present. The invention utilizes real-time effective meteorological parameters to calculate accurate and high-precision atmospheric delay. In addition, the invention simultaneously considers the distribution rules of the atmospheric turbulence effect and the layering effect, and utilizes the layering interpolation model to carry out interpolation calculation on the delay on the spatial height distribution, thereby improving the resolution ratio and obtaining the calculation result of the atmospheric delay with higher accuracy.
The invention also provides a system for correcting the InSAR troposphere atmosphere delay based on the ECMWF, and FIG. 2 is a system block diagram of the InSAR troposphere atmosphere delay correction based on the ECMWF provided by the embodiment of the invention. As shown in fig. 2, the present invention provides a system 200 comprising: the system comprises a data acquisition module 201, a total zenith delay calculation module 202, an interpolation model calculation module 203, an atmospheric delay phase value generation module 204 and an atmospheric delay correction module 205; wherein the content of the first and second substances,
the data acquisition module 201 is connected with the total zenith delay calculation module 202, and is configured to acquire ECMWF weather forecast data and DEM data of a research area, and send the ECMWF weather forecast data and the DEM data of the research area to the total zenith delay calculation module 202;
the total zenith delay calculation module 202 is connected with the interpolation model calculation module 203, and is used for calculating the dry delay of the troposphere zenith atmosphere and the wet delay of the zenith atmosphere in the research area according to the ECMWF weather forecast data and the DEM data of the research area, obtaining the total zenith delay of each SAR image in the research area corresponding to the coarse scale of each SAR image according to the dry delay of the zenith atmosphere and the wet delay of the zenith atmosphere in the research area, and sending the total zenith delay of each SAR image in the research area to the interpolation model calculation module 203;
optionally, the total zenith delay calculating module 202 includes an atmospheric dry delay calculating unit, an atmospheric wet delay calculating unit and a total delay calculating module; wherein the content of the first and second substances,
the atmosphere dry delay calculation unit is used for adopting a formula according to ECMWF weather forecast data and DEM data of a research area
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
and the atmosphere moisture delay calculation unit is used for calculating DEM data of the research area according to the ECMWF weather forecast data.
In one embodiment, a formula may be employed
Calculating zenith atmosphere wet retardation, wherein ZWD is zenith atmosphere wet retardation, k' and k are refractive constants, ρwIs liquid water density, TMIs the weighted average temperature, R, of the tropospherevIs a specific atmospheric constant, PWV is the reducible moisture content.
In another embodiment, the zenith atmospheric moisture delay may be calculated using the formula ZWD ═ Π · PWV, where ZWD is zenith atmospheric moisture delay, PWV is reducible moisture content, and Π is a scaling factor.
And the total delay calculation module is used for obtaining the total zenith delay of each SAR image of the research area according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area, wherein ZTD is ZHD + ZWD, and ZTD is the total zenith delay.
The interpolation model calculation module 203 is connected to the atmospheric delay phase value generation module 204, and configured to construct a hierarchical interpolation model according to the total zenith delay of each SAR image, and calculate a hierarchical delay value S of each point in the research area and a turbulence delay value T of each point according to the hierarchical interpolation model, including: a model construction unit 2031, a stratification delay value S calculation unit 2032, and a turbulence delay value T calculation unit 2033, wherein,
the model building unit 2031 is configured to build a hierarchical interpolation model, where the hierarchical interpolation model is: Δ L ═ S (h) + T (x) + w, Δ L is the total zenith delay difference of two views of SAR images, S is the stratified delay portion of atmospheric delay, T is the turbulent delay portion of atmospheric delay, w is the residual delay value, x is a certain point on the SAR image, and h is the elevation value at the corresponding point x;
the hierarchical delay value S calculating unit 2032 is configured to calculate the hierarchical delay value S of each point on the interferometric SAR image in the study region by interpolation, and first calculate the hierarchical delay value S according to an iterative function S (h) ═ qexp { -b (h-h)min)/(hmax-hmin) Finding the optimal parameter values of q and b, wherein q and b are index parameters of the area to be interpolated, h is the elevation value of the point to be interpolated, and h is the elevation value of the point to be interpolatedmaxIs the maximum elevation value, h, of the SAR image areaminThe minimum elevation value of the SAR image area is obtained; then, obtaining the layering delay value S of each point in the research area by adopting an interpolation method through the optimal parameter values of q and b, and sending the layering delay value S of each point in the research area to the atmosphere delay phase value generation module 204;
the turbulence delay value T calculation unit 2033 is configured to separate the hierarchical delay value S of each point in the research area according to the hierarchical interpolation model, obtain the turbulence delay value T and the residual delay value w of each point in the research area by using an inverse distance weighted interpolation method, and send the turbulence delay value T of each point in the research area to the atmospheric delay phase value generation module 204;
the atmospheric delay phase value generation module 204 is connected to the atmospheric delay correction module 205, and is configured to obtain an atmospheric delay phase value of each point on the interference SAR image according to the layered delay value S of each point in the research area and the turbulence delay value T of each point, and send the atmospheric delay phase value to the atmospheric delay correction module 205;
and the atmosphere delay correction module 205 is configured to remove the atmosphere delay phase value from the interference SAR unwrapped map to obtain an interferogram after the atmosphere delay correction.
By the embodiment, the method and the system for InSAR atmosphere delay correction based on ECMWF provided by the invention at least realize the following beneficial effects:
the method for correcting the atmospheric delay of the troposphere utilizes ECMWF data to calculate InSAR atmospheric delay, the used ECMWF data has the forecast data time interval of 3 hours and is updated once every 6 hours, the timeliness is strong, the spatial resolution can reach 0.125 degrees at most, and the method is the meteorological forecast data with the highest resolution at present. The invention utilizes real-time effective meteorological parameters to calculate accurate and high-precision atmospheric delay. In addition, the invention simultaneously considers the distribution rules of the atmospheric turbulence effect and the layering effect, and utilizes the layering interpolation model to carry out interpolation calculation on the delay on the spatial height distribution, thereby improving the resolution ratio and obtaining the calculation result of the atmospheric delay with higher accuracy.
Although some specific embodiments of the present invention have been described in detail by way of examples, it should be understood by those skilled in the art that the above examples are for illustrative purposes only and are not intended to limit the scope of the present invention. It will be appreciated by those skilled in the art that modifications may be made to the above embodiments without departing from the scope and spirit of the invention. The scope of the invention is defined by the appended claims.
Claims (4)
1. An InSAR atmospheric delay correction method based on ECMWF is characterized by comprising the following steps:
calculating the troposphere zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area according to the ECMWF weather forecast data and the DEM data of the research area, wherein the calculation comprises the following steps:
using a formula
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
using a formula
Calculating zenith atmosphere wet retardation, wherein ZWD is zenith atmosphere wet retardation, k' and k are refractive constants, ρwIs liquid water density, TMIs the weighted average temperature, R, of the tropospherevIs a specific atmospheric constant, PWV is the reducible moisture content;
obtaining the total zenith delay of the research area corresponding to the coarse scale of each SAR image according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area;
constructing a layered interpolation model according to the total zenith delay of each SAR image, wherein the layered interpolation model is as follows: Δ L ═ S (h) + T (x) + w, where Δ L is the total zenith delay difference of two views of SAR images, S is the stratified delay component of atmospheric delay, T is the turbulent delay component of atmospheric delay, w is the residual delay value, x is a certain point on the SAR image, and h is the elevation value at the corresponding point x;
calculating the layering delay value S of each point on the interference SAR image in the research area by adopting an interpolation mode, wherein the method comprises the following steps:
iterative function s (h) qexp { -b (h-h)min)/(hmax-hmin) Finding the optimal parameter values of q and b, wherein q and b are index parameters of the area to be interpolated, h is the elevation value of the point to be interpolated, and h is the elevation value of the point to be interpolatedmaxIs the maximum elevation value, h, of the SAR image areaminThe minimum elevation value of the SAR image area is obtained;
obtaining the layering delay value S of each point in the research area by using an interpolation method according to the optimal parameter values of q and b;
separating the layered delay values S of all points in the research area according to a layered interpolation model, and obtaining turbulence delay values T and residual delay values w of all points in the research area by adopting an inverse distance weight interpolation method;
obtaining an atmospheric delay phase value of each point on the interference SAR image according to the layering delay value S of each point in the research area and the turbulence delay value T of each point;
and removing the atmospheric delay phase value from the interference SAR unwrapping image to obtain an atmospheric delay corrected interference image.
2. The method of ECMWF-based InSAR atmospheric delay correction according to claim 1, characterized in that the total zenith delay of each SAR image of the study area is obtained from the zenith atmospheric dry delay and zenith atmospheric wet delay of the study area, comprising: ZTD is ZHD + ZWD, ZTD being the total zenith delay.
3. A system for ECMWF based InSAR atmospheric delay correction, comprising: the system comprises a data acquisition module, a total zenith delay calculation module, an interpolation model calculation module, an atmospheric delay phase value generation module and an atmospheric delay correction module; wherein the content of the first and second substances,
the data acquisition module is connected with the total zenith delay calculation module and is used for acquiring ECMWF weather forecast data and DEM data of a research area and sending the ECMWF weather forecast data and the DEM data of the research area to the total zenith delay calculation module;
the total zenith delay calculation module is connected with the difference model calculation module and used for calculating the dry delay of the zenith atmosphere and the wet delay of the zenith atmosphere of the troposphere of the research area according to the ECMWF weather forecast data and the DEM data of the research area, obtaining the total zenith delay of the coarse scale of each SAR image corresponding to the research area according to the dry delay of the zenith atmosphere and the wet delay of the zenith atmosphere of the research area, sending the total zenith delay of each SAR image of the research area to the interpolation model calculation module, wherein the total zenith delay calculation module comprises an atmosphere dry delay calculation unit and an atmosphere wet delay calculation unit,
the atmosphere dry delay calculation unit is used for adopting a formula according to ECMWF weather forecast data and DEM data of a research area
Calculating the zenith atmospheric drying delay, wherein ZHD is the zenith atmospheric drying delay and the unit is mm; ps is the surface air pressure, and the unit is hpa;
representing the dimension of a ground point in degrees; h represents the elevation of a ground point, and the unit is km;
the atmosphere moisture delay calculation unit is used for adopting a formula according to ECMWF weather forecast data and DEM data of a research area
Calculating zenith atmosphere wet retardation, wherein ZWD is zenith atmosphere wet retardation, k' and k are refractive constants, ρwIs liquid water density, TMIs the weighted average temperature, R, of the tropospherevIs a specific atmospheric constant, PWV is the reducible moisture content;
the interpolation model calculation module is connected with the atmospheric delay phase value generation module and used for constructing a layered interpolation model according to the total zenith delay of each SAR image and calculating the layered delay value S of each point in the research area and the turbulence delay value T of each point according to the layered interpolation model, and the method comprises the following steps: a model construction unit, a hierarchical delay value S calculation unit and a turbulence delay value T calculation unit, wherein,
the model building unit is used for building a layered interpolation model, and the layered interpolation model is as follows: Δ L ═ S (h) + T (x) + w, Δ L is the total zenith delay difference of two views of SAR images, S is the stratified delay portion of atmospheric delay, T is the turbulent delay portion of atmospheric delay, w is the residual delay value, x is a certain point on the SAR image, and h is the elevation value at the corresponding point x;
a hierarchical delay value S calculation unit for adoptingCalculating a layering delay value S of each point on an interference SAR image in a research area in an interpolation mode, and firstly calculating a layering delay value S according to an iteration function S (h) ═ qexp { -b (h-h)min)/(hmax-hmin) Finding the optimal parameter values of q and b, wherein q and b are index parameters of the area to be interpolated, h is the elevation value of the point to be interpolated, and h is the elevation value of the point to be interpolatedmaxIs the maximum elevation value, h, of the SAR image areaminThe minimum elevation value of the SAR image area is obtained; then obtaining the layered delay values S of all points in the research area by adopting an interpolation method through the optimal parameter values of q and b, and sending the layered delay values S of all points in the research area to an atmospheric delay phase value generation module;
the turbulence delay value T calculation unit is used for separating the layered delay value S of each point in the research area according to the layered interpolation model, obtaining the turbulence delay value T and the residual delay value w of each point in the research area by adopting an inverse distance weight interpolation method, and sending the turbulence delay value T of each point in the research area to the atmosphere delay phase value generation module;
the atmosphere delay phase value generation module is connected with the atmosphere delay correction module and used for obtaining an atmosphere delay phase value of each point on the interference SAR image according to the layering delay value S of each point in the research area and the turbulence delay value T of each point and sending the atmosphere delay phase value to the atmosphere delay correction module;
and the atmospheric delay correction module is used for removing the atmospheric delay phase value from the interference SAR unwrapping image to obtain an interference image after atmospheric delay correction.
4. The system for ECMWF based InSAR atmospheric delay correction according to claim 3, characterized in that the total zenith delay calculation module further comprises a total delay calculation module; and the total delay calculation module is used for obtaining the total zenith delay of each SAR image of the research area according to the zenith atmosphere dry delay and the zenith atmosphere wet delay of the research area, wherein ZTD is ZHD + ZWD, and ZTD is the total zenith delay.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910257685.0A CN110031841B (en) | 2019-04-01 | 2019-04-01 | ECMWF-based InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910257685.0A CN110031841B (en) | 2019-04-01 | 2019-04-01 | ECMWF-based InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110031841A CN110031841A (en) | 2019-07-19 |
| CN110031841B true CN110031841B (en) | 2021-07-23 |
Family
ID=67237134
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910257685.0A Expired - Fee Related CN110031841B (en) | 2019-04-01 | 2019-04-01 | ECMWF-based InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110031841B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112711022B (en) * | 2020-12-18 | 2022-08-30 | 中国矿业大学 | GNSS chromatography-assisted InSAR (interferometric synthetic aperture radar) atmospheric delay correction method |
| CN112816983B (en) * | 2021-01-06 | 2023-09-19 | 中南大学 | Time sequence InSAR turbulence atmosphere delay correction method based on optimized interference atlas |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103885046A (en) * | 2012-12-20 | 2014-06-25 | 河南省电力勘测设计院 | InSAR atmosphere delay correction method based on GPS |
| CN105842692A (en) * | 2016-03-17 | 2016-08-10 | 中国科学院遥感与数字地球研究所 | Atmospheric correction method during INSAR measurement |
| CN107367716A (en) * | 2017-07-04 | 2017-11-21 | 武汉大学 | A kind of high-precision satellite-borne SAR geometric calibration method |
| CN107389029A (en) * | 2017-08-24 | 2017-11-24 | 北京市水文地质工程地质大队 | A kind of surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6580998B2 (en) * | 1999-12-22 | 2003-06-17 | Rlm Software, Inc. | System and method for estimating aircraft flight delay |
| CN105629263B (en) * | 2015-12-21 | 2019-04-02 | 广州中海达卫星导航技术股份有限公司 | A kind of troposphere atmosphere delay estimation error correcting method and correction system |
| CN106526577B (en) * | 2016-10-09 | 2019-04-23 | 中国船舶重工集团公司第七一五研究所 | A kind of array shape estimation method using cooperation sound source information |
| JP6878982B2 (en) * | 2017-03-23 | 2021-06-02 | 株式会社デンソー | In-vehicle device |
| CN107839581B (en) * | 2017-09-05 | 2020-01-10 | 鲁东大学 | Automatic control feedback carriage and method based on temperature field three-dimensional imaging |
-
2019
- 2019-04-01 CN CN201910257685.0A patent/CN110031841B/en not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103885046A (en) * | 2012-12-20 | 2014-06-25 | 河南省电力勘测设计院 | InSAR atmosphere delay correction method based on GPS |
| CN105842692A (en) * | 2016-03-17 | 2016-08-10 | 中国科学院遥感与数字地球研究所 | Atmospheric correction method during INSAR measurement |
| CN107367716A (en) * | 2017-07-04 | 2017-11-21 | 武汉大学 | A kind of high-precision satellite-borne SAR geometric calibration method |
| CN107389029A (en) * | 2017-08-24 | 2017-11-24 | 北京市水文地质工程地质大队 | A kind of surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology |
Non-Patent Citations (2)
| Title |
|---|
| "不同全球对流层天顶延迟产品在中国区域的比较";李博峰 等;《同济大学学报（自 然科学版）》;20140831;第42卷(第8期);正文第2段 * |
| "综合时间序列与高程的天顶对流层延迟模型研究";周淼 等;《城市勘测》;20140430(第2期);正文第2.2节 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110031841A (en) | 2019-07-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Mateus et al. | Experimental study on the atmospheric delay based on GPS, SAR interferometry, and numerical weather model data | |
| Li et al. | Modeling atmospheric effects on InSAR with meteorological and continuous GPS observations: algorithms and some test results | |
| Tesfagiorgis et al. | Bias correction of satellite rainfall estimates using a radar-gauge product–a case study in Oklahoma (USA) | |
| CN109782282B (en) | Time series InSAR analysis method integrating troposphere atmospheric delay correction | |
| Adam et al. | Wide area persistent scatterer interferometry | |
| Benevides et al. | 4D wet refractivity estimation in the atmosphere using GNSS tomography initialized by radiosonde and AIRS measurements: results from a 1-week intensive campaign | |
| CN105182339A (en) | Method for correcting environmental influences at slope deformation monitoring on the basis of corner reflector | |
| CN110031841B (en) | ECMWF-based InSAR (interferometric synthetic aperture radar) atmospheric delay correction method and system | |
| Mateus et al. | Mapping precipitable water vapor time series from Sentinel-1 interferometric SAR | |
| Löfgren et al. | Tropospheric correction for InSAR using interpolated ECMWF data and GPS zenith total delay from the Southern California integrated GPS network | |
| Zhao et al. | An improved troposphere tomographic approach considering the signals coming from the side face of the tomographic area | |
| Yun et al. | Mitigating atmospheric effects in InSAR measurements through high-resolution data assimilation and numerical simulations with a weather prediction model | |
| Zhang et al. | A new integrated method of GNSS and MODIS measurements for tropospheric water vapor tomography | |
| Varade et al. | Snow depth in Dhundi: An estimate based on weighted bias corrected differential phase observations of dual polarimetric bi-temporal Sentinel-1 data | |
| Nitti et al. | On the use of COSMO/SkyMed data and Weather Models for interferometric DEM generation | |
| CN112711022B (en) | GNSS chromatography-assisted InSAR (interferometric synthetic aperture radar) atmospheric delay correction method | |
| Benevides et al. | Experimental GNSS tomography study in Lisbon (Portugal) | |
| CN115980751A (en) | Power law model InSAR troposphere delay correction method | |
| CN114252875B (en) | High-precision meshing method for imaging altitude data | |
| Miao et al. | Wet Tropospheric Correction Methods for Wide-Swath Altimeters | |
| Chang et al. | InSAR atmospheric distortions mitigation: GPS observations and NCEP FNL data | |
| Chang et al. | Remote sensing of atmospheric water vapor from synthetic aperture radar interferometry: case studies in Shanghai, China | |
| Du et al. | Real-time tropospheric delay map retrieval using sparse GNSS stations | |
| Wu et al. | Regression analysis of errors of sar-based dems and controlling factors | |
| Shimada et al. | Correction of atmospheric excess path delay appeared in repeat-pass SAR interferometry using objective analysis data |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210723 |