CN108717044B - 一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法 - Google Patents
一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法 Download PDFInfo
- Publication number
- CN108717044B CN108717044B CN201810505327.2A CN201810505327A CN108717044B CN 108717044 B CN108717044 B CN 108717044B CN 201810505327 A CN201810505327 A CN 201810505327A CN 108717044 B CN108717044 B CN 108717044B
- Authority
- CN
- China
- Prior art keywords
- albedo
- soil
- water content
- vegetation
- vegetation coverage
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N2021/1765—Method using an image detector and processing of image signal
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N2021/178—Methods for obtaining spatial resolution of the property being measured
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N2021/1793—Remote sensing
- G01N2021/1797—Remote sensing in landscape, e.g. crops
Abstract
本发明提供了一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法，本方法的步骤包括（1）材料及数据；（2）辅助参数计算（地表反照率、植被覆盖度计算）；（3）从地表反照率中分解出裸土反照率；（4）土壤含水量遥感估算模型与制图；（5）估算方法准确性评价。该方法直接运用裸土反照率来计算土壤含水量，遥感物理意义明确，凸显了卫星遥感估算土壤含水量的实质内涵；克服了传统方法（基于地表反照率的方法，如裴浩等.1999）需要计算测量中间参数的缺陷，避免了误差在参数间的传递，没有多余的中间环节，简单易行；突破了传统方法（如光学植被盖度方法）必须具备的区域条件限制。该方法实现了可靠、稳定和高精度地估算任意较大区域的表层土壤含水量。
Description
技术领域
本发明涉及卫星遥感的应用技术，具体地说是一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法。
背景技术
土壤含水量，是植物生长发育的基础条件，也是农牧业干旱监测的有效因子，对土壤的各种理化性质具有深刻影响。它参与生物圈、大气圈的物质能量循环，间接地决定与改变了陆地表面的太阳短波辐射和地表长波辐射吸收、反射和发射，进而影响着地球系统能量收支平衡。因此，土壤含水量在农业、陆地生态系统、气候等领域研究中扮演着重要的角色，其含量的准确，快速测定具有重要的作用和意义。
传统的土壤含水量测量方法虽然具有很高的精度，但费时、费力，样点有限，获取成本高，难以广泛应用到区域尺度。卫星遥感技术能够在短时间内，低成本地获取广阔地表的地物光谱信息，通过遥感模型，借助尺度推绎方法，实现像元尺度上的土壤含水量测定与空间分布格局，已成为区域尺度土壤含水量估算的主要技术手段，得到了广泛应用。
现有技术存在的问题：卫星遥感技术估算土壤含水量的基础是遥感像元。依据遥感技术成像原理，受空间分辨率影响，卫星遥感图像的像元是混合像元（即对于非人工干扰的自然景观，其遥感像元包含植被与土壤的反射率，也就是说，当有植被覆盖时，遥感像元的反射率不完全是纯土壤（裸土）的光谱反射率，而是土壤与植被反射率的加权和）。因此，植被（植被覆盖度）对土壤反射率的影响巨大，欲进一步提高遥感反演的精度，需将植被对混合像元反射率的贡献设法予以消除或尽可能地降低，“去除”植被光谱的影响，分离出裸土的光谱反射率。
目前，有三大类“去除”植被光谱，分离出裸土的光谱反射率的方法。
一是混合像元分解模型。主要有线性、非线性、几何光学、随机几何和模糊分析5种模型。线性模型的关键是要输入各种地物的参照光谱值，而实际地物的光谱值很难获得，应用困难。非线性模型，计算过程较复杂，其混合像元的分解效果通常好于线性模型，但受残差的影响，其结果仍不理想。几何（光学与随机几何）模型，需要大量参数，如植被的高度、形状参数以及太阳参数（入射与观测方向）等，这些参数通常难以全部获取。模糊分析模型，运用该方法的前提是数据必须符合或接近正态分布，但实际上，并不是所有区域的遥感数据满足这一前提条件。因此，混合像元分解模型虽然取得了较好的效果，但仍然存在着不足。
二是光学植被盖度方法。该方法成功地将混合像元的光谱信息纯化为裸土的光谱信息，但需要用统计方法确定裸土壤的虚生物量本底、全是植被光学信息时的近红外波段的光谱亮度等众多模型参数。这意味着区域内必须有全裸地（无植被）和植被全覆盖地区（植被覆盖度C达到100%），实际上，并非所有区域满足这一条件，因此应用受到限制。
三是梯形或三角形方法。研究表明，地表温度-植被指数、反照率-植被指数、反照率-植被盖度等的散点图能够形成梯形或三角形特征空间，因而得名。因其操作简便，与其他模型结合，可以对植被蒸腾、地表通量和土壤蒸发等进行估算，但对土壤含水量的估算较少。
为此，本发明申请旨在提供一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法。该发明申请从卫星遥感获取的地表反照率（属于混合像元）中提取裸露土壤的反照率，从而剔除了植被覆盖对土壤光谱的影响，恢复土壤本身的光谱信息。运用裸土反照率，能够显著地提高了区域尺度的土壤含水量估算与制图精度。该发明能够满足区域干旱监测、生态农业、精准农业，陆地生态系统和气候变化等领域的土壤含水量科学研究及实际应用需求。
发明内容
鉴于现有技术的不足，本发明提供了一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法，利用能够表达土壤含水量信息的地表反照率-植被覆盖度梯形特征空间，提出了土壤等湿度线及其斜率，基于混合像元的组份及其数学模型，通过严格的数学推导，从地表反照率混合像元中分解出遥感物理意义明确的裸土反照率，进而去除了植被覆盖的影响，实现了可靠、稳定和高精度地估算区域尺度的表层土壤含水量。
为了达到上述目的，本发明采用了如下的技术方案：
本发明提供了一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法，本方法的步骤包括（1）材料及数据；（2）辅助参数计算（地表反照率、植被覆盖度计算）；（3）从地表反照率中分解出裸土反照率；（4）土壤含水量遥感估算模型与制图；（5）估算方法准确性评价。
（1）材料及数据
（a）土壤含水量实测数据
2016年9月27日至10月4日进行野外土壤样点的采集。采用环刀法取地表 (0~10cm、10~20 cm、20~30 cm)的土样。样品带回实验室后，按照下列步骤测量：
准备工作：在室内将铝盒编号并称重，重量记为W0。
取样：取装有0~10 cm表层土壤环刀中的土壤样品约50 g，迅速装入铝盒内，利用电子天平，称量铝盒与新鲜土壤样品的质量，记为W2。
烘干：打开铝盒盖子（盖子放在铝盒旁边），放在105℃的恒温烘箱内烘8小时，盖好盖子，将铝盒置于干燥器内冷却30分钟，称重。
恒重：打开铝盒盖子，放在105℃的恒温烘箱内再次烘干3-5小时，盖好盖子，将铝盒置于干燥器内冷却30分钟，称重。若前后两次称重相差不超过3 mg，即可认为已达到恒重，记为W3；土壤含水量计算：
W%=（W2-W3）/（W2-W0）×100% (1)
式中，W指土壤含水量（%），W0指铝盒质量（g），W2指铝盒及新鲜土壤样品的质量（g），W3指铝盒及烘干土壤样品的质量（g）。
共获得22个表层土壤（0-10 cm）有效样点数据，其中，11个样点数据通过回归分析法，用于构建土壤含水量遥感估算模型，剩余的11个样点数据作为模型精度的验证数据，样点具体分布见附图1。
（b）遥感数据
中分辨率成像光谱仪MODIS(MODerate-resolution ImagingSpectroradiometer)是搭载在terra和aqua卫星上的一个重要的传感器，两颗卫星相互配合，每1-2天可重复观测整个地球表面，得到36个波段的观测数据，全球可以免费接收MODIS卫星遥感数据。MYD09A1，是MODIS 8天地表反射率合成产品，包含波段1至波段7的地表反射率，分辨率为500 m。为了确保遥感数据与实测数据在时间上的一致性，从NASA 网站下载了2016年10月7日的MYD09A1数据。
（2）辅助参数计算
（a）地表反照率计算
α = 0.16α1 + 0.29α2 + 0.243α3 + 0.116α4 + 0.112α5 + 0.08α7 - 0.0015 (2)
式中，a i (i=1,2,3,4,5 and 7) 是MODIS 产品MYD09A1的波段i。
（b）植被覆盖度计算
研究表明，植被覆盖度C和归一化植被指数NDVI之间存在极显著的线性相关关系，可以通过建立二者之间的转换关系，直接提取植被覆盖度信息，使用广泛应用的像元二分模型来计算。
C=（NDVI - NDVI s ）/（NDVI v - NDVI s ） (3)
式中，NDVI v 为植被覆盖部分的NDVI值，NDVI s 为土壤部分的 NDVI 值，C为植被覆盖度。
NDVI 按照下式计算：
NDVI=（NR - RED）/（NR + RED） (4)
式中，RED，NR分别为MODIS MYD09A1产品中第一波段（红光）和第二波段（近红外）的反射率。
在ENVI5.1软件平台下，计算NDVI值的累积频率，以 2% 的置信度截取 NDVI 的上下限阈值分别代表 NDVI v 和NDVI s 。
（3）从地表反照率中分解出裸土反射率
（a）反照率与土壤含水量关系
地表反照率与土壤含水量的关系可以表示为：
α = me np (5)
其中，a 是地表反照率，P 是土壤含水量质量百分比，m ，n是待定系数。
对方程（5）取自然对数，整理得到：
P = (lnα)/n –(lnm)/n (6)
如上文论述，为了消除植被对土壤光谱的影响，用裸土反照率a s 代替方程（5）中的地表反照率a，得到：
α s = m 1 e n1p (7)
类似地，对方程（7），取自然对数，得到：
P = (lnα s )/n 1 –(lnm 1 )/n 1 (8)
方程（6），（8）分别定义了利用地表反照率与裸土反照率估算土壤含水量的方法。地表反照率 a，由方程（2）计算得到；裸土反照率 a，可以由梯形法得到。
（b）梯形方法
研究表明，地表温度-植被指数、反照率-植被指数、反照率-植被盖度等的散点图在二维空间上形成梯形或三角形，因此得名为梯形或三角形方法。地表反照率是植被盖度与土壤含水量的函数，并且反照率-植被覆盖度散点图构成的包络线呈现梯形，如附图2所示。
四个顶点a sd 、a vd 、a vw ，a sw 构成了梯形特征空间。点a sd ，a sw 分别代表植被覆盖率为0时的裸土地表的最高与最低反照率，而点a vd ，a vw 分别表示全植被覆盖时植被表面的最高与最低反照率。由a sd ，a vd 两点决定的上边界线为“理论干边”，表示土壤水分极干的情况，是给定植被盖度条件下完全干旱土地对应的最高反照率的极限。现实条件下，由于研究区范围与地理位置、地表性状等的影响，“理论干边”往往无法准确地获得，通常用“实际干边”来代替，即对应覆盖度变化得到的地表反照率实际最高值组成的上边界线。由a sw ，a vw 两点决定的上边界线为最低反照率线，表示地表水分充足的状况，同样，在应用中通常用“实际湿边”来代替。位于“理论干边”，“理论湿边”，“实际干边”，“实际湿边”中间的过渡线（图中用虚线表示）近似为直线，且同一条线上的点具有相同的土壤湿度状况，即土壤等湿度线。当确定了梯形的4个顶点后，土壤等湿度线可以由线性内插方法获得。
（c）地表反照率分解
研究表明，遥感像元包含植被与土壤2类地物的反照率，混合像元地表反照率是土壤和植被反照率的加权和。
α =(1-C)α s + Cα v （9）
式中，a为混合像元地表反照率，C为植被覆盖度，a s 、a v 分别表示土壤、植被部分的反照率。以C对上式进行求导，获得新方程：
dα/dC =α v - α s （10）
由式（9）、（10），解得a s 、a v ：
α s = α– C × dα/dC （11）
α v = α + (1-C) dα/dC （12）
a，C可以由遥感数据获得，为已知项。欲求得a s ，需解得da/dC。da/dC是梯形框架中土壤等湿度线的斜率，令
k = dα/dC （13）
通过线性内插的假设，得到方程(14)：
k i =（α - α i,min ）（k d -k w ）/（α i,max - α i,min ）+ k w （14）
从（11）、（14）可以得到：
α=α s + C × [（α - α i,min ）（k d -k w ）/（α i,max - α i,min ）+ k w ] （15）
式中，k i 表示第i条土壤等湿度线的斜率，k d ，k w 分别表示“理论干边”和“理论湿边”的斜率值，实际应用时，通常用“实际干边”和“实际湿边”的斜率来代替，其余参数及其含义见附图2。
（d）地表反照率-植被覆盖度梯形空间绘制
为了得到方程（15）中k d ，k w ，必须借助于地表反照率-植被覆盖度散点图。而欲求地表反照率-植被覆盖度的散点图，需纳入植被覆盖度和反照率的极值（极大值，极小值）与其它值。因地表反照率的大小与地物类型密切相关，因此，结合植被覆盖度和具有特殊反照率的典型地物（如沙漠，水体），在MODIS假彩色合成图像（波段1、2、3）上选取采样区，以提取a和C。
沙漠可以视为十分干旱的地表，其反照率近似等于极干裸地（植被覆盖度近似等于0）的反照率。实验区青海湖湖东洱海周边沙漠广泛分布，选取了2个典型样区。大范围水体可以被看作是非常湿的地表，其地表反照率近似等于极湿裸地（植被覆盖度近似等于0）的反照率，其替代可行性已经得到了科学论证，在青海湖湖体中选取了2个典型样区。为了充分利用不同植被覆盖度条件下地物所蕴含的反照率信息，将植被覆盖度划分为无植被、较少、较多覆盖、基本全覆盖4个等级，然后在各个等级中选取若干个样区（如附图3）。
上述样区内，按照使用的MODIS数据空间分辨率，生成格网点（500 m×500 m），分别提取格网点地理位置上相应的地表反照率a、植被覆盖度C值。删除植被盖度异常值（负数以及大于1的值，由MODIS NDVI数据异常所致），然后绘制散点图，得到试验区实际地表反照率-植被覆盖度梯形特征空间（附图4）。
（e）裸土反照率计算
根据上述实际地表反照率-植被覆盖度梯形空间，得到其梯形的4个顶点坐标：a max （0.0321, 0.3220），B（0.4833, 0.1575），C（0.4833, 0.1155），D（0.0321, 0.0735）。据此，就很容易得到实际干边斜率K AB =-0.3646，即式（15）中的k d ，“实际湿边”斜率K DC =0.0931，即式（15）中的kw。依据方程（13）、（14）、（15），得到土壤等湿度线k i 与去除植被影响的裸土反照率a s。
k i = -1.8419α + 0.2285 (16)
α s = α (1 + 1.8419C) - 0.2285C (17)
式（17）表明，裸土反照率只与植被覆盖度和混合像元地表反照率有关，没有其它中间过渡参数。借助于斜率，实现了从混合像元地表反照率分解出裸土反照率，其分解方法巧妙，物理意义明确，这是本技术的一个重要创新。依据卫星遥感获取的地表反照率是植被盖度和土壤含水量的函数这一基本原理，提出了梯形特征空间，并构建了实际地表反照率-植被覆盖度梯形特征空间。由该梯形特征空间，获得了本技术发明的关键变量土壤等湿度线的斜率。通过该斜率，经过数学推导，巧妙地得到了裸土反照率。因考虑了植被覆盖度，进而去除了植被影响，恢复了土壤本身的反照率，有望提高土壤含水量估算与制图精度。
（4）土壤含水量遥感估算模型与制图
（a）估算模型构建
从方程（6）、（8）很容易发现，土壤含水量SM与反照率（地表反照率、裸土反照率）的自然对数呈线性关系。因此，应用回归分析法构建土壤含水量遥感估算模型(如附图5)，同时，用线性回归的决定系数R 2评价模型稳定性。
地表反照率的土壤含水量线性估算模型：
P = -7.6668lnα + 6.172, R2 = 0.3254 （18）
裸土反照率土壤含水量线性估算模型：
P = -9.9943lnα s + 2.8518, R2 = 0.807 （19）
式中，P为遥感估算的土壤含水量质量百分比，a为MODIS地表反照率，a s为裸土反照率。
通常，线性回归模型的优劣取决于两个方面，一是模型的决定系数R 2，R 2越大，模型越好。二是拟合模型点群的趋势，若点群中，更多的点通过或更加接近拟合直线，则模型优。因去除植被覆盖影响后的裸土反照率值小于混合像元反照率（即地表反照率）值，建模点在图5.a相对离散，而在图5.b中相对集中（图5.a，5. b的纵横坐标单位与刻度相同），因此，在图5.b中，更多的点通过或更加接近拟合直线。另外，从模型的判定系数R 2来看，裸土反照率估算模型的最大。因此，去除植被覆盖影响的裸土反照率估算模型，取得了较好的效果。
（b)土壤含水量制图
据（19）式，在ArcGIS 10软件平台下，借助栅格计算工具得到像元尺度的表层（0~10 cm）土壤含水量，并绘制成空间分布图（如图附6）。青海湖流域的东南部与西北部土壤含水量较低。东南部沙漠广布，其含水量较低。西北部多山，海拔较高，土壤类型为高寒草甸、寒漠土，其含水量较低。
（5）估算方法准确性评价
为了比较传统地表反照率与本方法提出的裸土反照率估算土壤含水量的优缺点，采用附图7所示均方根误差（Root Mean Squared Error ,RMSE），平均绝对百分比误差（Mean Absolute Percentage Error，MAPE）和希尔不等系数（Theil InequalityCoefficient，TIC）评价。附图7中，SM est,i， SM meas,i 分别为土壤含水量模型估算值，实测值，n为实测验证点的个数，本技术验证中，选取了11个验证点，即n为11。
这三个指标值越小，说明模型精度高，估算效果越好。从图8（表1）可以看出，裸土反照率模型的精度（MSE=4.20，MAPE=22.75%，TIC=0.67）明显高于传统的地表反照率模型的精度（RMSE=4.66，MAPE=25.46%，TIC=0.74），进一步证明了本方法的优越性。
有益效果：
由于植被（植被覆盖）影响了遥感传感器对地表土壤的可见性，从地表反照率（属于混合像元）中难以区分出裸土（纯土壤）的反照率。因此，传统的基于地表反照率的土壤含水量遥感估算方法，将植被与裸土的反照率视为土壤反照率，因而其估算的可靠性和稳定性难以保证。如裴浩等于1999年，利用气象卫星NOAA/AVHRR遥感数据，构建了土壤含水量线性估算模型：W=2.49+ 0.943P，W是土壤含水量，P是表观热惯量，P=(1-A)/(T max -T min )。其中，A是全反照率（由NOAA/AVHRR资料中午时次的第一、二通道的探测值计算得到，其含义与本技术中的地表反照率相同），最高最低温度T max 和T min 分别由中午、午夜或清晨两个时次的AVHRR卫星传感器的第四通道探测值推算得到。该方法未考虑植被的影响，将植被和裸土的反照率视为土壤的反照率，估算结果不可靠；地表温度T max 和T min 的推算受地表特征（主要由地物光谱比辐射率决定，地物不同，其光谱比辐射率不同）和大气状态（主要由大气透过率决定）的共同影响。在较大区域，地表特征在空间上变化很大，大气状态更是瞬息万变，大气透过率的不确定性较大，估算结果不稳定。因此，利用传统的地表反照率方法估算土壤含水量，其可靠性、稳定性难以保证，因而估算方法难以广泛地推广应用。
本方法通过植被覆盖-地表反照率构成的2维梯形空间，成功地分离出了裸土反照率，并依据反照率与土壤含水量的指数关系α s = m 1 e n1p ，构建了土壤含水量估算模型：P =(lnα s )/n 1 –(lnm 1 )/n 1 。该方法直接运用裸土反照率a s来计算土壤含水量P，物理意义明确，凸显了卫星遥感估算土壤含水量的实质内涵；克服了传统方法（如裴浩等.1999）需要计算其他中间参数，没有多余的中间环节，简单易行，避免了误差在参数间的传递。因此也显著提高了土壤含水量的估算与制图的精度，能够满足区域干旱监测、生态农业、精准农业，陆地生态系统、气候变化等领域土壤含水量的科学研究及实际应用需求。
附图说明
图1是本发明青海湖流域及地面实测土壤采样点位置图。
图2是本发明地表反照率与植被覆盖度形成的梯形特征空间示意图。
图3是本发明考虑了沙漠水体以及不同植被盖度，利用目视法选取典型采样区分布图。
图4是本发明试验区实际地表反照率-植被覆盖度梯形特征空间示意图。
图5是本发明土壤有含水量遥感估算模型图。
图6是本发明实验区表层（0~10 cm）土壤含水量遥感估算图。
图7均方根误差、平均绝对百分比误差和希尔不等系数的计算公式。
图8是本发明两种模型土壤含水量估算精度比较。
具体实施方式
为使发明的目的、技术方案和优点更加清楚，下面结合附图对本发明的具体实施方式进行详细说明。这些优选实施方式的示例在附图中进行了例示。附图中所示和根据附图描述的本发明的实施方式仅仅是示例性的，并且本发明并不限于这些实施方式。
在此，还需要说明的是，为了避免因不必要的细节而模糊了本发明的技术方案，在附图中仅仅示出了与根据本发明的方案密切相关的结构和/或处理步骤，而省略了关系不大的其他细节。
实施例1
本实施例提供了一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法，本方法的步骤包括（1）材料及数据；（2）辅助参数计算（地表反照率、植被覆盖度计算）；（3）从地表反照率中分解出裸土反照率；（4）土壤含水量遥感估算模型与制图；（5）估算方法准确性评价。
（1）材料及数据
（a）土壤含水量实测数据
2016年9月27日至10月4日进行野外土壤样点的采集。采用环刀法取地表 (0~10cm、10~20 cm、20~30 cm)的土样。样品带回实验室后，按照下列步骤测量：
准备工作：在室内将铝盒编号并称重，重量记为W0。
取样：取装有0~10 cm表层土壤环刀中的土壤样品约50 g，迅速装入铝盒内，利用电子天平，称量铝盒与新鲜土壤样品的质量，记为W2。
烘干：打开铝盒盖子（盖子放在铝盒旁边），放在105℃的恒温烘箱内烘8小时，盖好盖子，将铝盒置于干燥器内冷却30分钟，称重。
恒重：打开铝盒盖子，放在105℃的恒温烘箱内再次烘干3-5小时，盖好盖子，将铝盒置于干燥器内冷却30分钟，称重。若前后两次称重相差不超过3 mg，即可认为已达到恒重，记为W3。
土壤含水量计算：
W%=（W2-W3）/（W2-W0）×100% (1)
式中，W指土壤含水量（%），W0指铝盒质量（g），W2指铝盒及新鲜土壤样品的质量（g），W3指铝盒及烘干土壤样品的质量（g）。
共获得22个表层土壤（0-10 cm）有效样点数据，其中，11个样点数据通过回归分析法，用于构建土壤含水量遥感估算模型，剩余的11个样点数据作为模型精度的验证数据，样点具体分布见附图1。
（b）遥感数据
中分辨率成像光谱仪MODIS(MODerate-resolution ImagingSpectroradiometer)是搭载在terra和aqua卫星上的一个重要的传感器，两颗卫星相互配合，每1-2天可重复观测整个地球表面，得到36个波段的观测数据，全球可以免费接收MODIS卫星遥感数据。MYD09A1，是MODIS 8天地表反射率合成产品，包含波段1至波段7的地表反射率，分辨率为500 m。为了确保遥感数据与实测数据在时间上的一致性，从NASA 网站下载了2016年10月7日的MYD09A1数据。
（2）辅助参数计算
（a）地表反照率计算
α = 0.16α1 + 0.29α2 + 0.243α3 + 0.116α4 + 0.112α5 + 0.08α7 - 0.0015 (2)
式中，a i (i=1,2,3,4,5 and 7) 是MODIS 产品MYD09A1的波段i。
（b）植被覆盖度计算
研究表明，植被覆盖度C和归一化植被指数NDVI之间存在极显著的线性相关关系，可以通过建立二者之间的转换关系，直接提取植被覆盖度信息，使用广泛应用的像元二分模型来计算。
C=（NDVI - NDVI s ）/（NDVI v - NDVI s ） (3)
式中，NDVI v 为植被覆盖部分的NDVI值，NDVI s 为土壤部分的 NDVI 值，C为植被覆盖度。
NDVI 按照下式计算：
NDVI=（NR - RED）/（NR + RED） (4)
式中，RED，NR分别为MODIS MYD09A1产品中第一波段（红光）和第二波段（近红外）的反射率。
在ENVI5.1软件平台下，计算NDVI值的累积频率，以 2% 的置信度截取 NDVI 的上下限阈值分别代表 NDVI v 和NDVI s 。
（3）从地表反照率中分解出裸土反射率
（a）反照率与土壤含水量关系
地表反照率与土壤含水量的关系可以表示为：
α = me np (5)
其中，a 是地表反照率，P 是土壤含水量质量百分比，m ，n是待定系数。
对方程（5）取自然对数，整理得到：
P = (lnα)/n –(lnm)/n (6)
如上文论述，为了消除植被对土壤光谱的影响，用裸土反照率a s 代替方程（5）中的地表反照率a，得到：
α s = m 1 e n1p (7)
类似地，对方程（7），取自然对数，得到：
P = (lnα s )/n 1 –(lnm 1 )/n 1 (8)
方程（6），（8）分别定义了利用地表反照率与裸土反照率估算土壤含水量的方法。地表反照率 a，由方程（2）计算得到；裸土反照率 a，可以由梯形法得到。
（b）梯形方法
研究表明，地表温度-植被指数、反照率-植被指数、反照率-植被盖度等的散点图在二维空间上形成梯形或三角形，因此得名为梯形或三角形方法。地表反照率是植被盖度与土壤含水量的函数，并且反照率-植被覆盖度散点图构成的包络线呈现梯形，如附图2所示。
四个顶点a sd 、a vd 、a vw ，a sw 构成了梯形特征空间。点a sd ，a sw 分别代表植被覆盖率为0时的裸土地表的最高与最低反照率，而点a vd ，a vw 分别表示全植被覆盖时植被表面的最高与最低反照率。由a sd ，a vd 两点决定的上边界线为“理论干边”，表示土壤水分极干的情况，是给定植被盖度条件下完全干旱土地对应的最高反照率的极限。现实条件下，由于研究区范围与地理位置、地表性状等的影响，“理论干边”往往无法准确地获得，通常用“实际干边”来代替，即对应覆盖度变化得到的地表反照率实际最高值组成的上边界线。由a sw ，a vw 两点决定的上边界线为最低反照率线，表示地表水分充足的状况，同样，在应用中通常用“实际湿边”来代替。位于“理论干边”，“理论湿边”，“实际干边”，“实际湿边”中间的过渡线（图中用虚线表示）近似为直线，且同一条线上的点具有相同的土壤湿度状况，即土壤等湿度线。当确定了梯形的4个顶点后，土壤等湿度线可以由线性内插方法获得。
（c）地表反照率分解
研究表明，遥感像元包含植被与土壤2类地物的反照率，混合像元地表反照率是土壤和植被反照率的加权和。
α =(1-C)α s + Cα v （9）
式中，a为混合像元地表反照率，C为植被覆盖度，a s 、a v 分别表示土壤、植被部分的反照率。以C对上式进行求导，获得新方程：
dα/dC =α v - α s （10）
由式（9）、（10），解得a s 、a v ：
α s = α– C × dα/dC （11）
α v = α + (1-C) dα/dC （12）
a，C可以由遥感数据获得，为已知项。欲求得a s ，需解得da/dC。da/dC是梯形框架中土壤等湿度线的斜率，令
k = dα/dC （13）
通过线性内插的假设，得到方程(14)：
k i =（α - α i,min ）（k d -k w ）/（α i,max - α i,min ）+ k w （14）
从（11）、（14）可以得到：
α=α s + C × [（α - α i,min ）（k d -k w ）/（α i,max - α i,min ）+ k w ] （15）
式中，k i 表示第i条土壤等湿度线的斜率，k d ，k w 分别表示“理论干边”和“理论湿边”的斜率值，实际应用时，通常用“实际干边”和“实际湿边”的斜率来代替，其余参数及其含义见附图2。
（d）地表反照率-植被覆盖度梯形空间绘制
为了得到方程（15）中k d ，k w ，必须借助于地表反照率-植被覆盖度散点图。而欲求地表反照率-植被覆盖度的散点图，需纳入植被覆盖度和反照率的极值（极大值，极小值）与其它值。因地表反照率的大小与地物类型密切相关，因此，结合植被覆盖度和具有特殊反照率的典型地物（如沙漠，水体），在MODIS假彩色合成图像（波段1、2、3）上选取采样区，以提取a和C。
沙漠可以视为十分干旱的地表，其反照率近似等于极干裸地（植被覆盖度近似等于0）的反照率。实验区青海湖湖东洱海周边沙漠广泛分布，选取了2个典型样区。大范围水体可以被看作是非常湿的地表，其地表反照率近似等于极湿裸地（植被覆盖度近似等于0）的反照率，其替代可行性已经得到了科学论证，在青海湖湖体中选取了2个典型样区。为了充分利用不同植被覆盖度条件下地物所蕴含的反照率信息，将植被覆盖度划分为无植被、较少、较多覆盖、基本全覆盖4个等级，然后在各个等级中选取若干个样区（如附图3）。
上述样区内，按照使用的MODIS数据空间分辨率，生成格网点（500 m×500 m），分别提取格网点地理位置上相应的地表反照率a、植被覆盖度C值。删除植被盖度异常值（负数以及大于1的值，由MODIS NDVI数据异常所致），然后绘制散点图，得到试验区实际地表反照率-植被覆盖度梯形特征空间（附图4）。
（e）裸土反照率计算
根据上述实际地表反照率-植被覆盖度梯形空间，得到其梯形的4个顶点坐标：a max （0.0321, 0.3220），B（0.4833, 0.1575），C（0.4833, 0.1155），D（0.0321, 0.0735）。据此，就很容易得到实际干边斜率K AB =-0.3646，即式（15）中的k d ，“实际湿边”斜率K DC =0.0931，即式（15）中的kw。依据方程（13）、（14）、（15），得到土壤等湿度线k i 与去除植被影响的裸土反照率a s。
k i = -1.8419α + 0.2285 (16)
α s = α (1 + 1.8419C) - 0.2285C (17)
式（17）表明，裸土反照率只与植被覆盖度和混合像元地表反照率有关，没有其它中间过渡参数。借助于斜率，实现了从混合像元地表反照率分解出裸土反照率，其分解方法巧妙，物理意义明确，这是本技术的一个重要创新。依据卫星遥感获取的地表反照率是植被盖度和土壤含水量的函数这一基本原理，提出了梯形特征空间，并构建了实际地表反照率-植被覆盖度梯形特征空间。由该梯形特征空间，获得了本技术发明的关键变量土壤等湿度线的斜率。通过该斜率，经过数学推导，巧妙地得到了裸土反照率。因考虑了植被覆盖度，进而去除了植被影响，恢复了土壤本身的反照率，有望提高土壤含水量估算与制图精度。
（4）土壤含水量遥感估算模型与制图
（a）估算模型构建
从方程（6）、（8）很容易发现，土壤含水量SM与反照率（地表反照率、裸土反照率）的自然对数呈线性关系。因此，应用回归分析法构建土壤含水量遥感估算模型(如附图5)，同时，用线性回归的决定系数R 2评价模型稳定性。
地表反照率的土壤含水量线性估算模型：
P = -7.6668lnα + 6.172, R2 = 0.3254 （18）
裸土反照率土壤含水量线性估算模型：
P = -9.9943lnα s + 2.8518, R2 = 0.807 （19）
式中，P为遥感估算的土壤含水量质量百分比，a为MODIS地表反照率，a s为裸土反照率。
通常，线性回归模型的优劣取决于两个方面，一是模型的决定系数R 2，R 2越大，模型越好。二是拟合模型点群的趋势，若点群中，更多的点通过或更加接近拟合直线，则模型优。因去除植被覆盖影响后的裸土反照率值小于混合像元反照率（即地表反照率）值，建模点在图5.a相对离散，而在图5.b中相对集中（图5.a，5. b的纵横坐标单位与刻度相同），因此，在图5.b中，更多的点通过或更加接近拟合直线。另外，从模型的判定系数R 2来看，裸土反照率估算模型的最大。因此，去除植被覆盖影响的裸土反照率估算模型，取得了较好的效果。
（b)土壤含水量制图
据（19）式，在ArcGIS 10软件平台下，借助栅格计算工具得到像元尺度的表层（0~10 cm）土壤含水量，并绘制成空间分布图（如图附6）。青海湖流域的东南部与西北部土壤含水量较低。东南部沙漠广布，其含水量较低。西北部多山，海拔较高，土壤类型为高寒草甸、寒漠土，其含水量较低。
（5）估算方法准确性评价
为了比较传统地表反照率与本方法提出的裸土反照率估算土壤含水量的优缺点，采用附图7所示均方根误差（Root Mean Squared Error ,RMSE），平均绝对百分比误差（Mean Absolute Percentage Error，MAPE）和希尔不等系数（Theil InequalityCoefficient，TIC）评价。附图7中，SM est,i， SM meas,i 分别为土壤含水量模型估算值，实测值，n为实测验证点的个数，本技术验证中，选取了11个验证点，即n为11。
这三个指标值越小，说明模型精度高，估算效果越好。从图8（表1）可以看出，裸土反照率模型的精度（MSE=4.20，MAPE=22.75%，TIC=0.67）明显高于传统的地表反照率模型的精度（RMSE=4.66，MAPE=25.46%，TIC=0.74），进一步证明了本方法的优越性。
实施例2
实施例1是在遥感图像上，目视选取典型地物（需要遥感专业技术背景才能准确识别这些典型地物），进而得到这些典型地物地理位置上的植被盖度与地表反照率值，绘制了实际地表反照率-植被覆盖度梯形空间。为了使本方法更具普遍性，可以将研究区划定为一定数量的格网，然后提取每个格网位置上的植被盖度与地表反照率值，然后绘制实际地表反照率-植被覆盖度梯形空间。具体方法如下：
首先，通过ArcGIS 10软件的fishnet工具将研究区划分500 m*500 m的格网（如果研究区较大，可以划分成1 km*1 km或者10 km*10 km等格网）。
然后，通过ArcGIS 10软件的Extract MultiValues To Points工具，提取上述格网位置上的植被盖度与地表反照率值，再依据实施例1技术方案，进行土壤含水量的遥感估算。
有益效果：
卫星遥感技术是目前区域尺度上土壤含水量估算的主要技术手段。但由于卫星遥感估算实施的基本单元——像元，当有植被覆盖时，属于混合像元，包含了植被与土壤2类地物的光谱信息组份，像元的反照率是裸露土壤与植被反照率的加权和。传统的土壤含水量卫星遥感估算方法，将植被和裸土的混合反照率视为土壤的反照率，没有去除植被组份的光谱信息，遥感物理意义不明确，因而估算的精度与可靠性难以保证。
（1）基于混合像元的组份，本技术考虑了植被覆盖度的影响，将混合地表反照率抽象概括为植被覆盖度、植被反照率和裸土反照率的数学模型。以植被覆盖度对该数学模型进行导数运算，并从地表反照率-植被覆盖度梯形特征空间，引入了本技术的关键变量土壤等湿度线及其斜率，经过严格的数学推导，实现了混合地表反照率的分解，获得了裸土反照率。
（2）裸土反照率只与植被覆盖度和混合像元地表反照率有关。没有其它中间过渡参数，避免了传统方法（基于地表反照率法如裴浩等.1999和像元分解法）需要大量的参数，这些参数或是通过地面测量获取（通常难以全部获取），或是通过模型模拟得到（导致误差在多个参数间的传递）。而直接运用本技术获取的裸土反照率估算土壤含水量，遥感物理意义明确，且保证了估算结果的可靠性与稳定性。
（3）突破了某些传统方法的应用限制，提高了推广应用价值。如：“光学植被盖度”方法应用的前提是区域内有全裸地（无植被）和植被全覆盖地区（植被覆盖度C达到100%），本技术可以运用到任何较大的区域。
（4）简单易操作。本技术只需在遥感图像上选定典型地物或者根据实施例2的操作，再依据本专利的相应操作步骤，即使没有遥感专业背景的人员也可以实现土壤含水量的卫星遥感估算。
（5）应用领域广泛。本技术能够满足区域干旱监测、生态农业、精准农业，陆地生态系统、气候变化等领域的土壤含水量科学研究及实际应用需求。
以上所述仅是本申请的具体实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本申请原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本申请的保护范围。
Claims (1)
1.一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法，其特征在于，包括如下步骤：
（1）材料及数据收集，包括：土壤含水量实测数据和遥感数据收集；
（2）辅助参数计算，包括：
地表反照率计算：α = 0.16α1 + 0.29α2 + 0.243α3 + 0.116α4 + 0.112α5 + 0.08α7 -0.0015，式中，a i 是MODIS 产品MYD09A1的波段i，i=1，2，3，4，5，7；
植被覆盖度计算：植被覆盖度C和归一化植被指数NDVI之间存在极显著的线性相关关系，C=（NDVI - NDVI s ）/（NDVI v - NDVI s ），式中，NDVI v 为植被覆盖部分的NDVI值，NDVI s 为土壤部分的NDVI值，C为植被覆盖度；
（3）从地表反照率中分解出裸土反照率：
按照使用的MODIS数据空间分辨率，生成格网点，分别提取格网点地理位置上相应的地表反照率a、植被覆盖度C值，删除植被盖度异常值，然后绘制散点图，得到试验区实际地表反照率-植被覆盖度梯形特征空间，根据上述实际地表反照率-植被覆盖度梯形特征空间，得到其梯形的4个顶点坐标，进而得到实际干边斜率、实际湿边斜率，依据公式α=α s +C×[（α-α i,min ）（k d -k w ）/（α i,max -α i,min ）+ k w ]，得到去除植被影响的裸土反照率a s：α s =α(1 +1.8419C) - 0.2285C；
其中，kd为实际干边斜率，kw为实际湿边斜率，C为植被覆盖度，a为地表反照率，as为去除植被影响的裸土反照率，α i,max 和α i,min 分别表示为梯形特征空间中植被覆盖度最小时的地表反照率的最大值和最小值；
（4）土壤含水量遥感估算模型与制图：
土壤含水量与反照率的自然对数呈线性关系，应用回归分析法构建裸土反照率的土壤含水量线性估算模型：P = -9.9943lnα s + 2.8518, R2 = 0.807；其中，P为土壤含水量，as为去除植被影响的裸土反照率，R 2为线性回归的决定系数；
（5）估算模型准确性评价：采用均方根误差，平均绝对百分比误差和希尔不等系数对估算模型进行评价。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810505327.2A CN108717044B (zh) | 2018-05-24 | 2018-05-24 | 一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810505327.2A CN108717044B (zh) | 2018-05-24 | 2018-05-24 | 一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108717044A CN108717044A (zh) | 2018-10-30 |
| CN108717044B true CN108717044B (zh) | 2021-07-30 |
Family
ID=63900413
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810505327.2A Expired - Fee Related CN108717044B (zh) | 2018-05-24 | 2018-05-24 | 一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108717044B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109781635B (zh) * | 2018-12-29 | 2020-07-31 | 长沙天仪空间科技研究院有限公司 | 一种分布式遥感卫星系统 |
| CN110245327B (zh) * | 2019-03-28 | 2023-04-18 | 国智恒北斗好年景农业科技有限公司 | 一种基于gf-1数据重构的小麦单产遥感估算方法 |
| CN110610054B (zh) * | 2019-09-23 | 2021-03-23 | 北京师范大学 | 一种土壤湿度长方体反演模型构建方法及系统 |
| CN110688609B (zh) * | 2019-09-26 | 2021-08-31 | 中国水利水电科学研究院 | 一种黄土塬区浅层地下水补给-排泄单元划分方法 |
| CN110618144B (zh) * | 2019-09-26 | 2021-01-05 | 中国水利水电科学研究院 | 一种快速测定黄土塬区泉眼位置的方法 |
| CN111815102B (zh) * | 2020-04-15 | 2024-01-26 | 中国环境科学研究院 | 一种基于空间技术的生物多样性综合调查抽样方法 |
| CN114692378A (zh) * | 2020-12-31 | 2022-07-01 | 南开大学 | 一种Tm-fc梯形框架内极端温度组合的计算方法 |
| CN113553549B (zh) * | 2021-07-26 | 2023-04-14 | 中国科学院西北生态环境资源研究院 | 一种植被覆盖度反演方法、装置、电子设备及存储介质 |
| CN114397251A (zh) * | 2021-12-31 | 2022-04-26 | 核工业北京地质研究院 | 一种植被覆盖区铀矿化蚀变信息遥感识别方法 |
| CN115100609B (zh) * | 2022-08-26 | 2022-11-22 | 北京江河惠远科技有限公司 | 特高压施工扰动范围提取方法和系统 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101187630A (zh) * | 2007-12-05 | 2008-05-28 | 北京大学 | 一种农田干旱监测方法 |
| CN103675234A (zh) * | 2013-12-11 | 2014-03-26 | 中国科学院遥感与数字地球研究所 | 基于地表温度与植被指数特征空间的干旱指数快速监测方法 |
| US10564316B2 (en) * | 2014-09-12 | 2020-02-18 | The Climate Corporation | Forecasting national crop yield during the growing season |
| CN106779067B (zh) * | 2016-12-02 | 2019-04-05 | 清华大学 | 基于多源遥感数据的土壤湿度重建方法和系统 |
| CN107036968B (zh) * | 2016-12-27 | 2019-09-27 | 西安科技大学 | 一种土壤湿度实时监测方法 |
| CN107389895B (zh) * | 2017-06-08 | 2019-08-30 | 环境保护部卫星环境应用中心 | 土壤水分混合型遥感反演方法及系统 |
-
2018
- 2018-05-24 CN CN201810505327.2A patent/CN108717044B/zh not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN108717044A (zh) | 2018-10-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108717044B (zh) | 一种去除植被覆盖影响的表层土壤含水量卫星遥感估算方法 | |
| Ge et al. | Modeling alpine grassland cover based on MODIS data and support vector machine regression in the headwater region of the Huanghe River, China | |
| Mallick et al. | Land surface emissivity retrieval based on moisture index from LANDSAT TM satellite data over heterogeneous surfaces of Delhi city | |
| D'Odorico et al. | Intercomparison of fraction of absorbed photosynthetically active radiation products derived from satellite data over Europe | |
| Long et al. | A two-source trapezoid model for evapotranspiration (TTME) from satellite imagery | |
| Fang et al. | An inter-comparison of soil moisture data products from satellite remote sensing and a land surface model | |
| KR101914061B1 (ko) | 인공위성에 의한 열섬특성 분석방법 | |
| Restrepo-Coupe et al. | MODIS vegetation products as proxies of photosynthetic potential along a gradient of meteorologically and biologically driven ecosystem productivity | |
| WO2018107245A1 (en) | Detection of environmental conditions | |
| CN110927120B (zh) | 一种植被覆盖度预警方法 | |
| Foolad et al. | Comparison of the automatically calibrated Google Evapotranspiration Application—EEFlux and the manually calibrated METRIC application | |
| Blumberg et al. | Quantifying the accuracy and uncertainty of diurnal thermodynamic profiles and convection indices derived from the Atmospheric Emitted Radiance Interferometer | |
| WU et al. | Remotely sensed estimation and mapping of soil moisture by eliminating the effect of vegetation cover | |
| Rempel et al. | Object-based metrics for forecast verification of convective development with geostationary satellite data | |
| Alexandrov et al. | Remote sensing of atmospheric aerosols and trace gases by means of multifilter rotating shadowband radiometer. Part II: Climatological applications | |
| Bechtel et al. | Classification and modelling of urban micro-climates using multisensoral and multitemporal remote sensing data | |
| Lou et al. | An effective method for canopy chlorophyll content estimation of marsh vegetation based on multiscale remote sensing data | |
| Kelly et al. | Modelling and upscaling ecosystem respiration using thermal cameras and UAVs: Application to a peatland during and after a hot drought | |
| Li et al. | Derivation of the Green Vegetation Fraction of the Whole China from 2000 to 2010 from MODIS Data | |
| Yu et al. | Hyperspectral database prediction of ecological characteristics for grass species of alpine grasslands | |
| Boyd et al. | Exploring spatial and temporal variation in middle infrared reflectance (at 3.75@ m) measured from the tropical forests of west Africa | |
| Paul et al. | Evaluating surface energy balance system (SEBS) using aircraft data collected during BEAREX07 | |
| Jakobsson | The suitability of using Landsat TM-5 Images for estimating chromophoric dissolved organic matter in subarctic Lakes | |
| Cai | Vegetation Observation in the Big Data Era: Sentinel-2 data for mapping the seasonality of land vegetation | |
| Choudhury et al. | Potential role of landsat satellite data for the evaluation of land surface temperature and assessment of urban environment |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210730 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |