一种植被覆盖区域的土壤湿度估计方法及装置

摘要

本发明公开了一种植被覆盖区域的土壤湿度估计方法及装置，属于电子、信息及光学遥感等领域，首先，借助星载GNSS‑R反射计数据构建总体后向反射系数和DDM反射功率关系；利用光学影像数据求解植被覆盖区域植被反射系数和衰减系数；融合星载GNSS‑R反射计数据总体后向反射系数和DDM反射功率的关系以及光学数据中植被信息，求解植被覆盖区域纯粹土壤表面反射系数和DDM反射功率关系；利用多个控制点土壤湿度的实测数据建立土壤湿度和DDM反射功率回归模型；以某点GNSS‑R数据DDM反射功率代入回归方程反演土壤湿度。本发明的实施是一次星载GNSS‑R反射计数据与多光谱遥感数据影像共同反演植被覆盖区域土壤湿度的有效融合。

说明书

技术领域

本发明属于电子、信息及遥感技术领域，更具体地，涉及一种基于星载GNSS-R反射计与多光谱数据融合的土壤湿度估计方法及装置。

背景技术

土壤湿度是一种全球水循环的重要状态参数，在农业灌溉、环境研究、自然灾害预警及水土治理等方面都起着至关重要的作用。

目前测量土壤湿度的方法存在以下几种：传统的土壤水分的测量方法包括重量法、时域反射法、电容传感器、中子仪、电阻率仪、热脉冲传感器和光纤传感器等方法，这种实测方法优势在于测量精度大，并且可以除去植被影响。但是不足之处是定点测量，测量范围小，人工成本较高。

地基测量方面包括最新发出的基于GPS多星三频数据融合的GNSS-IR土壤湿度反演方法，利用不同轨道、不同频率的差异性和互补性，提出了一种将GPS多星的L1，L2，L5频段数据加权融合进行联合反演的方法，虽然提高了测量精度，但是具有测量范围局限性的弊病，测量范围小，并且只考虑了裸土情况下土壤湿度的测量，没有加入植被考虑因素。

在利用卫星遥感技术监测土壤湿度方面，主要方法有光学、热红外和微波遥感法等，而光学遥感法主要是利用光谱反射特性来估算土壤含水量，但容易受到天气影响。热红外遥感法是利用土壤的热特性来反演土壤湿度，但是在植被覆盖度高的地方，由于植被掩盖了土壤信息，会影响估算土壤含水量的准确性，因此该方法仅仅适用于裸土和植被稀疏地区土壤含水量的监测。利用合成孔径雷达SAR卫星的遥感土壤湿度法主要手段是通过土壤表面反射信号与发射信号功率比值反演出土壤水分。虽然解决了植被覆盖下土壤湿度的测量问题，但是该方法受限于SAR卫星的重返周期，无法保证数据连续性，同一地方需几日才能获得一场数据。同样，在专利CN103940834A中公开了一种采用合成孔径雷达技术测量土壤湿度的方法，同样也需要利用两个不同时刻t1和t2重复观测获取的观测区域的极化合成孔径雷达复图像来反演土壤湿度。

目前，利用导航信号的反射信号进行地面遥感是一个研究热点，当反射信号接收机安装在卫星上时，称为星载GNSS-R反射计。星载GNSS-R反射计利用的导航信号处于L波段，是一种微波信号，具有穿透性，可以不受云层，雨雪等天气条件的影响，同时还可以穿透植被，因此可以作为反演植被覆盖下土壤湿度的有利手段。但是植被参数很多情况下是未知的，导致植被对土壤湿度反演的结果造成干扰。

综上，亟需一种有效的植被覆盖区域土壤湿度估计方法。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提出了一种植被覆盖区域的土壤湿度估计方法及装置，由此解决现有植被覆盖区域土壤湿度估计方法不能有效的进行土壤湿度估计的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种植被覆盖区域的土壤湿度估计方法，包括：

(1)基于星载GNSS-R反射计数据构建土壤总体后向反射系数和DDM反射功率之间的关系；

(2)利用光学影像数据得到植被覆盖区域的植被反射系数和衰减系数；

(3)基于所述土壤总体后向反射系数和DDM反射功率之间的关系，消除所述植被反射系数和衰减系数的影响，得到植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系；

(4)根据所述植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系，结合土壤湿度与纯粹土壤表面反射系数的线性关系，构建植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系，再根据实测的若干控制点土壤湿度和相应点的DDM反射功率测量值，得出所述植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系中的未知参数，从而得到利用DDM反射功率测量值估算土壤湿度的回归模型；

(5)选取植被覆盖区域星载GNSS-R反射计数据中目标点的反射功率，基于所述回归模型，反演出植被覆盖区域的土壤湿度。

优选地，由构建土壤总体后向反射系数和DDM反射功率之间的关系，其中，δ_total表示土壤总体后向反射系数，Pr表示DDM反射功率，V_in表示微波入射镜反射点微波电压。

优选地，步骤(2)包括：

(2.1)对光学影像数据进行预处理，得到感兴趣区域的光谱影像数据；

(2.2)获取所述感兴趣区域的光谱影像数据中每个像素点的归一化水指数，根据各归一化水指数与第一预设阈值之间的关系，去除所述感兴趣区域的光谱影像数据中的水体部分，得到去除水体之后的光谱影像数据；

(2.3)获取所述去除水体之后的光谱影像数据中每个像素点的归一化植被指数，根据各归一化植被指数与第二预设阈值之间的关系确定植被覆盖区域；

(2.4)根据所述植被覆盖区域中每个像素点的归一化植被指数得到所述植被覆盖区域的植被含水量，进而根据所述植被含水量得到植被覆盖区域的植被反射系数和衰减系数。

优选地，由δ_veg(θ)＝A*M_vegcos(θ)[1-τ²(θ)]确定植被覆盖区域的植被反射系数，由确定植被覆盖区域的衰减系数，其中，δ_veg(θ)为植被反射系数，τ²(θ)为雷达波穿透植被层的衰减系数，θ表示雷达传感器的入射角，M_veg为植被含水量，A、B为植被类型经验参数。

优选地，由确定植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系，其中，δ_soil(θ)为纯粹土壤表面反射系数。

优选地，步骤(4)包括：

(4.1)由δ_soil(θ)＝a₁mv+a₂确定土壤湿度mv和纯粹土壤表面反射系数δ_soil(θ)之间的线性关系，其中，a₁、a₂为系数；

(4.2)由Pr＝20lgV_in[τ²(θ)(a₁mv+a₂)+δ_veg(θ)]确定植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系；

(4.3)实测m个植被覆盖区域控制点土壤湿度mv_i，再依次在GNSS-R反射计数据轨迹上读取相应控制点的反射功率Pr_i，然后得到相应植被覆盖区域植被反射系数δ_veg,i(θ)和衰减系数其中，i的取值范围为1到m；

(4.4)根据mv_i、Pr_i、δ_veg,i(θ)及采用最小二乘法得到所述植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系中的参数a₁、a₂和V_in，进而得到利用DDM反射功率测量值估算土壤湿度的回归模型。

优选地，由得到利用DDM反射功率测量值估算土壤湿度的回归模型。

优选地，步骤(5)包括：

(5.1)选取植被覆盖区域星载GNSS-R反射计数据中目标点的反射功率；

(5.2)利用光学影像数据得到植被覆盖区域中所述目标点的归一化植被指数，进而根据所述目标点的归一化植被指数得到所述目标点的植被反射系数和衰减系数；

(5.3)根据所述目标点的反射功率、所述目标点的植被反射系数和衰减系数，利用所述回归模型得到植被覆盖区域的土壤湿度。

按照本发明的另一个方面，提供了一种植被覆盖区域的土壤湿度估计装置，包括：

第一关系确定模块，用于基于星载GNSS-R反射计数据构建土壤总体后向反射系数和DDM反射功率之间的关系；

第二关系确定模块，用于利用光学影像数据得到植被覆盖区域的植被反射系数和衰减系数；

第三关系确定模块，用于基于所述土壤总体后向反射系数和DDM反射功率之间的关系，消除所述植被反射系数和衰减系数的影响，得到植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系；

回归模型确定模块，用于根据所述植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系，结合土壤湿度与纯粹土壤表面反射系数的线性关系，构建植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系，再根据实测的若干控制点土壤湿度和相应点的DDM反射功率测量值，得出所述植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系中的未知参数，从而得到利用DDM反射功率测量值估算土壤湿度的回归模型；

反演模块，用于选取植被覆盖区域星载GNSS-R反射计数据中目标点的反射功率，基于所述回归模型，反演出植被覆盖区域的土壤湿度。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

本发明提出了利用星载GNSS-R反射计数据和光学遥感影像融合的方式来反演植被覆盖区域土壤湿度。由于传统的SAR方法具有重返周期长等不足，而本发明的方法发挥了导航卫星数量多且全球覆盖的优势，同时针对植被遮盖的情况，利用多光谱影像信息消除植被影响，从而反演出区域土壤湿度。本发明因其范围覆盖广、高效的特点，可以应用于农业等其他领域，本发明的实施是一次星载GNSS-R反射计数据与多光谱遥感数据影像共同反演植被覆盖区域土壤湿度的有效融合。

附图说明

图1是本发明实施例提供的一种模型图；

图2是本发明实施例提供的一种方法流程示意图；

图3是本发明实施例提供的一种利用多光谱影像数据求解植被覆盖区域植被反射系数和衰减系数流程图；

图4是本发明实施例提供的一种利用星载GNSS-R反射计数据反演植被覆盖区域土壤湿度示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明充分利用星载GNSS-R反射计数据L波段微波信号具有穿透性、可以不受云层、雨雪等天气条件的影响，同时还可以穿透植被的特性和多光谱遥感数据波段信息求解植被参数的方法，再从反射计数据中扣除植被衰减，融合二者共同反演植被覆盖下土壤湿度。

本发明提出将星载GNSS-R反射计数据与多光谱高分数据影像结合，一方面利用星载GNSS-R反射计数据，构建总体后向反射系数和反射功率的计算模型，另一方面利用光学遥感图像中的光谱信息，区分出植被覆盖区域，计算植被反射系数和衰减系数，借助水云模型得到植被下纯粹土壤表面反射系数，最后利用土壤湿度对纯粹土壤表面反射系数的影响模型，构建植被覆盖下纯粹土壤湿度与反射功率的回归模型，并反演土壤湿度，本发明模型图如图1所示。

如图2所示，是本发明实施例提供的一种方法流程示意图，在图2所示的方法中包括以下步骤：

S1：利用星载GNSS-R反射计数据构建总体后向反射系数和DDM反射功率关系；

具体地，由星载GNSS-R反射计的延迟多普勒图(Delay-Doppler Map)数据得到DDM反射功率，设定直达电压为未知数，从而得到由DDM反射功率测量值计算总体后向反射系数的公式。

S2：利用多光谱影像数据求解植被覆盖区域植被反射系数和衰减系数；

具体地，先对光学影像进行预处理求得归一化植被指数，然后利用该植被指数区分出植被覆盖区域。借助光学影像光谱波段信息求解植被反射系数和衰减系数。

S3：建立纯粹土壤表面反射系数和DDM反射功率关系；

具体地，依据水云模型，从步骤S1中的总体后向反射系数中消除步骤S2中植被反射系数和衰减系数的影响，建立纯粹土壤表面反射系数和DDM反射功率测量值的关系。

S4：建立土壤湿度和DDM反射功率测量值的回归模型；

具体地，根据步骤S3中纯粹土壤表面反射系数和DDM反射功率测量值的关系，再带入土壤湿度对纯粹土壤表面反射系数的影响模型，从而建立土壤湿度和DDM反射功率测量值的回归模型。再代入实测控制点土壤湿度和相应点的DDM反射功率测量值，利用最小二乘法，解出模型中的未知参数，从而得到利用DDM反射功率测量值估算土壤湿度的回归模型。

S5：反演土壤湿度；

具体地，选取植被覆盖区域星载GNSS-R反射计数据里某一点的反射功率，代入步骤S4建立的回归模型，反演出土壤湿度。

本发明的特点是实现光学遥感数据与星载GNSS-R反射计数据的融合反演植被覆盖区域的土壤湿度，具有高效、便捷、大范围等特点。

下面以具体实施例，并结合附图对本发明作进一步说明。本发明所采用的技术方案包括以下步骤：

步骤1：利用星载GNSS-R反射计数据构建总体后向反射系数和DDM反射功率关系；

本步骤可以采用以下方式实现：

1)下载数据：从TDS-1官网：http://merrbys.c o.uk/data-access下载GNSS-R反射计数据GN_im。

2)获取反射功率Pr；

TDS-1数据文件GN_im会提供关于GNSS卫星运行轨迹的一些数据信息。将数据文件导入谷歌地图，可以看到多条轨迹，标记一条轨迹序号T，获取轨迹T上多个反射区域，得到每个反射区域上的DDM图，从DDM图中可以读取植被覆盖区反射区域的反射功率Pr，以下所讨论区域均为植被覆盖区域，所出现参数也为植被覆盖区域参数。

3)构建反射功率Pr与总体后向反射系数δ_total的关系；

i.由微波雷达反射功率Pr单位为db，求解总体后向反射系数δ_total，参数V_in、V_out分别为微波入射镜反射点微波电压与经反射点反射后输出微波电压，且V_in恒定可设为未知常数：

Pr＝10lg(V_out)²(1)

ii.由公式(1)和公式(2)，可以构建区域DDM反射功率Pr与总体后向反射系数δ_total关系：

步骤2：利用多光谱影像数据求解植被覆盖区域植被反射系数和衰减系数；

本步骤可以采用以下方式实现：

多光谱影像处理过程如图3所示。

1)光学影像预处理

在陆地卫星数据观测平台下载观测区域的高分1号光谱影像数据O_im[j]，其中，j代表光谱影像数据的每个像素点。在获取到光谱影像数据O_im[j]之后，利用ENVI软件对光谱影像数据的预处理包括以下几个步骤：投影转换、辐射定标、大气校正、几何校正、图像融合及图像裁剪，最终得到预处理之后感兴趣区域的光谱影像数据OR_im[j]，每个像素点的各个波段的光谱信息均可以得知，例如绿波段b_GREEN或近红外波段b_NIR等。

2)区分出植被覆盖区域

i.NDWI值是归一化水指数，作为区分水体的重要指数，先利用公式(4)求取NDWI的值，去除光谱影像数据OR_im[j]水体部分，公式(4)如下：

其中，b_GREEN表示绿波段，b_NIR表示近红外波段，根据求得的NDWI的值，设定阈值，优选为-0.05，当NDWI值大于-0.05时认为是水体部分，小于-0.05时为非水体。由此经图像裁剪，可以去除水体，设去除水体之后光谱影像数据为OW_im[j]。

ii.借用NDVI指数与光谱波段的关系，计算得到OW_im[j]中每个像素点的NDVI指数，计算方法如下：

其中，b_NIR为近红外波段，b_RED为红外波段。依据NDVI阈值设定，规定：当NDVI值大于0.4时为植被覆盖区域，当NDVI值小于0.4时为植被低覆盖区域或裸土区域，不予考虑植被影响。由此，划分出了植被覆盖区域OH_im[k]，k代表植被覆盖区域像素点。因此每个植被覆盖像素点NDVI值已知。

3)求解植被反射系数和衰减系数

根据植被覆盖区域OH_im[k]的NDVI值计算植被含水量M_veg，利用植被含水量计算出植被层反射系数δ_veg(θ)以及衰减系数τ²(θ)，最后得到纯粹的植被覆盖区域土壤后向反射系数δ_soil(θ)。具体实施步骤如下：

i.利用NDVI计算植被含水量M_veg，NDVI值在上面已经给出：

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

ii.由植被含水量M_veg计算植被后向反射系数δ_veg(θ)和衰减系数τ²(θ)，公式如下：

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

其中，δ_veg(θ)为植被层后向反射系数，τ²(θ)为雷达波穿透植被层的二次衰减因子，θ代表的是雷达传感器的入射角，可以在影像文件中获取，M_veg为植被含水量，A、B为植被类型经验参数，具体的取值根据研究区实际情况的不同而有所差异。

其中，对于经验参数A、B的取值，可以参考如下表1中的取值：

表1水云模型中的植被类型经验参数

经验参数放牧地草地所有植被冬小麦A0.00090.00140.00120.0018B0.0320.0840.0910.138

步骤3：建立纯粹土壤表面反射系数和DDM反射功率关系

作为一种可选的实施方式，可以通过以下方式实现：

1)植被高覆盖区域雷达后向反射系数由土壤反射系数经过两次衰减和植被层反射系数组成。根据步骤2求解的植被反射系数δ_veg(θ)和衰减系数τ²(θ)，依据水云模型，计算植被完全覆盖区域土壤表层的后向反射系数，公式如下：

植被完全覆盖下的情况：

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

其中，δ_total(θ)为反射计数据GN_im计算得到的轨迹T上与光谱影像相对应反射区域的总体后向反射系数，在步骤1中已经求得，δ_soil(θ)为纯粹土壤表面反射系数。

2)构建纯粹土壤表面反射系数和DDM反射功率的关系

代入公式(3)，则有：

其中，轨迹T上各植被覆盖区域的DDM反射功率Pr在步骤1中已经给出，植被反射系数δ_veg(θ)和衰减系数τ²(θ)在步骤2中已经得知。

步骤4：建立土壤湿度和DDM反射功率测量值的回归模型

1)有经验模型可知，土壤湿度mv和土壤表面纯粹反射系数δ_soil有线性关系如下：

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

2)建立土壤湿度和DDM反射功率测量值关系

联立式(11)和(12)，建立了植被覆盖区域土壤湿度和DDM反射功率测量值的关系：

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

3)求解V_in、a₁、a₂的值，建立回归模型

实测m个植被覆盖区域控制点土壤湿度mv_i，其中，i的取值范围为1到m，m为大于3的整数，再依次在GNSS-R反射计数据轨迹T上读取相应控制点反射功率Pr_i。由步骤2可以得到相应植被覆盖区域植被反射系数δ_veg,i(θ)和衰减系数

i.计算最小二乘残差平方和：

然后G(a)对参数a₁，a₂，V_in分别求偏导数并置为零

其中，m＝1,2。联立方程求出3个未知数a₁，a₂，V_in

ii.通过求出的参数a_m，V_in确定植被覆盖区域土壤湿度和DDM反射功率测量值的回归函数模型，其中植被反射系数δ_veg(θ)和衰减系数τ²(θ)可由多光谱影像相应区域点波段信息求得，步骤2已经给出详细解释：

步骤5：反演土壤湿度

如图4，为反演某一植被覆盖区域k的土壤湿度，具体步骤如下：

1)首先下载覆盖有此区域轨迹的GNSS-R反射计数据，标记此条轨迹序号T，读取轨迹上待求的植被覆盖区域k的DDM图，从DDM图上获取该区域的DDM反射功率Pr_k。

2)下载多光谱影像数据，裁剪出此植被覆盖区域k。得到NDVI_k值，结合多光谱波段信息按照步骤2计算得到植被反射系数δ_veg,k(θ)和衰减系数

3)将DDM反射功率Pr_k、植被反射系数δ_veg,k(θ)和衰减系数代入步骤4中最后给出的土壤湿度和DDM反射功率测量值的回归模型，也就是式子(15)，即可反演出此植被覆盖区域k的土壤湿度值。

应当理解的是，裸土表面或者植被低覆盖区域模型建立时植被反射系数为0，衰减系数值等于1，后面步骤与植被覆盖区域一致。

在本发明的另一实施例中，还提供了一种植被覆盖区域的土壤湿度估计装置，包括：

第一关系确定模块，用于基于星载GNSS-R反射计数据构建土壤总体后向反射系数和DDM反射功率之间的关系；

第二关系确定模块，用于利用光学影像数据得到植被覆盖区域的植被反射系数和衰减系数；

第三关系确定模块，用于基于所述土壤总体后向反射系数和DDM反射功率之间的关系，消除所述植被反射系数和衰减系数的影响，得到植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系；

回归模型确定模块，用于根据所述植被覆盖区域纯粹土壤表面反射系数与所述DDM反射功率之间的关系，结合土壤湿度与纯粹土壤表面反射系数的线性关系，构建植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系，再根据实测的若干控制点土壤湿度和相应点的DDM反射功率测量值，得出所述植被覆盖区域土壤湿度值与所述DDM反射功率之间的关系中的未知参数，从而得到利用DDM反射功率测量值估算土壤湿度的回归模型；

反演模块，用于选取植被覆盖区域星载GNSS-R反射计数据中目标点的反射功率，基于所述回归模型，反演出植被覆盖区域的土壤湿度。

其中，各模块的具体实施方式可以参考方法实施例的描述，本发明实施例将不再复述。

本发明的实现主要体现在以下两方面：

(1)融合星载GNSS-R反射计数据和多光谱信息去除植被覆盖层影响反演土壤湿度：利用星载GNSS-R反射计数据反演土壤湿度的过程中，为消除植被覆盖层影响，利用多光谱信息获取植被信息，例如植被反射系数和植被层对土壤反射系数的衰减系数。借助水云模型，用星载GNSS-R反射计数据读取的土壤总体后向反射系数与反射功率的关系融合植被层反射系数和植被衰减系数，可以得到植被覆盖情况下纯粹土壤反射系数与星载GNSS-R反射计数据中读取的DDM反射功率的函数关系。

(2)构建植被覆盖下土壤湿度和DDM反射功率的回归模型：根据特征1中的植被覆盖情况下纯粹土壤反射系数与星载GNSS-R反射计数据中读取的DDM反射功率的函数关系，结合经验模型中土壤湿度与纯粹土壤表面反射系数的线性关系，可以构建植被覆盖情况下土壤湿度值与DDM反射功率的关系模型。收集多个控制点实测区域的土壤湿度值与相对应区域的GNSS-R反射计数据中的DDM反射功率，利用回归方程最小二乘残差平方和建立植被覆盖情况下土壤湿度与反射功率的函数关系。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

需要指出，根据实施的需要，可将本申请中描述的各个步骤/部件拆分为更多步骤/部件，也可将两个或多个步骤/部件或者步骤/部件的部分操作组合成新的步骤/部件，以实现本发明的目的。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。