CN105842692A - Atmospheric correction method during INSAR measurement - Google Patents

Abstract

The invention discloses an atmospheric correction method during INSAR measurement. The method comprises the following steps of (1) using a WRF model to simulate a parameter needed during calculation of atmospheric delay, wherein the WRF model adopts GFS data as meteorological data; (2) calculating dry and wet atmospheric delay; (3) converting the atmospheric delay into phase delay; and (4) removing an atmospheric delay phase in an INSAR interference phase diagram. In the method, a mesoscale atmospheric model WRF is used to carry out atmospheric correction and the model can predict a weather condition and can reach a 1km resolution level, which is better than a prediction condition through using a method based on MERIS, MODIS and GPS data. By using the GFS meteorological data, a prediction data time interval is 3 hours and the data is updated for every 6b hours so that timeliness is high. A measurement result is accurate and the method is suitable for popularization and application.

Description

A Method of Atmospheric Correction in INSAR Measurement

technical field

The invention relates to the technical field of remote sensing, in particular to an atmospheric correction method in INSAR (Interferometric SAR, Synthetic Aperture Radar Interferometry) measurement.

Background technique

Atmospheric delay is one of the most important factors affecting the interferometric phase accuracy. Generally, the flight altitude of polar-orbiting SAR satellites is generally 500-800km. SAR electromagnetic wave propagation needs to pass through the ionosphere (the air layer from about 80km to 85-800km above the earth's surface) and the troposphere (the air layer from the earth's surface to 7-12km), thus being affected by the ionosphere. and tropospheric effects. In the ionosphere, the propagation of electromagnetic waves is delayed mainly due to the scattering effect of the ionized atmosphere; at the same time, due to changes in the total electron content (TEC) of the ionosphere, the propagation path of electromagnetic waves changes. In the troposphere, since the temperature, pressure and humidity of the atmosphere change with height, the atmosphere behaves as a layered medium, causing the refractive index of the atmosphere to change with height, which changes the propagation path of electromagnetic waves; in addition, due to Electromagnetic waves are bent, absorbed, reflected, and scattered by liquid and solid particles such as clouds, rainfall, and suspended particles, which can also cause signal propagation delays and path bending.

Atmospheric parameters are usually divided into wet atmospheric parameters (referring to the atmospheric pressure of water vapor in the atmosphere) and dry atmospheric parameters (that is, static atmospheric parameters, including dry atmospheric pressure and temperature), and the atmospheric delays caused by them are respectively atmospheric wet delay and atmospheric dry delay . The dry delay is relatively stable in the time domain and has large-scale variation characteristics in the air domain, and the wet delay in the troposphere dominates the overall atmospheric delay. For example, for the L-band, the radar wavelength λ is 22cm, the incident angle θ ranges from 20° to 30°, and a 10mm ZWD (Zenith Wet Delay, delay in the zenith direction caused by water vapor) error can cause 0.16 to 0.23 interferograms. Phase delay; deformation error caused by ZWD error of 10mm is 4.8mm~7.1mm. When the slope distance R is 800km and the vertical baseline is 251m, the elevation error caused by the 10mm ZWD error is 23m~36m.

At present, the more effective atmospheric correction model is atmospheric correction based on MERIS water vapor products (Xu Ji, Xie Chong et al., 2007). Atmospheric phase correction using MERIS data mainly includes three parts: obtaining the settling water vapor content, calculating the corresponding zenith wet delay when the satellite transits, and analyzing the atmospheric phase on the interference pair. Precipitable water vapor content (PerceptibleWaveVapor, PWV) and cloud information were obtained from MERIS data. The wavelengths of channels 14 and 15 of the MERIS sensor are 0.89 μm and 0.90 μm respectively, of which 0.90 μm is in the atmospheric absorption band and 0.89 μm is the atmospheric window. The ratio of reflected radiance between these two channels can be used as an indicator of the total amount of atmospheric water vapor for the MERIS sensor. The MERIS atmospheric inversion algorithm usually establishes the polynomial relationship between the ratio of channels 15 and 14 and the integrated water vapor content (Integrated Water Vapor, IWV), namely

$\mathrm{<span\; class="notranslate">\; I</span>}\mathrm{<span\; class="notranslate">\; W</span>}\mathrm{<span\; class="notranslate">\; V</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\mathrm{<span\; class="notranslate">\; lg</span>}<span\; class="notranslate">\; (</span>\frac{{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; 15</span>}}}{{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; 14</span>}}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; lg</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\frac{{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; 15</span>}}}{{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; 14</span>}}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

In the formula, I ₁₅ and I ₁₄ represent the radiation values of channels 15 and 14 of the MERIS sensor respectively; k ₀ , k ₁ and k ₂ are regression coefficients. Since the density of water vapor is 1.0×10 ³ kg/m ³ , IWV and settleable water vapor content PWV have the same value. In cloudless conditions, the theoretical accuracy of the algorithm on land is 1.6mm. MER_RR_2P provides data on atmospheric water vapor content in g/ ^cm2 . The cloud information provided by MER_RR_2P includes the type of cloud and the optical thickness of the cloud. Based on this information, it can be determined whether the cloud cover in this area is too large and whether the MERIS data of this scene is suitable for atmospheric correction. For the jth point on the SAR data of the i-th scene, according to the radar imaging geometry, the spatial three-dimensional coordinate vector Xj of the point corresponding to the ground target can be calculated. The incident angle θ ⁱ _j , according to Xj, through two-dimensional interpolation, from the MERIS data corresponding to the i-th scene ASAR data, can get the settling water vapor content PWV corresponding to the j-th point when the i-th scene ASAR data is obtained ⁱ _j and cloud information. For the pixels on the SAR image, the delay in the zenith direction (ZenithWetDelay, ZWD) caused by water vapor can be expressed by PWV, that is,

ZWD＝∏ ^-1 PWV (2)

In the formula, Π is the conversion coefficient, which is related to the real surface temperature of the area where the data is acquired. When it is determined that there is no cloud or the cloud cover is small, according to formula (2), combined with the real surface temperature T ⁱ of the area where the data is acquired, the water vapor at the jth point of ASAR data of the i-th scene can be obtained from the settling water vapor content The resulting delay in the direction of the zenith. According to formula (2), due to the measurement error of PWV, the calculation error of atmospheric phase in SAR data can be expressed as

In the formula, λ is the wavelength of the SAR sensor; θ _inc is the incident angle of the ground point corresponding to the pixel on the SAR image during imaging.

There are following deficiencies in the above-mentioned method:

1. The water vapor content data obtained by MERIS is affected by clouds, so in cloudy areas, the water vapor content of the data is not accurate.

2. This method is limited by the time resolution of MERIS data, the data receiving time is fixed, and the data time and the actual interferogram time will deviate, which will affect the calculation results.

It can be seen that the above-mentioned existing atmospheric correction method obviously still has inconvenience and defects, and needs to be further improved urgently. How to create a new atmospheric correction method in INSAR measurement with reliable results has become a goal that needs to be improved in the current industry.

Contents of the invention

The technical problem to be solved by the invention is to provide an atmospheric correction method in INSAR measurement with accurate and reliable results.

In order to solve the problems of the technologies described above, the present invention adopts the following technical solutions:

A kind of atmospheric correction method in INSAR measurement, comprises: (1) utilizes WRF model to simulate the parameter needed to calculate atmospheric delay, the meteorological data that described WRF model adopts is GFS data; (2) calculate dry, wet atmospheric delay; (3) Convert the atmospheric delay to phase delay; (4) Remove the atmospheric delay phase in the INSAR interferogram.

Further, the WRF model includes a WPS processing step and a WRF processing step, and the parameters required for calculating the atmospheric delay include temperature, humidity, and air pressure.

Further, while calculating dry and wet atmospheric delays, input DEM data, incident angles, wavelengths, resolutions, longitude and latitude ranges, and interferogram files at the same time, perform spatial interpolation, and convert the calculated delays in the zenith direction into oblique Then convert the slant distance delay (unit: cm) into phase delay (unit: rad), and draw the atmospheric delay effect map through the drawing software; subtract the processed INSAR interferometric phase map Atmospheric phase delay, the interferogram after atmospheric phase correction can be obtained.

The present invention uses the WRF model to perform InSAR atmospheric correction to obtain more real-time and effective meteorological parameters, and the calculated atmospheric delay phase is more accurate; the method uses elevation information to interpolate the delay phase on the spatial height distribution, and the result is more reliable. Calibration methods contribute significantly to the accuracy of the InSAR technique.

Description of drawings

The above is only an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

Fig. 1 is the WPS, the WRF processing program and the relationship diagram among them in the WRF model of Fig. 1;

Figure 2 is a flowchart of atmospheric correction using the WRF model;

Fig. 3 is the atmospheric delay result figure that adopts the method of the present invention to carry out atmospheric correction to the data in the Yellow River Delta region;

Figure 4 is the interferogram before atmospheric correction in the Yellow River Delta region;

Figure 5 is the interferogram after atmospheric correction in the Yellow River Delta region.

detailed description

As shown in Figure 2, for the atmospheric correction in the INSAR measurement, the present invention mainly uses the WRF model to simulate the temperature, humidity, air pressure, and geopotential height needed to calculate the atmospheric delay; use these atmospheric parameters to calculate the atmospheric dry and wet delay.

The WRF (Weather Research Forecast) model system is a new generation of mesoscale forecast model and assimilation system jointly developed and researched by scientists from many American research departments and universities. It is a flexible and perfect atmospheric simulation system, which is easy to carry, efficient, and The characteristics of parallel computing are widely used from meters to thousands of kilometers. Including: real-time numerical weather prediction, forecast research, parameterization research, etc. WRF model supports multiple meteorological data, and the meteorological data that the present invention adopts is GFS data, and the GFS (Global Forecasting System) of U.S. National Environmental Forecasting Center, its forecast data can predict the weather of 192 hours altogether in next 8 days, forecast data time The interval is 3 hours, and the resolution is 1°*1° or 0.5°*0.5°. It is updated every 6 hours, four times a day, at 06:00, 12:00, 18:00, and 00:00, respectively at 03:30, 09:30, 15:30, and 21:30 UTC. The WRF Preprocessing System (WPS) is used for real-time data processing. Its functions include: defining the simulation area; interpolating terrain data (such as terrain, soil surface and soil type) to the simulation area; interpolating data from other models ( Such as meteorological elements, etc.) to the simulation area and model coordinates.

As shown in Figure 1, the three steps of WPS include: use the geogrid module to determine the rough area (the outermost range) of a model; use ungrib to extract the meteorological element field required during the simulation from the grib data set; use metgrid to extract the The above meteorological element fields are horizontally interpolated to the model area. Two steps of WRF: run WRF data program real.exe; run WRF mode main program wrf.exe. The result files wrfplev_d* and wrfout_d* generated by the WRF model will be used as input files for subsequent atmospheric correction. The temperature, humidity, air pressure and other parameters of the meteorological parameters are extracted through the MATLAB program, and the atmospheric delay is calculated.

As shown in Figure 2, the atmospheric correction method in the INSAR measurement of the present invention, while calculating the dry and wet atmospheric delay, at the same time, it is necessary to input DEM data, incident angle, wavelength, resolution, latitude and longitude range, and interferogram files for spatial interpolation , and convert the calculated zenith-direction delay (unit: cm) into slant-range delay, and then convert the slant-range delay into phase delay (unit: rad), and draw the atmospheric delay effect map by drawing software. The atmospheric phase-corrected interferogram can be obtained by subtracting the atmospheric delay phase from the processed interferogram.

The above-mentioned atmospheric correction method of the present invention utilizes the mesoscale atmospheric model WRF to perform atmospheric correction, and the WRF model can simulate various meteorological elements required for calculating atmospheric delay, including temperature, humidity, air pressure, etc., and the model can predict weather conditions up to The resolution level of 1 km is far superior to methods based on MERIS, MODIS and GPS data. For the GFS meteorological data used, the time interval of the forecast data is 3 hours, and it is updated every 6 hours, which is more time-sensitive. This method utilizes the elevation information to interpolate the delayed phase on the spatial height distribution, and the result is more reliable.

Atmospheric corrections have been made to the data in the Yellow River Delta region using this method. The data used is ALOS-1, the data dates are 20070628, 20070813, and the resolution is 15 meters. The result of the total atmospheric delay is shown in Figure 3, and the atmospheric delay value is -1-7cm. The interferogram before atmospheric correction is shown in Figure 4, and the interferogram after atmospheric correction is shown in Figure 5.

The above is only a preferred embodiment of the present invention, and does not limit the present invention in any form. Those skilled in the art make some simple modifications, equivalent changes or modifications by using the technical content disclosed above, all of which fall within the scope of the present invention. within the scope of protection of the invention.

Claims (3)

1. the atmospheric correction method during an INSAR measures, it is characterised in that including:

(1) WRF modeling is utilized out to calculate the parameter required for atmosphere delay, the meteorological data that described WRF model uses It is GFS data；

(2) dry, damp atmosphere delay is calculated；

(3) atmosphere delay is converted to Phase delay；

(4) in INSAR interferometric phase image, atmosphere delay phase place is removed.

Atmospheric correction method in a kind of INSAR the most according to claim 1 measurement, it is characterised in that described WRF model Processing step including WPS and WRF processes step, the parameter required for calculating atmosphere delay includes temperature, humidity, air pressure, potential Highly.

Atmospheric correction method in a kind of INSAR the most according to claim 1 and 2 measurement, it is characterised in that calculating When dry, damp atmosphere postpones, it is simultaneously entered dem data, angle of incidence, wavelength, resolution, longitude and latitude scope, interferogram file, carries out Space interpolation, and the delay of the zenith direction calculated is converted into oblique distance delay, then postpone to be changed into phase place by oblique distance Postpone, draw atmosphere delay design sketch by drawing software；Processed good INSAR interferometric phase image will cut atmosphere delay Phase place, i.e. can get the interferogram after atmospheric phase correction.