阅读说明：本技术 一种植被覆盖区域的土壤湿度估计方法及装置 (Soil humidity estimation method and device for vegetation coverage area ) 是由 严颂华 谢媛 陈能成 于 2019-11-13 设计创作，主要内容包括：本发明公开了一种植被覆盖区域的土壤湿度估计方法及装置,属于电子、信息及光学遥感等领域,首先,借助星载GNSS-R反射计数据构建总体后向反射系数和DDM反射功率关系；利用光学影像数据求解植被覆盖区域植被反射系数和衰减系数；融合星载GNSS-R反射计数据总体后向反射系数和DDM反射功率的关系以及光学数据中植被信息,求解植被覆盖区域纯粹土壤表面反射系数和DDM反射功率关系；利用多个控制点土壤湿度的实测数据建立土壤湿度和DDM反射功率回归模型；以某点GNSS-R数据DDM反射功率代入回归方程反演土壤湿度。本发明的实施是一次星载GNSS-R反射计数据与多光谱遥感数据影像共同反演植被覆盖区域土壤湿度的有效融合。(The invention discloses a soil humidity estimation method and a soil humidity estimation device for a vegetation coverage area, which belong to the fields of electronics, information, optical remote sensing and the like, and firstly, a relation between an overall backward reflection coefficient and DDM (distributed data model) reflection power is constructed by means of satellite-borne GNSS-R reflectometer data; solving vegetation reflection coefficients and attenuation coefficients of the vegetation coverage area by using the optical image data; the method comprises the steps of fusing the relation between the total retroreflection coefficient and the DDM reflected power of satellite-borne GNSS-R reflectometer data and vegetation information in optical data, and solving the relation between the pure soil surface reflection coefficient and the DDM reflected power of a vegetation coverage area; establishing a soil humidity and DDM reflected power regression model by using measured data of soil humidity of a plurality of control points; and substituting the DDM reflected power of certain point GNSS-R data into a regression equation to invert the soil humidity. The implementation of the method is effective fusion of the data of the satellite-borne GNSS-R reflectometer and the multispectral remote sensing data image together for inverting the soil humidity of the vegetation coverage area.)

1. A method for estimating soil moisture in a vegetation covered area, comprising:

(1) constructing a relation between a soil total back reflection coefficient and DDM reflection power based on satellite-borne GNSS-R reflectometer data;

(2) obtaining vegetation reflection coefficients and attenuation coefficients of vegetation coverage areas by using optical image data;

(3) based on the relation between the soil total backward reflection coefficient and the DDM reflection power, eliminating the influence of the vegetation reflection coefficient and the attenuation coefficient to obtain the relation between the pure soil surface reflection coefficient of the vegetation coverage area and the DDM reflection power;

(4) according to the relation between the pure soil surface reflection coefficient and the DDM reflection power of the vegetation covered area and the linear relation between soil humidity and the pure soil surface reflection coefficient, the relation between the soil humidity value and the DDM reflection power of the vegetation covered area is constructed, and then according to the actually measured soil humidity of a plurality of control points and the DDM reflection power measured value of the corresponding point, an unknown parameter in the relation between the soil humidity value and the DDM reflection power of the vegetation covered area is obtained, so that a regression model for estimating the soil humidity by using the DDM reflection power measured value is obtained;

(5) and selecting the reflection power of a target point in satellite-borne GNSS-R reflectometer data of the vegetation coverage area, and performing inversion to the soil humidity of the vegetation coverage area based on the regression model.

2. The method of claim 1, wherein the method is performed by

3. The method of claim 1 or 2, wherein step (2) comprises:

(2.1) preprocessing the optical image data to obtain spectral image data of the region of interest;

(2.2) acquiring a normalized water index of each pixel point in the spectral image data of the region of interest, and removing a water body part in the spectral image data of the region of interest according to the relation between each normalized water index and a first preset threshold value to obtain spectral image data after the water body is removed;

(2.3) acquiring a normalized vegetation index of each pixel point in the spectral image data after the water body is removed, and determining a vegetation coverage area according to the relation between each normalized vegetation index and a second preset threshold;

and (2.4) obtaining the vegetation water content of the vegetation coverage area according to the normalized vegetation index of each pixel point in the vegetation coverage area, and further obtaining the vegetation reflection coefficient and the attenuation coefficient of the vegetation coverage area according to the vegetation water content.

4. The method of claim 2, wherein δ is defined by_veg(θ)＝A*M_vegcos(θ)[1-τ²(θ)]Determining vegetation reflectance of a vegetation coverage areaFrom

5. The method of claim 4, wherein the method is performed by

6. The method of claim 5, wherein step (4) comprises:

(4.1) preparation of a copolymer of_soil(θ)＝a₁mv+a₂Determination of soil moisture mv and pure soil surface reflection coefficient delta_soil(theta) linear relationship therebetween, wherein a₁、a₂Is a coefficient;

(4.2) by Pr ═ 20lgV_in[τ²(θ)(a₁mv+a₂)+δ_veg(θ)]Determining a relationship between a soil humidity value of a vegetation covered area and the DDM reflected power;

(4.3) actually measuring the soil humidity mv of the control points of the m vegetation coverage areas_iAnd then sequentially reading the reflection power Pr of the corresponding control point on the GNSS-R reflection data track_iThen the vegetation reflection coefficient delta of the corresponding vegetation coverage area is obtained_veg,i(theta) and attenuation coefficient

(4.4) according to mv_i、Pr_i、δ_veg,i(theta) and

7. The method of claim 6, wherein the method is performed by

8. The method of claim 6 or 7, wherein step (5) comprises:

(5.1) selecting the reflection power of a target point in satellite-borne GNSS-R reflectometer data of a vegetation coverage area;

(5.2) obtaining a normalized vegetation index of the target point in a vegetation coverage area by using the optical image data, and further obtaining a vegetation reflection coefficient and an attenuation coefficient of the target point according to the normalized vegetation index of the target point;

and (5.3) obtaining the soil humidity of the vegetation coverage area by utilizing the regression model according to the reflection power of the target point, the vegetation reflection coefficient and the attenuation coefficient of the target point.

9. An apparatus for estimating soil moisture in a vegetation covered area, comprising:

the first relation determining module is used for constructing the relation between the soil total retroreflection coefficient and the DDM reflected power based on the satellite-borne GNSS-R reflectometer data;

the second relation determination module is used for obtaining vegetation reflection coefficients and attenuation coefficients of the vegetation coverage area by using the optical image data;

the third relation determination module is used for eliminating the influence of the vegetation reflection coefficient and the attenuation coefficient based on the relation between the soil total backward reflection coefficient and the DDM reflection power to obtain the relation between the pure soil surface reflection coefficient of the vegetation covered area and the DDM reflection power;

the regression model determining module is used for constructing a relation between a soil humidity value of the vegetation coverage area and the DDM reflected power according to the relation between the pure soil surface reflection coefficient and the DDM reflected power of the vegetation coverage area and the linear relation between the soil humidity and the pure soil surface reflection coefficient, and obtaining an unknown parameter in the relation between the soil humidity value of the vegetation coverage area and the DDM reflected power according to the actually measured soil humidity of a plurality of control points and the DDM reflected power measured value of the corresponding point, so that a regression model for estimating the soil humidity by using the DDM reflected power measured value is obtained;

and the inversion module is used for selecting the reflection power of a target point in the satellite-borne GNSS-R reflectometer data of the vegetation coverage area and inverting the soil humidity of the vegetation coverage area based on the regression model.

Technical Field

The invention belongs to the technical field of electronics, information and remote sensing, and particularly relates to a soil humidity estimation method and device based on fusion of a satellite-borne GNSS-R reflectometer and multispectral data.

Background

Soil humidity is an important state parameter of global water circulation, and plays an important role in agricultural irrigation, environmental research, natural disaster early warning, water and soil treatment and the like.

The current methods for measuring soil moisture include the following: the traditional soil moisture measuring method comprises a gravimetric method, a time domain reflection method, a capacitance sensor, a neutron instrument, a resistivity instrument, a heat pulse sensor, an optical fiber sensor and the like. But the defects are fixed-point measurement, small measurement range and higher labor cost.

The ground measurement aspect comprises a newly-sent GNSS-IR soil humidity inversion method based on GPS multi-satellite three-frequency data fusion, and a method for performing joint inversion by weighting and fusing GPS multi-satellite L1, L2 and L5 frequency band data is provided by utilizing the difference and complementarity of different tracks and different frequencies.

In the aspect of monitoring soil humidity by using a satellite remote sensing technology, optical, thermal infrared and microwave remote sensing methods are mainly used, and the optical remote sensing method is mainly used for estimating the water content of soil by using spectral reflection characteristics, but is easily influenced by weather. The thermal infrared remote sensing method is used for inverting the soil humidity by utilizing the thermal characteristics of soil, but in places with high vegetation coverage, because vegetation covers soil information, the accuracy of estimating the soil water content can be influenced, and therefore the method is only suitable for monitoring the soil water content in bare soil and vegetation sparse areas. The remote sensing soil humidity method utilizing the synthetic aperture radar SAR satellite mainly adopts the method that the soil moisture is inverted through the power ratio of a soil surface reflection signal to a transmission signal. Although the problem of measuring the soil humidity under vegetation coverage is solved, the method is limited by the re-return period of the SAR satellite, the data continuity cannot be guaranteed, and one field of data can be obtained in the same place for several days. Similarly, in patent CN103940834A, a method for measuring soil humidity by using synthetic aperture radar technology is disclosed, and soil humidity also needs to be inverted by using polarized synthetic aperture radar complex images of an observation area obtained by repeatedly observing at two different times t1 and t 2.

At present, ground remote sensing by using a reflection signal of a navigation signal is a research hotspot, and when a reflection signal receiver is installed on a satellite, the reflection signal receiver is called as an on-board GNSS-R reflectometer. The navigation signal utilized by the satellite-borne GNSS-R reflectometer is in an L wave band, is a microwave signal, has penetrability, can not be influenced by weather conditions such as cloud layers, rain, snow and the like, and can penetrate through vegetation, so that the satellite-borne GNSS-R reflectometer can be used as a favorable means for inverting the soil humidity under vegetation coverage. However, the vegetation parameters are in many cases unknown, which causes vegetation to interfere with the results of soil moisture inversion.

In summary, an effective method for estimating soil moisture in vegetation covered areas is needed.

Disclosure of Invention

Aiming at the defects or improvement requirements of the prior art, the invention provides a method and a device for estimating soil humidity of a vegetation covered area, so that the technical problem that the existing method for estimating the soil humidity of the vegetation covered area cannot effectively estimate the soil humidity is solved.

To achieve the above object, according to one aspect of the present invention, there is provided a soil moisture estimation method of a vegetation coverage area, including:

(1) constructing a relation between a soil total back reflection coefficient and DDM reflection power based on satellite-borne GNSS-R reflectometer data;

(2) obtaining vegetation reflection coefficients and attenuation coefficients of vegetation coverage areas by using optical image data;

(3) based on the relation between the soil total backward reflection coefficient and the DDM reflection power, eliminating the influence of the vegetation reflection coefficient and the attenuation coefficient to obtain the relation between the pure soil surface reflection coefficient of the vegetation coverage area and the DDM reflection power;

(4) according to the relation between the pure soil surface reflection coefficient and the DDM reflection power of the vegetation covered area and the linear relation between soil humidity and the pure soil surface reflection coefficient, the relation between the soil humidity value and the DDM reflection power of the vegetation covered area is constructed, and then according to the actually measured soil humidity of a plurality of control points and the DDM reflection power measured value of the corresponding point, an unknown parameter in the relation between the soil humidity value and the DDM reflection power of the vegetation covered area is obtained, so that a regression model for estimating the soil humidity by using the DDM reflection power measured value is obtained;

(5) and selecting the reflection power of a target point in satellite-borne GNSS-R reflectometer data of the vegetation coverage area, and performing inversion to the soil humidity of the vegetation coverage area based on the regression model.

Preferably, is prepared from

Constructing a relation between soil total retroreflection coefficient and DDM reflected power, wherein delta_totalRepresenting the total retroreflectance coefficient of the soil, Pr representing the DDM reflected power, V_inThe microwave voltage at the reflection point of the microwave incidence mirror is shown.

Preferably, step (2) comprises:

(2.1) preprocessing the optical image data to obtain spectral image data of the region of interest;

(2.2) acquiring a normalized water index of each pixel point in the spectral image data of the region of interest, and removing a water body part in the spectral image data of the region of interest according to the relation between each normalized water index and a first preset threshold value to obtain spectral image data after the water body is removed;

(2.3) acquiring a normalized vegetation index of each pixel point in the spectral image data after the water body is removed, and determining a vegetation coverage area according to the relation between each normalized vegetation index and a second preset threshold;

and (2.4) obtaining the vegetation water content of the vegetation coverage area according to the normalized vegetation index of each pixel point in the vegetation coverage area, and further obtaining the vegetation reflection coefficient and the attenuation coefficient of the vegetation coverage area according to the vegetation water content.

Preferably, by_veg(θ)＝A*M_vegcos(θ)[1-τ²(θ)]Determining vegetation coverageVegetation reflection coefficient of area composed of

Determining the attenuation coefficient of the vegetation coverage area, wherein delta_veg(theta) is the vegetation reflection coefficient, tau²(theta) is the attenuation coefficient of radar waves penetrating the vegetation layer, theta represents the incident angle of the radar sensor, M_vegFor vegetation moisture content, A, B is a vegetation type empirical parameter.

Preferably, is prepared from

Determining a relationship between a pure soil surface reflection coefficient of a vegetation covered area and the DDM reflected power, wherein δ_soilAnd (theta) is the pure soil surface reflection coefficient.

Preferably, step (4) comprises:

(4.1) preparation of a copolymer of_soil(θ)＝a₁mv+a₂Determination of soil moisture mv and pure soil surface reflection coefficient delta_soil(theta) linear relationship therebetween, wherein a₁、a₂Is a coefficient;

(4.2) by Pr ═ 20lgV_in[τ²(θ)(a₁mv+a₂)+δ_veg(θ)]Determining a relationship between a soil humidity value of a vegetation covered area and the DDM reflected power;

(4.3) actually measuring the soil humidity mv of the control points of the m vegetation coverage areas_iAnd then sequentially reading the reflection power Pr of the corresponding control point on the GNSS-R reflection data track_iThen the vegetation reflection coefficient delta of the corresponding vegetation coverage area is obtained_veg,i(theta) and attenuation coefficient

Wherein the value range of i is 1 to m;

(4.4) according to mv_i、Pr_i、δ_veg,i(theta) and

obtaining the vegetation coverage by adopting a least square methodParameter a in the relationship between regional soil moisture value and the DDM reflected power₁、a₂And V_inAnd further obtaining a regression model for estimating the soil humidity by using the DDM reflected power measured value.

Preferably, is prepared fromA regression model for estimating soil moisture using DDM reflected power measurements is obtained.

Preferably, step (5) comprises:

(5.1) selecting the reflection power of a target point in satellite-borne GNSS-R reflectometer data of a vegetation coverage area;

(5.2) obtaining a normalized vegetation index of the target point in a vegetation coverage area by using the optical image data, and further obtaining a vegetation reflection coefficient and an attenuation coefficient of the target point according to the normalized vegetation index of the target point;

and (5.3) obtaining the soil humidity of the vegetation coverage area by utilizing the regression model according to the reflection power of the target point, the vegetation reflection coefficient and the attenuation coefficient of the target point.

According to another aspect of the present invention, there is provided a soil moisture estimation device for a vegetation coverage area, comprising:

the first relation determining module is used for constructing the relation between the soil total retroreflection coefficient and the DDM reflected power based on the satellite-borne GNSS-R reflectometer data;

the second relation determination module is used for obtaining vegetation reflection coefficients and attenuation coefficients of the vegetation coverage area by using the optical image data;

the third relation determination module is used for eliminating the influence of the vegetation reflection coefficient and the attenuation coefficient based on the relation between the soil total backward reflection coefficient and the DDM reflection power to obtain the relation between the pure soil surface reflection coefficient of the vegetation covered area and the DDM reflection power;

the regression model determining module is used for constructing a relation between a soil humidity value of the vegetation coverage area and the DDM reflected power according to the relation between the pure soil surface reflection coefficient and the DDM reflected power of the vegetation coverage area and the linear relation between the soil humidity and the pure soil surface reflection coefficient, and obtaining an unknown parameter in the relation between the soil humidity value of the vegetation coverage area and the DDM reflected power according to the actually measured soil humidity of a plurality of control points and the DDM reflected power measured value of the corresponding point, so that a regression model for estimating the soil humidity by using the DDM reflected power measured value is obtained;

and the inversion module is used for selecting the reflection power of a target point in the satellite-borne GNSS-R reflectometer data of the vegetation coverage area and inverting the soil humidity of the vegetation coverage area based on the regression model.

In general, compared with the prior art, the above technical solution contemplated by the present invention can achieve the following beneficial effects:

the invention provides a method for inverting the soil humidity of a vegetation coverage area by utilizing a satellite-borne GNSS-R reflectometer data and optical remote sensing image fusion mode. Because the traditional SAR method has the defects of long return cycle and the like, the method of the invention exerts the advantages of large number of navigation satellites and global coverage, simultaneously eliminates the vegetation influence by utilizing multispectral image information under the condition that the vegetation is covered, and thereby reversely shows the regional soil humidity. The invention can be applied to other fields such as agriculture and the like due to the characteristics of wide range coverage and high efficiency, and the implementation of the invention is effective fusion of one-time satellite-borne GNSS-R reflectometer data and multispectral remote sensing data image for jointly inverting the soil humidity of a vegetation coverage area.

Drawings

FIG. 1 is a diagram of a model provided by an embodiment of the present invention;

FIG. 2 is a schematic flow chart of a method provided by an embodiment of the present invention;

fig. 3 is a flowchart illustrating a method for solving vegetation reflection coefficients and attenuation coefficients in a vegetation coverage area according to an embodiment of the present invention;

fig. 4 is a schematic diagram of soil humidity of a vegetation coverage area obtained by utilizing data of a satellite-borne GNSS-R reflectometer according to an embodiment of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In addition, the technical features involved in the embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

The invention fully utilizes the L wave band microwave signal of the satellite-borne GNSS-R reflectometer data to have penetrability, can not be influenced by weather conditions such as cloud layers, rain, snow and the like, can also penetrate the vegetation characteristic and the multispectral remote sensing data wave band information to solve the vegetation parameter method, deducts vegetation attenuation from the reflectometer data, and combines the two to jointly invert the soil humidity under the vegetation coverage.

The invention provides a method for combining satellite-borne GNSS-R reflectometer data with a multispectral high-score data image, on one hand, a calculation model of total backward reflection coefficient and reflection power is constructed by utilizing the satellite-borne GNSS-R reflectometer data, on the other hand, a vegetation coverage area is distinguished by utilizing spectral information in an optical remote sensing image, the vegetation reflection coefficient and attenuation coefficient are calculated, a pure soil surface reflection coefficient under vegetation is obtained by means of a water cloud model, finally, a regression model of the pure soil humidity and the reflection power under vegetation coverage is constructed by utilizing an influence model of the soil humidity on the pure soil surface reflection coefficient, and the soil humidity is inverted, and a model diagram of the method is shown in figure 1.

As shown in fig. 2, which is a schematic flow chart of a method provided in an embodiment of the present invention, the method shown in fig. 2 includes the following steps:

s1: constructing a relation between an overall backward reflection coefficient and DDM reflection power by using satellite-borne GNSS-R reflectometer data;

specifically, DDM reflected power is obtained from Delay-Doppler Map (Delay-Doppler Map) data of a satellite-borne GNSS-R reflectometer, and direct voltage is set to be an unknown number, so that a formula for calculating an overall backward reflection coefficient from a DDM reflected power measurement value is obtained.

S2: utilizing the multispectral image data to solve vegetation reflection coefficients and attenuation coefficients of vegetation coverage areas;

specifically, the optical image is preprocessed to obtain a normalized vegetation index, and then the vegetation index is used for distinguishing a vegetation coverage area. And solving the vegetation reflection coefficient and the attenuation coefficient by means of the optical image spectral band information.

S3: establishing a relation between a pure soil surface reflection coefficient and a DDM reflection power;

specifically, according to the water cloud model, the influence of the vegetation reflection coefficient and the attenuation coefficient in step S2 is eliminated from the overall retroreflection coefficient in step S1, and the relationship between the pure soil surface reflection coefficient and the DDM reflected power measurement value is established.

S4: establishing a regression model of the soil humidity and DDM reflected power measurement values;

specifically, according to the relationship between the reflection coefficient of the pure soil surface and the DDM reflected power measurement value in step S3, a model of the influence of the soil humidity on the reflection coefficient of the pure soil surface is introduced, so as to establish a regression model of the soil humidity and the DDM reflected power measurement value. And substituting the measured soil humidity of the control point and the DDM reflected power measured value of the corresponding point, solving unknown parameters in the model by using a least square method, and thus obtaining a regression model for estimating the soil humidity by using the DDM reflected power measured value.

S5: inverting the soil humidity;

specifically, the reflection power of a certain point in the data of the satellite-borne GNSS-R reflectometer in the vegetation coverage area is selected and substituted into the regression model established in the step S4 to reverse the soil humidity.

The method has the characteristics of realizing the fusion and inversion of the optical remote sensing data and the satellite-borne GNSS-R reflectometer data for the soil humidity of the vegetation coverage area, along with high efficiency, convenience, large range and the like.

The invention is further illustrated by the following specific examples in combination with the accompanying drawings. The technical scheme adopted by the invention comprises the following steps:

step 1: constructing a relation between an overall backward reflection coefficient and DDM reflection power by using satellite-borne GNSS-R reflectometer data;

this step can be implemented in the following way:

1) downloading data: from TDS-1 official website: http:// merbys. c o.uk/data-access downloads GNSS-R reflectometer data GN_im。

2) Acquiring reflected power Pr;

TDS-1 data file GN_imSome data information about the trajectory of the GNSS satellites is provided. And importing the data file into a Google map, seeing a plurality of tracks, marking a track serial number T, acquiring a plurality of reflection areas on the track T, acquiring a DDM (data division multiplexing) map on each reflection area, and reading the reflection power Pr of the reflection areas of the vegetation coverage area from the DDM map, wherein the areas discussed below are the vegetation coverage areas, and the appearing parameters are also the parameters of the vegetation coverage areas.

3) Constructing the reflection power Pr and the total retroreflection coefficient delta_totalThe relationship of (1);

i. solving the total back reflection coefficient delta by using the reflection power Pr unit of the microwave radar as db_totalParameter V_in、V_outRespectively microwave voltage at reflection point of microwave incidence mirror and microwave voltage output after reflection by reflection point, and V_inConstant can be set as an unknown constant:

Pr＝10lg(V_out)²(1)

from equations (1) and (2), the regional DDM reflected power Pr and the overall retroreflection coefficient δ can be constructed_totalThe relationship is as follows:

step 2: utilizing the multispectral image data to solve vegetation reflection coefficients and attenuation coefficients of vegetation coverage areas;

this step can be implemented in the following way:

the multispectral image processing is shown in fig. 3.

1) Optical image preprocessing

High-resolution No. 1 spectral image data O of observation area downloaded on land satellite data observation platform_im[j]Wherein j represents each pixel of the spectral image data. After acquiring the spectral image data O_im[j]Then, the preprocessing of the spectral image data by utilizing the ENVI software comprises the following steps: projection conversion, radiometric calibration, atmospheric correction, geometric correction, image fusion and image cutting are carried out to finally obtain the spectral image data OR of the region of interest after pretreatment_im[j]The spectral information of each band of each pixel can be known, for example, the green band b_GREENOr near infrared band b_NIRAnd the like.

2) Differentiating vegetation coverage areas

NDWI value is normalized water index, and is used as important index for distinguishing water body, firstly, the NDWI value is obtained by using formula (4), and the spectral image data OR is removed_im[j]The water body part, equation (4) is as follows:

wherein, b_GREENDenotes the green band, b_NIRAnd representing a near infrared band, setting a threshold value, preferably-0.05, according to the obtained NDWI value, and considering the water body part when the NDWI value is more than-0.05 and considering the water body part as a non-water body when the NDWI value is less than-0.05. Therefore, the water body can be removed through image cutting, and the spectral image data after the water body is removed is set as OW_im[j]。

ii, calculating to obtain OW by using the relation between the NDVI index and the spectral band_im[j]The NDVI index of each pixel point is calculated by the following method:

wherein, b_NIRIn the near infrared band, b_REDIs in the infrared band. In terms of NDVI threshold settings, provision is made for: when the NDVI value is more than 0.4, the vegetation coverage area is determined, when the NDVI value is less than 0.4, the vegetation coverage area is determined to be vegetation low coverage area or bare soil area, and the vegetation coverage area is not considered to be plantedIs affected. Thereby, a vegetation coverage area OH is divided_im[k]And k represents a pixel point of the vegetation coverage area. Therefore, the NDVI value of each vegetation covered pixel is known.

3) Solving vegetation reflection coefficient and attenuation coefficient

According to vegetation coverage area OH_im[k]NDVI value of (D) calculating the water content M of vegetation_vegCalculating the reflection coefficient delta of the vegetation layer by using the water content of the vegetation_veg(theta) and attenuation coefficient tau²(theta) and finally obtaining the soil retroreflection coefficient delta of the pure vegetation covered area_soil(theta). The specific implementation steps are as follows:

i. calculating vegetation water content M by using NDVI_vegNDVI values have been given above:

M_veg＝1.9134NDVI²-0.3215*NDVI (6)

water content M of vegetation_vegCalculating vegetation retroreflection coefficient delta_veg(theta) and attenuation coefficient tau²(θ), the formula is as follows:

δ_veg(θ)＝A*M_vegcos(θ)[1-τ²(θ)](7)

wherein, delta_veg(theta) is the layer retroreflectivity, τ²(theta) is a secondary attenuation factor of radar waves penetrating through a vegetation layer, theta represents an incident angle of the radar sensor and can be obtained from an image file, and M is_vegFor vegetation water content, A, B is an empirical parameter of vegetation type, and the specific value varies according to the actual situation of the research area.

For the values of the empirical parameter A, B, the following values in table 1 may be referred to:

TABLE 1 empirical parameters of vegetation type in Water cloud model

Empirical parameters	Pasture land	Grass land	All vegetation	Winter wheat
					A	0.0009	0.0014	0.0012	0.0018
B	0.032	0.084	0.091	0.138

And step 3: establishing the relation between the reflection coefficient of the pure soil surface and the DDM reflection power

As an alternative embodiment, the following method can be used:

1) the radar retroreflection coefficient of the vegetation high coverage area is composed of the soil reflection coefficient subjected to twice attenuation and the vegetation layer reflection coefficient. Vegetation reflection coefficient delta solved according to step 2_veg(theta) and attenuation coefficient tau²(theta), calculating the retroreflection coefficient of the soil surface layer of the vegetation completely covered area according to the water cloud model, wherein the formula is as follows:

case of vegetation full coverage:

δ_total(θ)＝δ_veg(θ)+τ²(θ)*δ_soil(θ) (9)

wherein, delta_total(theta) reflectometer data GN_imThe overall retroreflection coefficient of the reflection area corresponding to the spectral image on the calculated trajectory T is already found in step 1, δ_soilAnd (theta) is the pure soil surface reflection coefficient.

2) Construction of relation between pure soil surface reflection coefficient and DDM reflection power

Substituting equation (3) with:

the DDM reflected power Pr of each vegetation coverage area on the track T is given in step 1, and the vegetation reflection coefficient delta_veg(theta) and attenuation coefficient tau²(θ) is already known in step 2.

And 4, step 4: establishing regression model of soil humidity and DDM reflected power measurement value

1) According to an empirical model, the soil humidity mv and the soil surface pure reflection coefficient delta_soilThe linear relationship is as follows:

δ_soil＝a₁mv+a₂(12)

2) establishing the relation between the soil humidity and the DDM reflected power measurement value

The united type (11) and (12) establish the relation between the soil humidity of the vegetation coverage area and the measured value of the DDM reflected power:

Pr＝20lgV_in[τ²(a₁mv+a₂)+δ_veg](13)

3) solving for V_in、a₁、a₂To establish a regression model

Actually measuring soil humidity mv of control points of m vegetation coverage areas_iWherein, the value range of i is 1 to m, m is an integer larger than 3, and then the corresponding control point reflection power Pr is read on the GNSS-R reflection data track T in sequence_i. The vegetation reflection coefficient delta of the corresponding vegetation coverage area can be obtained by the step 2_veg,i(theta) and attenuation coefficient

i. Calculate the least squares residual sum of squares:

g (a) then₁，a₂，V_inRespectively calculating partial derivatives and setting them to zero

Wherein m is 1, 2. Solving 3 unknowns a by simultaneous equations₁，a₂，V_in

By the parameter a found_m，V_inDetermining a regression function model of soil moisture and DDM reflected power measurements for a vegetation covered area, wherein the vegetation reflection coefficient δ_veg(theta) and attenuation coefficient tau²(θ) can be obtained from the information of the corresponding region point wave band of the multispectral image, and the detailed explanation of step 2 has been given:

and 5: inversion of soil moisture

As shown in fig. 4, in order to invert the soil moisture of a certain vegetation coverage area k, the specific steps are as follows:

1) firstly, downloading GNSS-R reflectometer data covered with the track of the area, marking the serial number T of the track, reading a DDM (distributed data model) graph of a vegetation coverage area k to be solved on the track, and acquiring the DDM reflection power Pr of the area from the DDM graph_k。

2) And downloading the multispectral image data, and cutting out the vegetation coverage area k. Obtaining NDVI_kCalculating the value of the vegetation reflection coefficient delta according to the step 2 by combining the multispectral wave band information_veg,k(theta) and attenuation coefficient

3) Reflecting power Pr of DDM_kVegetation reflection coefficient delta_veg,k(theta) and attenuation coefficient

And (5) substituting the regression model of the soil humidity and the DDM reflected power measurement value finally given in the step (4), namely the formula (15), so that the soil humidity value of the vegetation coverage area k can be inverted.

It should be understood that the vegetation reflection coefficient when the model of the bare soil surface or vegetation low coverage area is established

Is 0, attenuation coefficientThe value is equal to 1 and the following steps coincide with vegetation coverage.

In another embodiment of the present invention, there is also provided a soil moisture estimation device for a vegetation coverage area, including:

the first relation determining module is used for constructing the relation between the soil total retroreflection coefficient and the DDM reflected power based on the satellite-borne GNSS-R reflectometer data;

the second relation determination module is used for obtaining vegetation reflection coefficients and attenuation coefficients of the vegetation coverage area by using the optical image data;

the third relation determination module is used for eliminating the influence of the vegetation reflection coefficient and the attenuation coefficient based on the relation between the soil total backward reflection coefficient and the DDM reflection power to obtain the relation between the pure soil surface reflection coefficient of the vegetation covered area and the DDM reflection power;

the regression model determining module is used for constructing a relation between a soil humidity value of the vegetation coverage area and the DDM reflected power according to the relation between the pure soil surface reflection coefficient and the DDM reflected power of the vegetation coverage area and the linear relation between the soil humidity and the pure soil surface reflection coefficient, and obtaining an unknown parameter in the relation between the soil humidity value of the vegetation coverage area and the DDM reflected power according to the actually measured soil humidity of a plurality of control points and the DDM reflected power measured value of the corresponding point, so that a regression model for estimating the soil humidity by using the DDM reflected power measured value is obtained;

and the inversion module is used for selecting the reflection power of a target point in the satellite-borne GNSS-R reflectometer data of the vegetation coverage area and inverting the soil humidity of the vegetation coverage area based on the regression model.

The specific implementation of each module may refer to the description of the method embodiment, and the embodiment of the present invention will not be repeated.

The invention is realized mainly in the following two aspects:

(1) the method comprises the following steps of fusing satellite-borne GNSS-R reflectometer data and multispectral information to remove vegetation coverage to influence inversion of soil humidity: in the process of inverting the soil humidity by using the satellite-borne GNSS-R reflection count data, in order to eliminate the influence of vegetation coverage, the vegetation information, such as a vegetation reflection coefficient and an attenuation coefficient of the vegetation layer to the soil reflection coefficient, is obtained by using multispectral information. By means of a water cloud model, the relation between the total soil retroreflection coefficient read by satellite-borne GNSS-R reflectometry data and the reflection power is fused with the vegetation layer reflection coefficient and the vegetation attenuation coefficient, and the functional relation between the pure soil reflectance coefficient under the vegetation coverage condition and the DDM reflection power read by the satellite-borne GNSS-R reflectometry data can be obtained.

(2) Constructing a regression model of soil humidity and DDM reflected power under vegetation coverage: according to the function relation between the pure soil reflection coefficient under the vegetation cover condition in the characteristic 1 and the DDM reflection power read from the satellite-borne GNSS-R reflectometer data, and by combining the linear relation between the soil humidity and the pure soil surface reflection coefficient in the empirical model, a relation model between the soil humidity value and the DDM reflection power under the vegetation cover condition can be constructed. And collecting soil humidity values of actually measured areas of a plurality of control points and DDM (distributed data model) reflected power in GNSS-R reflectometer data of corresponding areas, and establishing a functional relation between the soil humidity and the reflected power under the vegetation coverage condition by using the least square residual sum of regression equation.

It should be understood that parts of the specification not set forth in detail are well within the prior art.

It should be noted that, according to the implementation requirement, each step/component described in the present application can be divided into more steps/components, and two or more steps/components or partial operations of the steps/components can be combined into new steps/components to achieve the purpose of the present invention.

It will be understood by those skilled in the art that the foregoing is only a preferred embodiment of the present invention, and is not intended to limit the invention, and that any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the present invention.