基于Landsat8和Sentinel-1A数据的焦作广利灌区夏玉米土壤墒情监测方法研究

郭二旺，郭乙霏，罗蔚然，王文婷

(1.焦作市水利勘测设计院，河南 焦作 454003；2.东北农业大学水利与土木工程学院，哈尔滨 150003；3.郑州大学水利与环境学院，郑州 450001；4.焦作市抗旱防汛通讯站，河南 焦作 454150)

0 引 言

土壤墒情是水文学、气候学、农业学、生态环境学等各个研究领域中的重要指标，在生态系统水循环中扮演着非常重要的角色，是农作物旱情监测、估产等研究的重要参数[1-4]。土壤墒情在时间、空间范围上变化较大，传统的土壤墒情监测方法是通过设立监测站点进行人工监测，不能满足较大范围上对土壤墒情进行高精度、实时监测的要求[5]。随着微波遥感理论基础经过多年不断的发展，利用微波数据对土壤水分监测的研究已取得了突破性进展，特别是合成孔径雷达(Synthetic Aperture Radar, 简称SAR)技术不断的改进，以其多频率、多极化(或者全极化)、多角度、可变工作模式等能力，在土壤水分研究中变得日益重要[6]。根据微波遥感反演对象不同可分为裸露地表区和植被覆盖区。国内外研究现状将针对国内外裸露地表的反演和植被覆盖区的反演进行陈述。对于裸露地表土壤墒情反演，Fung[7]等构建了IEM 模型(integral equation model)。IEM 模型是基于电磁波辐射传输方程的地表散射模型，能在一个很宽的地表粗糙度范围内再现真实地表后向散射情况，已经被广泛应用于微波地表散射、辐射的模拟和分析。Chen 等[8]对 IEM 模型存在的不足之处进行分析，从而提出AIEM 模型。对于微波遥感反演植被覆盖区的研究中，Ulaby的贡献是毋庸置疑的，其主要突出贡献为反演植被覆盖区域奠定了扎实的基础：一是水云模型[9](Water-Cloud model)，该模型是Ulaby和Attema等以植被覆盖区作为研究对象，通过简化植被层的散射机制，得到的适用于低矮植被覆盖地表的模型；一是密歇根微波植被散射模型(Michigan Microwave Canopy Scattering model, MIMICS)[10]，该模型是 Ulaby等人于 1990 年基于微波辐射传输式子的一阶解建立的描述植被覆盖地区散射特性的模型。

本研究采用Sentinel-1A作为雷达遥感监测的数据源，结合Landsat8数据来进行土壤墒情反演，基于改进的水云模型研究植被的后向散射作用，尝试消除植被效应得到地表后向散射系数，根据地表反射率和地表参数关系实现研究区土壤墒情的反演。

1 研究区概况与数据选取

1.1 研究区简介

(1)地理位置。广利灌区位于河南省焦作市，地理位置为东经112°37′～113°13′，北纬34°55′～35°11′，灌区范围在济源市、沁阳市沁南地区、温县城、武陟沁南几个地区都有覆盖。灌区的南面是黄河，北面挨着沁河，属于黄河流域，引沁河水灌溉。灌区总面积3.4 万hm2，如图1所示。灌区内主要种植小麦、玉米等，夏玉米的播种期为5月，5月中到7月为幼苗期，9月底成熟。

图1 广利灌区总体布置图Fig.1 Distribution map of guangli irrigation area

(2)水文特征。规划面积34 000 hm2，其中灌溉面积20 666 hm2、补源面积13 333 hm2，设计流量28.0 m3/s。现有总干渠1条，长29.1 km；灌溉干支渠34条，长214.25 km；补水干支渠19条，长134.06 km。灌区多年平均降水量590.66 mm，降水量年内分配极不均匀，多集中在汛期夏秋两季之间，其中汛期6-9月份降水量占全年降水量的70%以上。由于受到降雨的影响，地下水位从4月中下旬开始由3 m以下逐渐上升，到9月中旬地下水位逐渐降至3 m以下。年内蒸发量的变化基本与气温变化一致，变化规律呈正态分布，一年中以6月份最大，平均306.2 mm，11月份最小，平均86.7 mm。

1.2 遥感数据选取

本文采用的是2017年7月的Sentinel-1A数据和Landsat8数据，用SNAP软件对Sentinel-1A影像的VV和VH影像进行辐射定标，将影像上的像元亮度值转化为后向散射系数；对得到的结果去噪声处理后，做投影转换和坐标系定义，将该影像设置为和Landsat8影像相同的投影坐标系，即UTM投影、WGS-84坐标系。研究中所用到的是Landsat8的OLI传感器影像，对分辨率30 m的Landsat8数据利用广利灌区的矢量边界裁剪得到研究范围的影像数据；利用ENVI5.2软件对得到的研究区范围的Landsat8影像进行辐射定标和Flaash大气校正，并利用遥感影像计算NDVI，根据阈值分割法提取夏玉米分布，结果如图2所示。

图2 夏玉米空间分布图Fig.2 Spatial distribution map of summer corn

1.3 实测野外样品数据

选择卫星过境的同一时间进行实地采样，使用GPS记录样点坐标。本研究选取2017 年7 月15 日的位于广利灌区的24个实测点的土壤平均含水量数据、植被含水量数据及地表粗糙度等数据，图3所示。①土壤墒情数据 ：在焦作市广利灌区进行实测土壤墒情时，采用方法为土钻取土，然后用烘箱进行烘干。广利灌区的土壤为粉质壤土，土质分布均匀，容重为1.40 g/m3，田间持水率为28%(重量含水率)。在野外采样的过程中，每个采样点取深度为0～10 cm处的土壤，用铝盒收集样品并及时记录总重量，对样本进行烘干处理获得烘干后的样本铝盒的重量并记录，计算土壤含水量数据。图3是研究区野外采样点的分布情况。②植被含水量数据：测量采样点植被的高度，叶片长度、宽度，分别记录采样点植株叶和根烘干前和烘干后的重量，计算鲜重含水量和干重含水量。

图3 研究区样点分布Fig.3 Sample point distribution in research area

(3)地表粗糙度:本研究中的粗糙度数据是采用剖面板测量法测量采样点的数据，随后利用离散数据的均方差高度以及表面相关长度计算公式来计算S(均方根高度)和CL(相关长度)，对于一维离散数据，表面高度的标准偏差(均方根高度)s为：

在针对城乡进行规划和具体设计的时候，为了保证城乡规划的科学性和合理性，首先要做的一点就是要树立具有生态文明特征的规划理念，并且在实践中可以将该理念落实到实处。也就是说在实践中，要采取有针对性的措施，尽可能保证遵循人与自然环境相互协调的基础发展原则。首先，要采取有针对性的措施，尽可能降低CO2的排放力度，同时还要使用一些具有清洁性功能特征的能源。这样不仅能够针对城乡规划相关措施具体落实过程中的一些污染问题进行有针对性的处理，而且还能够保证生态文明城市理念在实践中的有效落实。

式中：N为采样数目。

对于离散数据，相距x′=(j-1)x(j为自然数)的归一化自相关函数由下式给出：

当相关函数ρ(x′)=1/e时,间隔x′称为表面相关长度。

为了在一定程度上避免采样误差，在进行样本采集时，对于每一个样点附近分3个子样点，将子样本点的均值作为该点土壤含水量，同时考虑遥感影像的与实测数据之间的差异性,将测点对应的影像点的邻近像元值进行均值化处理，使得采样点与影像点能够更好的匹配。

2 植被覆盖地表土壤墒情反演模型

微波传感器接收的回波信号中包含了土壤水分和植被水分的信息，地表植被对微波遥感反演土壤水分的精度有较大的影响。因此在进行土壤墒情反演时，消除植被的影响是研究的重点。

2.1 植被散射模型

Attema和Ulaby等在1987年提出了以农作物为研究对象的水云模型。在水云模型中，研究区植被层被理想化为云层一样均匀分布在地表上方，只考虑植被散射和地表散射而忽略了植被和土壤表面的多次散射。可用下式表示：

γ2(θ)=exp[-2Bmv/cos(θ)]

式中：γ2(θ)为植被双层衰减因子；θ指入射角；A和B的值取决植被类型及入射电磁波频率,A=0.001 2，B=0.091；mv表示植被含水量。

本文将采取归一化水指数(NDWI)来估算植被含水量，Jackson(Jacson et al.，2004)[11]等人的研究证明了NDWI估算植被指数比NDVI有优势。 公式如下：

从实测的植被含水量数据中选取一部分来建立NDWI和植被含水量的关系，如图4所示：

y=-49.027x+4.413 1

(5)

利用公式(5)经波段运算得到的研究区植被含水量分布。根据实测植被含水量数据检测NDWI估算的VWC值的准确度，结果表明模拟植被含水量和实测值R2=0.535 7，可以用来估算植被含水量。

图4 植被含水量和NDWI之间的关系图Fig4.Map of the relationship between vegetation water content and NDWI

根据公式(1)～(3)运用ENVI5.2中的波段运算结合Sentinel-1A影像和Landsat8影像可以计算去除植被影响的后向散射。利用AIEM模型，研究在VV极化下地表土壤墒情对后向散射的影响。设置AIEM模型中的参数：雷达入射频率采用Sentinel1A固有频率f=5.4 GHz，VV极化，相关长度cl取值范围5～20 cm，均方根高度s取值范围0.3～3.3 cm，由于研究区地势比较平坦，雷达入射角波动范围较小，整体入射角在θ=40°左右变化，因此在参数模拟过程中，入射角固定为40°。模拟结果如图5所示。

图5 不同地表粗糙度下土壤含水量对后向散射系数的影响Fig.5 Effect of soil water content on backscattering coefficient under different surface roughness

2.2 联合地表参数反演土壤墒情

利用AIEM模型模拟的过程中需要输入均方根高度、相关长度、入射角、介电常数等参数，采用Zribi[12]等提出的组合参数(Zs=s2/cl)来表征地表粗糙度参数。其中，用土壤墒情Mv和地表粗糙度Zs的函数来表示同极化雷达后向散射系数[13]：

联合以上公式，组合粗糙度的反演公式为：

结合土壤水分和土壤粗糙度参数的式子：

得到结果：

图6 实测值和模拟值关系曲线Fig.6 Curve of relationship between measured and simulated values

3 反演结果与精度分析

3.1 土壤水分反演的结果

利用回归公式对交叉极化组合反演的结果进行计算，得到VV-VH极化下的土壤墒情空间分布图，如图7所示。在图中可以看到土壤水分大多集中在20%～30%之间，并且土壤水分在空间上的分布也比较均匀。

图7 交叉极化VV-VH的土壤水分反演图Fig.7 Soil moisture inversion map of cross-polarized VV-VH

3.2 精度分析

对模拟得到的土壤水分数据选取9个野外实测数据进行精度验证，图8是交叉极化VV-VH组合反演的土壤水分值与野外实测值的拟合关系，相关系数为0.7。

图8 模拟数据与实测数据的相关关系Fig.8 The correlation between simulated data and measured data

通过模拟数据与实测数据的相关关系可以得知本研究中使用的基于水云模型、AIEM模型和组合粗糙度参数建立的反演模型在广利灌区有一定的适用性。但由于实测数据数量的限制和模型本身存在的问题使得精度不是很高，得到的土壤墒情分布图有一定的参考作用。

4 结 语

本文研究区是河南省焦作市的广利灌区，该研究区的植被类型主要是玉米、大豆、小麦等低矮植被。使用的数据是Sentinel-1A和Landsat8遥感数据，结合光学遥感和微波遥感的特性，考虑到灌区植被的生长和分布状况，选用水云模型去除地表植被对后向散射的影响。研究利用水云模型得到地表后向散射系数之后，通过VV、VH交叉极化组合方式得到交叉极化差，结合组合粗糙度参数建立研究区的土壤水分反演模型，结果较为满意。本文研究了植被体散射之间单次散射关系，并未考虑土壤--植被层之间的多次交互散射的影响。为进一步提高土壤水分反演精度，应考虑实际植被枝叶形态、朝向、体积以及介电常数等因素，进而来获取更高的反演精度。