一种综合植被和生境特征的草原区土壤水分快速反演方法与流程

1.本发明涉及表层土壤水分遥感应用技术领域，具体的说是一种基于植被冠层参量和生境特征的干旱半干旱草原区表层土壤水分含量快速反演的方法。


背景技术：

2.土壤水分通常指存储在非饱和土壤中的水量，又称土壤含水量或土壤湿度，一般采用土壤重量含水量和土壤体积含水量表示。土壤水分是地表和大气能量交换以及水分交换过程中的重要环节，可以调节空气的温度和湿度。土壤水分是影响地表降雨与径流的重要参量，也是理解地气间以水和能量为媒介的交互反馈作用的关键参量，对该参量的监测对于生态系统健康状况、农业发展等具有重要意义。
3.土壤水分监测传统的方式一般是通过便于携带且易于安装的仪器进行土壤水分数据的采集，但是此方法在获取土壤水分数据的过程中需要投入大量的人力和物力，同时测量的数据代表范围小，难以满足大范围土壤水分的持续监测。遥感技术具有信息量大、快速易行、无污染、无破坏等优点，尤其是高光谱遥感技术能快速获取土壤的连续光谱信息，探测土壤水分变化的细微差异，为土壤水分的定量反演提供了可能。在分析过程中，通过分析不同土壤类型区供试样方中的光谱形状特征，利用相关系数比较，发现1400nm及1900nm附近的吸收峰是随着土壤水分变化显著的波段。由此说明，通过挖掘高光谱曲线的特征可以找出其与土壤水分之间的规律并实现土壤水分的反演。


技术实现要素：

4.本发明针对传统土壤水分测量方法测定周期长、成本高、不能实现大范围快速监测等问题，提出一种基于植被冠层参量和生境特征的干旱半干旱草原区表层土壤水分含量快速反演的方法，该方法可进一步应用到卫星影像，为实现区域快速重复监测土壤水分含量，提供理论和技术支持。
5.这种基于植被冠层参量和生境特征的干旱半干旱草原区表层土壤水分含量快速反演的方法，包括以下步骤：
6.步骤s1、对研究区进行野外实地调查，采集土壤水分数据、植被冠层反射光谱数据、叶面积指数数据、叶片叶绿素含量数据、土壤温度数据；
7.步骤s2、对野外实测冠层参量数据以及原始反射率数据进行预处理，随后对反射率数据进行savitzky-golay滤波以消除光谱噪声；
8.步骤s3、利用滤波后的反射率数据进行相关指数(归一化植被指数ndvi、增强型植被指数evi、比值植被指数sr、红边归一化植被指数ndre、水分胁迫指数msi、归一化水分指数ndwi、简单比值水分指数srwi)的计算；
9.步骤s4、采用竞争性自适应重加权算法(competitive adaptive reweighted sampling，cars)对经过滤波后的反射率数据进行波段选择；
10.步骤s5、利用相关性分析的方法，分析植被生物理化参量(叶面积指数、冠层叶绿
素含量、ndvi、evi、sr、ndre)、冠层含水量指标(msi、ndwi、srwi)以及生境特征(土壤温度)与土壤水分之间的相关性，筛选出相关性好的几种参量；
11.步骤s6、将筛选出的数据(含特征波段和优选指数)进行训练集和验证集的划分；确定bp神经网络结构，构建bp神经网络模型，从而实现土壤水分的反演。
12.其中，本发明的技术方案还进一步包括：
13.执行步骤s1时，主要通过野外实地测量得到研究区的土壤水分数据、植被冠层反射光谱数据、叶面积指数数据、叶片叶绿素含量数据、土壤温度数据，具体步骤如下；
14.步骤s1.1、土壤水分数据采用tdr-350进行测量，该仪器可以快速简便地测量表层土壤水分，通过新功能的改进，能够对土壤电导率进行测量，修正土壤水分读数。测量样方土壤水分的过程中，同时采集样方四角以及中心点位置的0-10cm深度土壤水分数据及土壤温度数据，随后分别将采集的样方四角以及中心点位置的土壤水分数据以及土壤温度数据进行算数平均作为该样方土壤水分数据以及土壤温度数据并加以记录，其中土壤水分数据单位为％vwc，土壤温度数据单位为℃；
15.步骤s1.2、植被冠层反射光谱数据采用svc光谱仪进行采集，其光谱范围为350nm-2500nm。在实测过程中，选择在天气晴朗、无风或风速很小的情况下进行。在每次观测前先进行白板校正，测量过程中，探头对准样方内植被冠层，并保持1m的距离，每个样本点记录5条光谱数据并取平均值；
16.步骤s1.3、叶面积指数数据采用lai-2200进行采集，在各样方点观测过程中，一般每个采集点测量一次天空光和4次植被下数据，尽量选择天空光变化较小的时段进行。同时在观测lai时要变化角度，以确保观测值的精确性，最终对每个像元内的观测值取平均值；
17.步骤s1.4、叶片叶绿素含量数据使用spad-502进行采集，其测量面积为2mm
×
3mm(厚度不超过1.2mm)，在测量前，先对仪器进行校准，校准后开始对样方内叶片的叶绿素含量进行测量。在测量过程中，尽量避免选择过厚的样品，如果测定有较多叶脉的叶片，可进行多次测量并取其平均值，此时记录的叶片叶绿素含量为spad值。
18.执行步骤s2时，具体步骤如下；
19.步骤s2.1、将采集的植被冠层反射光谱数据“*.sig”文件转换为操作性强的“*.xls”文件，每个xls文件中包含波段信息及其对应的反射率值；
20.步骤s2.2、对同一样方内的多条光谱数据进行平均计算，用平均后的光谱数值代表该样方的光谱属性。
21.步骤s2.3、采用savitzky-golay滤波的方法对光谱数据进行多点平滑处理，在较好地保持原始反射率光谱信息的同时，有效地降低和消除数据中的噪声。
22.步骤s2.4、通过换算，将叶片叶绿素含量换算为冠层叶绿素含量。其计算公式为：
23.ccc＝lcc
×
lai
ꢀꢀ
(1)
24.其中，ccc表示冠层叶绿素含量，单位为μg/cm2；lcc表示叶片叶绿素含量，单位为μg/cm2；lai表示叶面积指数，单位为m2/m2；
25.执行步骤s3时，主要包括几种指数的计算。在指数计算过程前，针对每一种指数对波段进行组合以计算该指数对应的所有波段组合，选出与土壤水分相关性最好的波段组合作为后续指数计算的依据。
26.指数计算的具体过程为：
27.步骤s3.1、对归一化植被指数的计算，其计算公式为：
[0028][0029]
其中，ndvi表示归一化植被指数；ρ
741
表示在波长741nm处的反射率；ρ
676
表示在波长676nm处的反射率；
[0030]
步骤s3.2、对归一化水分指数的计算，其计算公式为：
[0031][0032][0033][0034]
其中，ndwi表示归一化水分指数；ρ
800
表示在波长800nm处的反射率；ρ
1137
表示在波长1137nm处的反射率；ρ
500
表示在波长500nm处的反射率；ρ
741
表示在波长741nm处的反射率；ρ
1000
表示在波长1000nm处的反射率；
[0035]
步骤s3.3、对简单比值水分指数的计算，其计算公式为：
[0036][0037][0038][0039]
其中，srwi表示简单比值水分指数；ρ
800
表示在波长800nm处的反射率；ρ
1137
表示在波长1137nm处的反射率；ρ
500
表示在波长500nm处的反射率；ρ
741
表示在波长741nm处的反射率；ρ
1000
表示在波长1000nm处的反射率；
[0040]
步骤s3.4、对增强型植被指数的计算，其计算公式为：
[0041][0042]
其中，evi表示增强型植被指数；ρ
741
表示在波长741nm出的反射率；ρ
676
表示在波长676nm处的反射率；ρ
408
表示在波长408nm处的反射率；
[0043]
步骤s3.5、对比值植被指数sr的计算，其计算公式为：
[0044][0045]
其中，sr表示比值植被指数；ρ
741
表示在波长741nm处的反射率；ρ
721
表示在波长721nm处的反射率；
[0046]
步骤s3.6、对红边归一化植被指数的计算，其计算公式为：
[0047][0048]
其中，ndre表示红边归一化植被指数；ρ
741
表示在波长741nm处的反射率；ρ
690
表示
在波长690nm处的反射率；
[0049]
步骤s3.7、对水分胁迫指数的计算，其计算公式为：
[0050][0051]
其中，msi表示水分胁迫指数；ρ
1464
表示在波长1464nm处的反射率；ρ
800
表示在波长800nm处的反射率；
[0052]
执行步骤s4时，主要通过cars算法进行波段的选择，具体过程如下；
[0053]
步骤s4.1、利用matlab软件，将测得的m
×
1的土壤水分数据以及m
×
n的光谱数据作为输入数据，进行cars计算；
[0054]
cars算法是一种结合蒙特卡洛与pls模型回归系数的特征变量选择方法。cars算法中，每次通过自适应加权采样保留pls模型中回归系数绝对值权重较大的点作为新的子集，去掉权值较小的点，然后基于新的子集建立pls模型，经过多次运算，选择pls模型交互验证均方根误差(rmsecv)最小的子集中的波长作为特征波长。
[0055]
步骤s4.2、将cars算法筛选出的波段进行记录并保存。执行步骤s5时，其具体步骤如下；
[0056]
步骤s5.1、将测得的土壤水分数据、叶面积指数数据、土壤温度数据，以及计算得到的冠层叶绿素含量数据、ndvi、ndwi、srwi、evi、sr、ndre、msi数据以样方为单位进行整合；
[0057]
步骤s5.2、对土壤水分和植被生物理化参量(叶面积指数、冠层叶绿素含量、ndvi、evi、sr、ndre、冠层含水量指标(msi、ndwi、srwi)以及土壤温度进行相关性分析。
[0058]
步骤s5.3、筛选出相关系数绝对值较高的指数数据进行记录并保存。
[0059][0060]
其中，σ
x，y
表示变量x，y的相关系数；μ
x
和μy分别表示变量x，y的均值(期望)，σ
x
、σy分别表示变量x，y的标准差。σ
x，y
取值范围为[-1，1]，绝对值越大，相关性越强。执行步骤s6时，主要进行数据集的构建并划分，随后构建bp神经网络结构并进行土壤水分的反演，具体步骤如下：
[0061]
步骤s6.1、构建含有上述筛选出的光谱波段数据以及各种指数数据的数据集，并按照7：3的比例分别划分测试集和验证集。
[0062]
步骤s6.2、确定bp神经网络结构，以筛选的波段数据和指数数据作为输入值，以土壤水分作为输出值，构建bp神经网络模型。
[0063]
步骤s6.3、利用构建的bp神经网络模型进行土壤水分的反演，将反演结果同实测值进行rmse的计算评价精度，rmse越小，反演效果越好。
[0064][0065]
本发明的优点及有益效果在于：本发明利用cars算法对savitzky-golay滤波后的光谱数据进行特征波段选择，以选出能够反映表层土壤水分的波段数据；同时利用滤波后的光谱数据计算ndvi、ndwi、srwi、evi、sr、ndre、msi指数，随后采用相关性分析的方法，计
算植被生物理化参量(叶面积指数、冠层叶绿素含量、ndvi、evi、sr、ndre)、冠层含水量指标(msi、ndwi、srwi)以及土壤温度与表层土壤水分的相关系数，选出相关系数绝对值较高的几种指数。最后通过特征波段与优选指数结合表层土壤水分数据构建bp神经网络模型对表层土壤水分进行反演。本发明快速高效的进行干旱半干旱草原区表层土壤水分的反演，能够大幅提高干旱半干旱草原区表层土壤水分的监测效率，降低成本；同时，本发明涉及的植被冠层参量均可通过波段插值等方法在星载高光谱数据中得以实现，结合热红外遥感技术，可以实现大尺度干旱半干旱草原区表层土壤水分的快速反演。
附图说明
[0066]
图1为本发明的实施流程图。
[0067]
图2为savitzky-golay滤波前后的对比示例。
[0068]
图3为土壤水分预测情况。
具体实施方式
[0069]
为使本发明的技术方案、解决的技术问题更加清楚明白，以下结合具体实施例图1，对本发明做进一步的描述。
[0070]
本发明提供了一种基于植被冠层参量和生境特征的干旱半干旱草原区表层土壤水分含量快速反演的方法，高光谱技术具有光谱分辨率高、无损等特点，能够快速获取土壤的连续光谱信息，探测土壤水分变化的细微差异，为综合反演土壤水分提供了可能。
[0071]
步骤s1、主要通过野外实地测量得到研究区的采集土壤水分数据、植被冠层反射光谱数据、叶面积指数数据、叶绿素含量数据、土壤温度数据，具体步骤如下；
[0072]
步骤s1.1、土壤水分数据采用tdr-350进行测量，该仪器可以快速简便地测量土壤水分，通过新功能的改进，能够对土壤电导率进行测量，修正土壤水分读数。测量样方土壤水分的过程中，采集样方四角以及中心点位置的土壤水分数据及土壤温度数据，随后将采集的5个数据进行平均作为该样方土壤水分数据以及土壤温度数据并加以记录，示例如下：
[0073]
测量得到的样方土壤水分数据为(％vwc)：12.40、9.60、11.80、9.30、9.70；
[0074]
样方土壤水分算数平均值计算(％vwc)：(12.40 9.60 11.80 9.30 9.70)
÷
5＝10.56；
[0075]
测量得到的样方土壤温度数据为(℃)：34.5、33.9、34.2、34.6、34.8；
[0076]
样方土壤温度算术平均值计算(℃)：(34.5 33.9 34.2 34.6 34.8)
÷
5＝34.4；
[0077]
步骤s1.2、利用svc进行测量前，先对仪器进行预热，将手簿和光谱仪预备好后，手持光谱仪探头面向太阳，保证样方内没有阴影，测量时光谱仪探头距植被冠层1m。每次测量前后用白板进行校正，每个样方测量5次重复数据，每条光谱不低于10次采样。
[0078]
步骤s1.3、叶面积指数数据采用lai-2200进行采集，在各样方点观测过程中，一般每个采集点测量一次天空光和4次植被下数据，尽量选择天空光变化较小的时段进行。同时在观测lai时要变化角度，以确保观测值的精确性，最终对每个像元内的观测值取平均值；
[0079]
在测量过程中，尽量避免直射阳光，如果天空有明显的明亮和阴暗区域，测量不均匀冠层时易出现错误；天空最亮方向下的冠层会使测量数据增大。叶面积指数计算示例如下：
[0080]
叶面积指数算数平均计算：(1.34 1.38 1.32 1.30 1.41)
÷
5＝1.35；
[0081]
步骤s1.4、叶绿素使用spad-502进行采集，其测量面积为2mm*3mm(厚度不超过1.2mm)，在测量前，先对仪器进行校准，校准后开始对样方内叶片的叶绿素含量进行测量。在测量过程中，尽量避免选择过厚的样品，如果测定有较多叶脉的叶片，可进行多次测量并取其平均值，此时记录的数据为叶片叶绿素含量。叶片叶绿素含量算数平均计算(spad)：(35.20 33.40 39.30 34.80 39.90)
÷
5＝36.5；
[0082]
步骤s2、野外实测地物光谱由于受环境、仪器和目标自身光谱特性等各种因素的影响，获得的光谱曲线中往往会含有噪声，通过光谱数据的预处理以及savitzky-golay滤波可以在一定程度上消除噪声，并突出地物光谱的某些细微差别。
[0083]
步骤s2.1、野外植被冠层反射光谱数据由svc采集，采集后的植被冠层反射光谱数据为“*.sig”格式，该格式数据可以通过svc配套的pc端软件进行查看以及编辑，但是其操作性不够强。为方便后续光谱数据的应用，将svc光谱仪采集的光谱数据“*.sig”文件转换为操作性强的“*.xls”文件，每个xls文件中包含波段信息及其对应的反射率值；
[0084]
步骤s2.2、高光谱数据的测量过程存在一定的随机误差，在数据处理的过程中利用平均的方法可以在一定程度上削弱随机误差的影响。平均的方法是指在每个样本的不同位置进行多次重复测量，每次重复采集多条光谱，随后对其数据进行平均计算，用其平均后的光谱值代表该目标的光谱属性。
[0085]
步骤s2.3、采用savitzky-golay滤波的方法对光谱数据进行多点平滑处理，也就是对窗口中奇数个点根据最小二乘方法进行拟合，再将求得的中心点的拟合值作为平滑之后该点的值。使用这种平滑方法能够较好地保持原始反射率光谱信息的同时，有效地降低和消除数据中的噪声。
[0086]
步骤s2.4、通过换算，将叶片叶绿素含量spad值转换为单位为μg/cm2，之后通过下面公式换算为冠层叶绿素含量，单位为μg/cm2。其计算公式为：
[0087]
ccc＝lcc
×
lai
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0088]
其中，ccc表示冠层叶绿素含量，单位为μg/cm2；lcc表示叶片叶绿素含量，单位为μg/cm2；lai表示叶面积指数，单位为m2/m2；如：
[0089]
ccc＝lcc
×
lai＝31.4
×
1.04＝32.66
[0090]
步骤s3、利用滤波后的反射率数据进行相关指数(归一化植被指数ndvi、增强型植被指数evi、比值植被指数sr、红边归一化植被指数rededgndvi、水分胁迫指数msi、归一化水分指数ndwi、简单比值水分指数srwi)的计算，在指数计算过程前，针对每一种指数对波段进行组合以计算该指数对应的所有波段组合，选出与土壤水分相关性最高的波段组合作为后续指数计算的依据。
[0091]
步骤s3.1、对归一化植被指数的计算，其计算公式为：
[0092][0093]
其中，ndvi表示归一化植被指数；ρ
741
表示在波长741nm处的反射率；ρ
676
表示在波长676nm处的反射率；如：
[0094]
[0095]
步骤s3.2、对归一化水分指数的计算，其计算公式为：
[0096][0097][0098][0099]
其中，ndwi表示归一化水分指数；ρ
800
表示在波长800nm处的反射率；ρ
1137
表示在波长1137nm处的反射率；ρ
500
表示在波长500nm处的反射率；ρ
741
表示在波长741nm处的反射率；ρ
1000
表示在波长1000nm处的反射率；如：
[0100][0101][0102][0103]
步骤s3.3、对简单比值水分指数的计算，其计算公式为：
[0104][0105][0106][0107]
其中，srwi表示简单比值水分指数；ρ
800
表示在波长800nm处的反射率；ρ
1137
表示在波长1137nm处的反射率；ρ
500
表示在波长500nm处的反射率；ρ
741
表示在波长741nm处的反射率；ρ
1000
表示在波长1000nm处的反射率；如：
[0108][0109][0110][0111]
步骤s3.4、对增强型植被指数的计算，其计算公式为：
[0112][0113]
其中，evi表示增强型植被指数；ρ
741
表示在波长741nm出的反射率；ρ
676
表示在波长676nm处的反射率；ρ
408
表示在波长408nm处的反射率；如：
[0114][0115]
步骤s3.5、对比值植被指数sr的计算，其计算公式为：
[0116][0117]
其中，sr表示比值植被指数；ρ
741
表示在波长741nm处的反射率；ρ
721
表示在波长721nm处的反射率；如：
[0118][0119]
步骤s3.6、对红边归一化植被指数的计算，其计算公式为：
[0120][0121]
其中，ndre表示红边归一化植被指数；ρ
741
表示在波长741nm处的反射率；ρ
690
表示在波长690nm处的反射率；如：
[0122][0123]
步骤s3.7、对水分胁迫指数的计算，其计算公式为：
[0124][0125]
其中，msi表示水分胁迫指数；ρ
1464
表示在波长1464nm处的反射率；ρ
800
表示在波长800nm处的反射率；如：
[0126][0127]
步骤s4、采用竞争性自适应重加权算法(competitive adaptive reweighted sampling，cars)对经过滤波后的反射率数据进行波段选择。
[0128]
步骤s4.1、利用matlab软件，将测得的m
×
1的土壤水分数据以及m
×
n的光谱数据作为输入数据，进行cars计算；
[0129]
cars算法是一种结合蒙特卡洛与pls模型回归系数的特征变量选择方法。cars算法中，每次通过自适应加权采样保留pls模型中回归系数绝对值权重较大的点作为新的子集，去掉权值较小的点，然后基于新的子集建立pls模型，经过多次运算，选择pls模型交互验证均方根误差(rmsecv)最小的子集中的波长作为特征波长。
[0130]
首先采用蒙特卡洛采样法，每次随机从校正集中选择一定数量(一般为80％)的样本进入建模集，剩余的20％作为预测集建立pls模型。蒙特卡洛的采样次数(n)需要提前设定。记录每一次采样过程pls模型中的回归系数的绝对值权重：
[0131][0132]
其中，|bi|表示第i个变量的回归系数绝对值；m表示每次采样中剩余的变量数；ωi表示第i个变量的回归系数绝对值权重。
[0133]
随后利用指数衰减函数强行去除回归系数绝对值权重相对较小的波长。在第i次基于mc采样建立pls模型时，根据edf得到保留的波长点的比例ri为：
[0134]ri
＝μe-ki
ꢀꢀ
(14)
[0135]
其中，μ和k实常数。
[0136]
在每次采样时，都从上次采样时的变量数中采用自适应加权采样选择数量为ri*n个的波长变量，进行pls建模，计算rmsecv。在n次采样完成之后，cars算法得到n组候选的特征波长子集，选择rmsecv最小值所对应的波长变量自己为特征波长。
[0137]
步骤s4.2、将cars算法筛选出的波段进行记录并保存。
[0138]
步骤s5、利用相关性分析的方法，分析植被生物理化参量(叶面积指数、冠层叶绿素含量、ndvi、evi、sr、ndre)、冠层含水量指标(msi、ndwi、srwi)以及土壤温度与土壤水分之间的相关性，筛选出相关性好的几种指数。步骤s5.1、将测得的土壤水分数据、叶面积指数数据、土壤温度数据，以及计算得到的冠层叶绿素含量、ndvi、ndwi、srwi、evi、sr、ndre、msi数据以样方为单位进行整合，示例如下：
[0139]
表1参数整理示例
[0140][0141]
步骤s5.2、对土壤水分和植被生物理化参量(叶面积指数、冠层叶绿素含量、ndvi、evi、sr、ndre)、冠层含水量指标(msi、ndwi、srwi)以及土壤温度进行相关性分析。
[0142][0143]
其中，σ
x，y
表示变量x，y的相关系数；μ
x
和μy分别表示变量x，y的均值(期望)，σ
x
、σy分别表示变量x，y的标准差。σ
x，y
取值范围为[-1，1]，绝对值越大，相关性越强。
[0144]
步骤s5.3、筛选出相关系数绝对值较大的指数数据进行记录并保存。
[0145]
步骤s6、将筛选出的数据(含特征波段和优选指数)进行训练集和验证集的划分；确定bp神经网络结构，构建bp神经网络模型，从而实现土壤水分的反演。
[0146]
步骤s6.1、构建含有上述筛选出的光谱波段数据以及各种指数数据的数据集，并按照7:3的比例分别划分测试集和验证集。
[0147]
步骤s6.2、确定bp神经网络结构，以筛选的波段数据和指数数据作为输入值，以土壤水分作为输出值，构建bp神经网络模型。
[0148]
步骤s6.3、利用构建的bp神经网络模型进行土壤水分的反演，将反演结果同实测值进行rmse的计算评价精度，rmse越小，反演效果越好。
[0149][0150]
综上可知，采用本发明一种基于植被冠层参量和生境特征的干旱半干旱草原区表层土壤水分含量快速反演的方法，能够充分挖掘高光谱数据中土壤水分相关的特征，结合热红外遥感技术，能够高效地进行干旱半干旱草原区表层土壤水分含量的反演。
[0151]
本发明着重强调了一种基于植被冠层参量和生境特征的干旱半干旱草原区表层土壤水分含量快速反演的方法，也可以被认为是针对特定场景下一套完整的操作流程，包括野外数据采集、光谱数据预处理、cars算法波段选择、相关指数计算及优选、构建土壤水分bp神经网络反演模型。基于本发明的上述具体实施例，本技术领域的技术人员在不脱离本发明原理的前提下，对本发明所作出的任何改进和修饰(包括更换应用场景)，皆应纳入本发明的专利保护范围。