CN114966681A - A Soil Moisture Estimation Method Based on Atmospheric Correction C-band InSAR Data - Google Patents

Abstract

The invention discloses a soil moisture estimation method based on atmospheric correction C-band InSAR data, comprising: denoising a SAR image; taking the denoised SAR image as an input image to calculate an interference image; As input images, APD phase estimation was performed based on Sentinel‑1 images and other external data, and APD phase calibration was performed; an analytical model was built to provide a direct relationship between soil moisture changes and interferometric phase and coherence; based on step output SAR The difference between the image interferometric phase data and the APD phase is used to obtain the remaining phase of the InSAR data, which is substituted into the analysis model as an input to obtain the soil moisture; the measured value of the soil moisture and the remaining phase are compared to evaluate the accuracy. . The invention can effectively improve the estimation accuracy of soil moisture.

Description

A Soil Moisture Estimation Method Based on Atmospheric Correction C-band InSAR Data

technical field

The invention belongs to the technical field of image processing, and in particular relates to a soil moisture estimation method based on atmospheric correction C-band InSAR data.

Background technique

Synthetic Aperture Radar (SAR) is an active earth observation system that observes the earth all day, all weather, and has a certain surface penetration capability. Therefore, the SAR system has unique advantages in disaster monitoring, environmental monitoring, ocean monitoring, resource exploration, crop yield estimation, mapping and military applications. Among them, soil moisture retrieval is a hot spot in SAR image application.

The traditional soil moisture measurement method is to use instruments or drying methods to sample and measure on the spot. This method is not only time-consuming and labor-intensive, but also has high cost and small scope of application, so it is difficult to extend to large areas. The emission of SAR (Synthetic Aperture Radar) solves these problems. The research results of soil moisture content in recent decades have shown that soil moisture content will have a great impact on the signal scattering of radar, and the two are related. Through SAR Extracting soil moisture can greatly improve the reliability and accuracy of the inversion.

The closed phase or phase triplet is the algebraic sum of the interferometric phases of the three interferograms that can be calculated from the three SAR images. Phase triples are computed for each pixel of the three co-registered SAR images. Notably, the closure phase is insensitive to surface deformation, atmospheric delays, or topographic errors, thus providing a method for estimating soil moisture. Phase disparity can be used to estimate soil moisture values. However, the model inversion of closed phase is uncertain, and even the combination of closed phase and coherent magnitude produces different results. This means that different soil moisture values combined in a triplet will produce similar closure stages.

In the absence of terrain deformation and proper removal of atmospheric phase delay (APD), direct inversion of soil moisture from interferometric phase works better than using closed phase and coherence. The main limitations of differential interferometry in measuring soil moisture are the temporal decorrelation of synthetic aperture radar (SAR) signals, especially in vegetated and agricultural areas, and the imprint of atmospheric path delay (APD).

SUMMARY OF THE INVENTION

Purpose of the invention: In order to overcome the deficiencies in the prior art, a soil moisture estimation method based on atmospheric correction C-band InSAR data is provided, the InSAR data is quantitatively linked with the soil moisture, and the atmospheric phase delay is obtained by using the SAR interference image, thereby The influence of the atmosphere on the InSAR phase is removed, and only the influence of soil moisture on the InSAR phase is retained, thereby effectively improving the inversion accuracy of soil moisture.

Technical solution: In order to achieve the above purpose, the present invention provides a soil moisture estimation method based on atmospheric correction C-band InSAR data, comprising the following steps:

S1: Denoise the SAR image;

S2: using the denoised SAR image in step S1 as the input image, process the time series of N SAR images, and calculate N-1 interference images;

S3: take the interference image output in step S2 as the input image, perform APD phase estimation according to the Sentinel-1 image and other external data, and perform APD phase calibration;

It should be noted here that Sentinel-1 consists of two polar-orbiting satellites, A and B. The sensors carried by the two satellites are Synthetic Aperture Radar (SAR), which are active microwave remote sensing satellites, and the sensors carry C-band.

S4: Build an analytical model to provide a direct relationship between soil moisture variation and interference phase and coherence;

S5: According to the difference between the SAR image interference phase data output in step S2 and the APD phase output in step S3, the remaining phase of the InSAR data is obtained, and it is substituted into the soil moisture change and interference phase in step S4 as an input Between the analytical models, the soil moisture is derived.

The present invention also includes step S6, the step S6 is specifically: converting the measured value of soil moisture into phase and coherence values according to the analysis model in step S4, and then removing the APD phase result in the interference phase of the SAR data to obtain the residual phase, and compare the two for accuracy evaluation.

Further, the step S1 is specifically: filtering the SAR image through a filter before extracting the soil moisture.

Further, the calculation method of N-1 interference images in the step S2 is:

A1: Use the daisy-chain strategy to select SAR image pairs for stack generation, and take two SAR images at time t and time t+Δt;

A2: Registration using precise orbit and external DEM;

A3: Interferogram obtained by calculation: The time difference Δt of two SAR images in the same area can be minimized to remove the phase effect of terrain changes on the SAR image. In this case, the phase difference between the two SAR images is mainly the same as that of the SAR image. The time change of APD between the acquisition times of , so the interferometric image can be obtained by relying on the phase difference between the two SAR images, thereby obtaining the interferometric phase map.

Further, the estimation method of APD in the step S3 is:

In the absence of significant surface deformation between the acquisitions of the two SAR images, the effect of APD on the interference phase φ can be assumed:

where Δφ _APD is the phase contribution caused by atmospheric delay, and λ is the radar wavelength;

Given a time series of N interferograms produced to minimize the time baseline, the APD to time t can be calculated:

和

为日期t+1和t的大气相位延迟。where i,j are pixel coordinates,

and

Atmospheric phase delay for dates t+1 and t.

Further, the method for APD phase calibration in the step S3 is:

The present invention proposes a four-parameter linear regression model, introduces the low-frequency spatial distribution information of APD, and calibrates it with a group of global navigation satellite system stations;

The system equation is:

是像素(i,j)的测地坐标。where a ₁ ,...,a ₄ are the parameters to be estimated,

are the geodesic coordinates of pixel (i,j).

Further, the expression of the analysis model in the step S4 is as follows:

Assuming that the differential propagation of electromagnetic waves in soil is related to the effect of soil moisture level on the vertical wavenumber of the surface layer, the soil is modeled as a uniform lossy dielectric layer with complex permittivity ε′, and the vertical wavenumber k′ _Z is:

where ω and μ are the angular velocity and dielectric constant, and k _x is the wave number in the line of sight;

The interference phase is:

The coherence is:

Further, the acquisition method of soil moisture in the step S5 is:

The estimated soil moisture is obtained by minimizing a cost function representing the misfit between the measured interferometric phase and the phase predicted by the model, as a function of soil moisture:

Among them, N ₁ is the number of interferograms, and φ(mυ _i , mυ _j ) and γ(mυ _i , mυ _j ) are the phase and coherence values predicted by the model, respectively.

Further, the specific method of precision evaluation in the step S6 is:

Firstly, according to the analysis model between soil moisture change and interference phase and coherence, the measured value of soil moisture is converted into phase and coherence value, so as to avoid obtaining wrong soil moisture when inverting the phase; The APD results were removed from the interference phase to obtain the residual phase, and the soil moisture inversion accuracy was obtained by comparing the residual phase with the phase and coherence values of the soil moisture measured data.

The invention establishes the connection between InSAR data and soil moisture, uses the C-band Synthetic Aperture Radar (SAR) Interferometry (InSAR) technology to process the C-band Synthetic Aperture Radar (SAR) image, generates a corrected soil moisture map, and uses The Atmospheric Phase Delay (APD) map obtained from the time series of Sentinel-1 interferograms is used to decompose the images of APD and soil moisture on Sentinel-1 interferograms, thereby effectively improving the accuracy of soil moisture inversion.

Beneficial effects: Compared with the prior art, the present invention quantitatively links InSAR data with soil moisture, and uses SAR interference images to obtain atmospheric phase delay, thereby removing the influence of atmosphere on InSAR phase, and only retaining soil moisture on InSAR phase Therefore, the inversion accuracy of soil moisture can be effectively improved.

Description of drawings

FIG. 1 is a schematic flow chart of the method of the present invention.

Detailed ways

Below in conjunction with the accompanying drawings and specific embodiments, the present invention will be further clarified. It should be understood that these embodiments are only used to illustrate the present invention and not to limit the scope of the present invention. Modifications of equivalent forms all fall within the scope defined by the appended claims of this application.

The present invention provides a soil moisture estimation method based on atmospheric correction C-band InSAR data, as shown in FIG. 1 , which includes the following steps:

S1: Denoise the SAR image:

The SAR image is filtered through a filter before extracting soil moisture.

S2: Take the denoised SAR image in step S1 as the input image, process the time series of N SAR images, and calculate N-1 interference images:

The calculation method of N-1 interference images is:

A1: Use the daisy-chain strategy to select SAR image pairs for stack generation, and take two SAR images at time t and time t+Δt;

A2: Registration using precise orbit and external DEM;

A3: Interferogram obtained by calculation: The time difference Δt of two SAR images in the same area can be minimized to remove the phase effect of terrain changes on the SAR image. In this case, the phase difference between the two SAR images is mainly the same as that of the SAR image. The time change of APD between the acquisition times of , so the interferometric image can be obtained by relying on the phase difference between the two SAR images, thereby obtaining the interferometric phase map.

S3: Take the interference image output in step S2 as the input image, perform APD phase estimation according to Sentinel-1 image and other external data, and perform APD phase calibration:

The estimation method of APD is:

In the absence of significant surface deformation between the acquisitions of the two SAR images, the effect of APD on the interference phase φ can be assumed:

where Δφ _APD is the phase contribution caused by atmospheric delay, and λ is the radar wavelength;

Given a time series of N interferograms produced to minimize the time baseline, the APD to time t can be calculated:

和

为日期t+1和t的大气相位延迟。where i,j are pixel coordinates,

and

Atmospheric phase delay for dates t+1 and t.

The method of APD phase calibration is:

The present invention proposes a four-parameter linear regression model, introduces the low-frequency spatial distribution information of APD, and calibrates it with a group of global navigation satellite system stations;

The system equation is:

是像素(i,j)的测地坐标。where a ₁ ,...,a ₄ are the parameters to be estimated,

are the geodesic coordinates of pixel (i,j).

S4: Build an analytical model to provide a direct relationship between soil moisture changes and interference phase and coherence:

The analytical model is expressed as follows:

Assuming that the differential propagation of electromagnetic waves in soil is related to the effect of soil moisture level on the vertical wavenumber of the surface layer, the soil is modeled as a uniform lossy dielectric layer with complex permittivity ε′, and the vertical wavenumber k′ _Z is:

where ω and μ are the angular velocity and dielectric constant, and k _x is the wave number in the line of sight;

The interference phase is:

The coherence is:

S5: According to the difference between the SAR image interference phase data output in step S2 and the APD phase output in step S3, the remaining phase of the InSAR data is obtained, and it is substituted into the soil moisture change and interference phase in step S4 as an input Between the analytical models, the soil moisture is derived:

The method of obtaining soil moisture is:

The estimated soil moisture is obtained by minimizing a cost function representing the misfit between the measured interferometric phase and the phase predicted by the model, as a function of soil moisture:

Among them, N ₁ is the number of interferograms, and φ(mυ _i , mυ _j ) and γ(mυ _i , mυ _j ) are the phase and coherence values predicted by the model, respectively.

S6: Convert the measured value of soil moisture into phase and coherence values according to the analysis model in step S4, and then remove the APD phase result in the interferometric phase of the SAR data to obtain the remaining phase, and compare the two for accuracy evaluation.

The specific method of accuracy evaluation is as follows:

Firstly, according to the analysis model between soil moisture change and interference phase and coherence, the measured value of soil moisture is converted into phase and coherence value, so as to avoid obtaining wrong soil moisture when inverting the phase; The APD results were removed from the interference phase to obtain the residual phase, and the soil moisture inversion accuracy was obtained by comparing the residual phase with the phase and coherence values of the soil moisture measured data.

Based on the above scheme, the above scheme is applied by example in this embodiment, and the specific process is as follows:

Step 1: Select a sliding window of size N*N, and perform mean filtering on the input SAR image according to the principle of spatial filtering.

Step 2: Use daisy-chain strategy to select SAR image pairs for stack generation, take two SAR images at time t and time t+Δt, Δt is 6 days; use precise orbit and external DEM for registration; interferogram calculation; Earth curvature and terrain effect removal; coherence estimation using azimuth 3, distance 3 windows; interferogram and coherent terrain correction and geocoding for WGS84 UTM coordinate reference system. The output pixel resolution is about 20m.

Step 3: Estimate Atmospheric Phase Delay (APD) from Sentinel-1 images and other external data:

In the absence of significant surface deformation between the acquisitions of the two SAR images, it can be assumed that APD is the dominant effect on the interferometric phase φ:

where Δφ _APD is the phase contribution caused by atmospheric delay, and λ is the radar wavelength;

Given a time series of N interferograms produced to minimize the time baseline, the APD at time t can be calculated:

和

为日期t+1和t的大气相位延迟。where i,j are pixel coordinates,

and

Atmospheric phase delay for dates t+1 and t.

From this, a four-parameter linear regression model is proposed, which introduces the low-frequency spatial distribution information of APD and calibrates it with a set of GNSS stations.

The system equation is:

是像素(i,j)的测地坐标。这四个参数是用已知的APD值估计的，所以至少需要四个全球导航卫星系统台站。where a ₁ ,...,a ₄ are the parameters to be estimated,

are the geodesic coordinates of pixel (i,j). These four parameters are estimated with known APD values, so at least four GNSS stations are required.

Step 4: Assuming that the differential propagation of electromagnetic waves in the soil is related to the effect of soil moisture level on the vertical wavenumber of the surface layer, the soil is modeled as a uniform lossy dielectric layer with complex permittivity ε′.

The vertical wave number k′ _Z is:

Among them, ω and μ are the angular velocity and dielectric constant, and k _x is the wave number in the line of sight.

The interference phase is:

The coherence is:

Step 5: The estimated soil moisture is obtained by minimizing a cost function representing the misfit between the measured interference phase and the model predicted phase, which is a function of soil moisture:

Among them, N ₁ is the number of interferograms, and φ(mυ _i , mυ _j ) and γ(mυ _i , mυ _j ) are the phase and coherence values predicted by the model, respectively.

Step 6: Accuracy evaluation: Compare the derived residual phase with the soil moisture measurement as follows, first convert the measured soil moisture into phase and coherence values according to the analytical model between soil moisture changes and the interferometric phase and coherence, This avoids getting the wrong soil moisture when inverting the phase; computes APD results with different multi-windows and removes them from the interferometric phase to get the residual phase, which uses a coherence-based space In this way, pixels with coherence less than 0.2 are not considered in the multi-window calculation; the soil moisture inversion accuracy is obtained by comparing the residual phase and the phase and coherence values of the soil moisture measured data.

This embodiment also provides a soil moisture estimation system based on atmospheric corrected C-band InSAR data, the system includes a network interface, a memory and a processor; wherein the network interface is used for sending and receiving information with other external network elements. , realize the reception and transmission of signals; the memory is used to store the computer program instructions that can be run on the processor; the processor is used to execute the steps of the above consensus method when running the computer program instructions.

This embodiment also provides a computer storage medium, where a computer program is stored in the computer storage medium, and the method described above can be implemented when a processor executes the computer program. The computer-readable medium may be considered tangible and non-transitory. Non-limiting examples of non-transitory tangible computer-readable media include non-volatile memory circuits (eg, flash memory circuits, erasable programmable read-only memory circuits, or masked read-only memory circuits), volatile memory circuits (eg, static random access memory circuits or dynamic random access memory circuits), magnetic storage media such as analog or digital magnetic tapes or hard drives, and optical storage media such as CD, DVD or Blu-ray discs, among others. A computer program includes processor-executable instructions stored on at least one non-transitory tangible computer-readable medium. The computer program may also include or rely on stored data. Computer programs may include a basic input/output system (BIOS) that interacts with the hardware of the special purpose computer, device drivers that interact with specific devices of the special purpose computer, one or more operating systems, user applications, background services, background applications, etc. .

As will be appreciated by those skilled in the art, the embodiments of the present application may be provided as a method, a system, or a computer program product. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) having computer-usable program code embodied therein.

The present application is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the present application. It will be understood that each flow and/or block in the flowchart illustrations and/or block diagrams, and combinations of flows and/or blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to the processor of a general purpose computer, special purpose computer, embedded processor or other programmable data processing device to produce a machine such that the instructions executed by the processor of the computer or other programmable data processing device produce Means for implementing the functions specified in a flow or flow of a flowchart and/or a block or blocks of a block diagram.

These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory result in an article of manufacture comprising instruction means, the instructions The apparatus implements the functions specified in the flow or flow of the flowcharts and/or the block or blocks of the block diagrams.

These computer program instructions can also be loaded on a computer or other programmable data processing device to cause a series of operational steps to be performed on the computer or other programmable device to produce a computer-implemented process such that The instructions provide steps for implementing the functions specified in the flow or blocks of the flowcharts and/or the block or blocks of the block diagrams.

Claims (9)

1. A soil humidity estimation method based on atmospheric correction C-band InSAR data is characterized by comprising the following steps:

s1: denoising the SAR image;

s2: processing the time sequence of the N SAR images by taking the denoised SAR image in the step S1 as an input image, and calculating N-1 interference images;

s3: with the interference image output in step S2 as an input image, performing APD phase estimation from the Sentinel-1 image and other external data, and performing APD phase calibration;

s4: establishing an analysis model for providing a direct relation between soil moisture change and interference phase and coherence;

s5: and obtaining the residual phase of the InSAR data according to the difference between the SAR image interference phase data output in the step S2 and the APD phase output in the step S3, and substituting the residual phase of the InSAR data as input into an analysis model between the soil moisture change and the interference phase in the step S4 to obtain the soil humidity.

2. The method for estimating soil humidity based on atmospheric correction C-band InSAR data according to claim 1, further comprising a step S6, wherein the step S6 specifically comprises: and converting the measured value of the soil moisture into a phase and a coherent value according to the analysis model in the step S4, removing an APD phase result in the SAR data interference phase to obtain a residual phase, comparing the residual phase and the residual phase, and performing precision evaluation.

3. The method for estimating soil humidity based on atmospheric correction C-band InSAR data according to claim 1, wherein the step S1 is specifically as follows: and filtering the SAR image by a filter before extracting the soil humidity.

4. The soil humidity estimation method based on atmospheric correction C-band InSAR data according to claim 1, characterized in that the calculation method of the N-1 interference images in the step S2 is as follows:

a1: selecting SAR image pairs by using a daisy chain strategy to perform stack generation, and taking two SAR images at the time t and the time t + delta t;

a2: registration using precision orbit and external DEM;

a3: and calculating to obtain an interference pattern.

5. The soil humidity estimation method based on atmospheric correction C-band InSAR data according to claim 1, wherein the estimation method of APD in step S3 is:

assuming the effect of APD on the interference phase phi:

wherein, is _APD For the phase contribution caused by atmospheric delay, λ is the radar wavelength;

given a time series of N interferograms generated to minimize the time base, the APD to time t can be calculated:

where i, j are the pixel coordinates,

and

atmospheric phase delay for dates t +1 and t.

6. The soil humidity estimation method based on atmospheric correction C-band InSAR data according to claim 5, wherein the APD phase calibration method in step S3 is:

providing a four-parameter linear regression model, introducing low-frequency spatial distribution information of APD, and calibrating the APD by using a group of global navigation satellite system stations;

the system equation is:

wherein a is ₁ ,...,a ₄ Is the parameter to be estimated and is,

is the geodetic coordinates of pixel (i, j).

7. The method for estimating soil humidity based on atmospheric correction C-band InSAR data as claimed in claim 1, wherein the expression of the analytical model in step S4 is as follows:

assuming that differential propagation of electromagnetic waves in soil is related to the effect of soil moisture level on the vertical wave number of surface layers, the soil is modeled as a uniform lossy dielectric layer having a complex dielectric constant ε ', the vertical wave number k' _Z Comprises the following steps:

wherein ω and μ are angular velocity, dielectric constant, k _x The wave number in the direction of the line of sight;

the interference phase is:

the coherence is as follows:

8. the soil humidity estimation method based on the atmospheric correction C-band InSAR data as claimed in claim 1, wherein the soil humidity acquisition method in the step S5 is as follows:

the estimated soil moisture is obtained by minimizing a cost function representing the mismatch between the measured interferometric phase and the phase predicted by the model, the soil moisture function being:

wherein N is ₁ Is the number of interferograms, phi (m upsilon) _i ,mυ _j ) And gamma (m upsilon) _i ,mυ _j ) Respectively, the phase value and coherence value predicted by the model.

9. The soil humidity estimation method based on atmospheric correction C-band InSAR data according to claim 2, characterized in that the concrete method of precision evaluation in step S6 is:

firstly, converting measured values of soil moisture into phases and coherence values according to an analysis model between soil moisture change and interference phases and coherence, calculating APD (avalanche photo diode) results with different multiple windows, removing the APD results from the interference phases to obtain residual phases, and comparing the phases and the coherence values of the residual phases and measured data of the soil moisture to obtain soil moisture inversion accuracy.

Images