1 引言

土壤含水量(SWC)是土壤理化特性的一个重要指标,可以用质量含水量、体积含水量等方式表示(文中后续含水量都是以含水量质量百分比来表示)^[1-2]。不论是对水资源管理还是对土壤墒情进行预报,土壤水分的遥感监测都有着至关重要的意义。因此,及时监测和掌握土壤水分状况具有极其重要的实际应用价值^[3]。

传统的土壤含水量测定方法--取土烘干法、中子水分仪法等均是以点测量、实验室测定为基础,精度高但范围有限、工作量大、效率低,难以满足土壤含水量监测的实际需要^[4-5]。近年来随着遥感高光谱技术的发展,利用高光谱遥感监测土壤水分技术能够快速获取土壤反射光谱信息,实现实时高效的大范围精准的土壤含水量信息提取,弥补了传统方法的不足,因此成为研究的热点^[6-7]。利用遥感技术进行土壤含水量测定的基础是确定土壤含水量与土壤反射光谱间的关系,并进行反演。这方面国内外学者已有大量研究,研究表明当土壤含水量在一定范围内时,土壤光谱反射率随土壤水分的增加而降低^[8-11]。

目前利用土壤光谱特性反演土壤含水量的方法主要分为两类,一类是基于经验的统计模型法,主要包括多元逐步回归、偏最小二乘法回归、反向传播(BP)神经网络;另一类是基于机理信息的机理模型法^[12-15]。但在机理模型法中,由于不同含水量造成的土壤反射率的定量变化机理非常复杂,目前还未获得反演的最佳参数值,其反演效果及结果的适应性有限,现在多采用基于经验的统计模型。其中BP神经网络因具有较强的自主学习能力,能较好地处理土壤含水量反演中的非线性情况^[16],已有研究表明该模型土壤水分反演值与实测值的误差在9%左右,反演效果优于其他大多数模型^[17]。

在反演中,原始土壤光谱反射率数据包含噪音及大量的冗余信息,直接用于土壤含水量的反演计算复杂且效果不好。因此,绝大多数情况下都需对原始土壤光谱反射数据进行光谱变换提取特征波段或参数获得反演因子^[18-20]。传统对光谱数据进行预处理的方法主要包括:倒数、对数、一阶微分(FD)等^[21-23],但这些方法降噪功能低,无法很好地处理数据背景、噪声等,这直接影响了后续反演精度。

土壤含水量反演精度取决于反演模型和反演因子的选择。在反演模型确定时,反演因子的质量直接决定了反演效果的好坏。为得到好的反演因子以提高反演精度,本文引入谐波分析(HA)理论对传统方法处理后的光谱数据做时频空间转换^[24-25]。通过谐波分解抑制或消除地物背景噪声,达到数据压缩的效果^[26-28]。本文利用谐波分解后获取的最佳谐波分量作为反演因子,利用反演效果更佳的BP神经网络,构建土壤含水量反演模型。

综上,从反演因子的构建着手研究,在利用光谱FD、小波包变换(WPT)等光谱预处理方法的基础上,引进了HA,对原始光谱(OS)数据进行光谱变换处理提取敏感特征波段、光谱吸收参数及谐波分解参数。再分别进行主成分分析(PCA)得到三种不同的反演因子,最终基于三种反演因子建立BP神经网络模型来反演土壤含水量,对三种模型的反演精度进行比较,最终得到最适用于土壤含水量反演的反演因子获取方法。

2 材料与方法

2.1 研究区概况

实验区为横山县(北纬37°22'~38°14',东经109°14'~110°20'),地处陕北黄土梁峁丘陵区北端,实验区内总体地势西南高东北低,平均海拔1250 m,年均气温8.6 ℃,年降水量399 mm,属温带干旱、半干旱大陆性季风气候,受北部西伯利亚冷高压气候的长期影响,出露的地层岩石逐渐被剥蚀、搬运、沉积而使境内呈现出“三沙七丘”的地貌特征,形成了包含黄绵土、棉沙土、风沙土等众多类型的土壤,不同类型的土壤含水量不同,相同类型土壤因地域差异含水量也不同。横山县河道属黄河流域,主要支流有无定河、芦河、大理河等,靠近河流干道的土壤含水量较多。正是由于这些地理因素造成该实验区不但土壤类型丰富,且含水量也有较大差异,这对于含水量反演研究非常有利。选择横山县为研究区,可以避免由于单一土壤、单一含水状况造成的偶然性和局限性,使结果更具推广性和普遍性。

2.2 样品采集与制备

在采集土样时,选择地势平坦、土壤裸露区作为采样区。考虑到不同的土地利用类型以及土壤类型,每个样区内选择4~5个具有代表性的测点,每个测点采集一个表层土(约5 cm)土样,共采集84个土壤样品,如图1所示。将采集的样品带回实验室进行封装,制备成待测土壤样品。

图 1. 土壤样本采集区域

Fig. 1. The area of collecting soil samples

下载图片 查看所有图片

2.3 土壤反射光谱及含水量测定

2.3.1 土壤光谱测定

将采集的84个土壤样品磨碎,然后用直径为2 mm的孔径筛进行过筛处理,作为待测光谱反射率的土壤样品,在室内采用ASD Field Spec FR地物光谱仪对制备好的2 mm粒径土壤样品进行光谱测定,光谱仪测量波长范围为350~2500 nm,其中350~1000 nm之间的光谱分辨率为3 nm,1000~2500 nm之间的光谱分辨率为10 nm。室内光谱测定条件包括:光源为50 W的卤素灯,光源离待测土样约为50 cm。采用3°视场角的探头,探头直径约为2 cm,探头垂直对准样品,距离待测土样约为20 cm。测定前先用75 cm×75 cm的白色参考板定标,使待测土样充满整个视场,每条光谱曲线的扫描时间设定为5 s,每个土壤样品的光谱测量4次,最后取4次反射率平均值作为土壤样品的光谱反射率,土壤样品PCA曲线如图2所示。

2.3.2 土壤含水量测定

在室内采用称重法(烘干法)测定土壤样品中的含水量,该方法(重量含水量法)为国际上认可的测定土壤含水量标准方法^[29],测定过程为:首先将采集的土壤样品都分别装入铝盒,然后用0.1 g精度的天平称取铝盒和土壤的总重量M₁,然后在105 ℃烤箱内将土壤样品烘6~8 h至恒重,用天平测量烘干后的铝盒及土壤样品的总重量M₂,再单独测量烘干后每个铝盒的重量M₃,最后计算土壤样品中的含水量为

VSWC=M1-M2M2-M3×100%,(1)

完成土壤样品的制备。

图2为测定的其中7条不同含水量土壤光谱曲线,这7条反射光谱曲线波形基本相似,随着土壤含水量数值的增加,光谱反射率逐渐下降。

图 2. 土壤含水量OS反射率

Fig. 2. OS reflectance of soil water content

下载图片 查看所有图片

3 光谱数据预处理与土壤含水量反演流程

为消除仪器的系统误差及光谱中部分线性背景和噪声,对高频信息进行更精细的分解与重构,去除高频信息的噪声,本研究首先对OS数据进行预处理,采用的方法包括:五点加权平均、光谱FD、WPT和HA。最终为相关性分析做准备。预处理流程包括:五点加权平均、光谱FD、WPT、HA、PCA。

3.1 五点加权平均

为消除随机误差,首先采用五点加权法处理PCA数据,将具有5个光谱数据点的窗口在每个样本的光谱数据中移动,窗口中间数据点值用其前面和后面的两个数据平滑,离中间的越远,分配的权重就越小,窗口中间点的值即为这五点的加权平均值,如此计算直到最后即为五点加权平均^[30]。计算公式如下:

S=(sn-2/4+sn-1/2+sn+sn+1/2+sn+2/4)/2.5,(2)

式中S为光谱的五点加权平均值,s_n为平滑前光谱中点值。

3.2 光谱FD

为消除目标光谱中部分线性以及接近线性的背景、噪声等影响。在五点加权后再进行光谱FD处理,得到FD光谱数据。计算如下式所示:

R'(λi)=R(λi+1)-R(λi-1)/(λi+1-λi-1),(3)

式中λ_i-1、λ_i、λ_i+1为相邻波长,R'(λ_i)为波长λ_i对应的FD光谱值,R(λ_i)为波长λ_i对应的光谱反射率。

3.3 WPT

为去除高频信息的噪声,采用WPT对高频信息进行更加精细的分解与重构。WPT是在小波变换的基础上引入^[31-32]的。WPT在分解和重构高频信息方面优于小波变换,信息处理结果无冗余也无疏漏,因此更有利于光谱去噪及原始信息保留。本研究按以下4个步骤实现对光谱数据的分解与重构:1) 光谱的小波包分析,研究中所用的小波母函数为Db10^[32],对含噪光谱进行多层WPT分解;2) 确定最佳小波包基,根据最小代价原理计算出WPT分解的最优小波包基;3) 小波包系数阈值化,选择连续性好的软阈值对小波包系数进行量化;4) WPT光谱信息重构,根据最优小波包基和量化后的最优小波包分解系数进行光谱重构,得到WPT降噪处理后的光谱。

根据以上步骤进行WPT去噪的处理结果如图3所示。图3为84组光谱数据中某一条光谱曲线去噪前后的对比,由于原始曲线在红框标注的波段区间内含噪声较多,而WPT恰好能达到在此波段区间的平滑作用。然后对WPT数据进行FD处理得到WPT-FD数据。

图 3. 小波包(a)去噪前和(b)去噪后的光谱曲线

Fig. 3. Spectral curves (a) before and (b) after wavelet packet denoising

下载图片 查看所有图片

3.4 HA原理

经过上面的数据预处理后依然存在大量噪声、冗余等问题,针对此问题,本研究引入HA对处理后的光谱数据作为序列信号进行时频空间转换。HA最早由Jakubauskas等^[24]提出,主要应用于电力系统谐波监测等研究领域,光谱数据经过谐波分解得到一系列的谐波能量特征参数。其谐波理论是通过正(余)弦波(谐波)相叠加的形式来表示时间序列f(t),即任何一条关于时间t的时间序列函数f(t)都能以若干个正(余)弦和的波形叠加表示^[33]。光谱仪采集的每个样本的光谱数据可表达为一条连续的曲线。利用HA处理光谱数据时,可将由N个波段组成的光谱曲线看作周期为N的函数,HA分解就是将每个样本的光谱曲线表达为一系列由谐波余项(A₀/2)、振幅(A_h、B_h、C_h)和相位(φ_h)等能量谱特征分量组成的正(余)弦波叠加之和。若将由N个波段组成的一组光谱表示为V(k)=(v₁,v₂,…,v_s,…,v_N),每个波段的光谱值记为v_k,k为波段序列号(k=1,2,…,N)。N为波段总数(周期),h次HA得到的谐波分解展开式为

V(k)=A02+∑h=1∞[Ahcos(2πhk/N)+Bhsin(2πhk/N)]=A02+∑h=1∞Chsin(2πhk/N+φh),(4)

V(k)的h次谐波分解的谐波特征分量计算为

A0/2=1N∑k=1Nvk,(5)Ah=2N∑k=1Nvkcos(2πhk/N),(6)Bh=2N∑k=1Nvksin(2πhk/N),(7)Ch=(Ah2+Bh2),(8)φh=tan-1(Ah/Bh),(9)

式中h(h=1,2,3,…) 为谐波分解次数;A₀/2为谐波余项;C_hsin(2πhk/N+φ_h)为第h次谐波分量,当h=1时表示基波分量;A_h、B_h、C_h和φ_h分别为第h次谐波分解的余弦振幅、正弦振幅、谐波分量振幅和谐波分量的相位。其中,A₀/2、C_h和φ_h反映了光谱在各个波段的能量均值、不同波段的能量变化及能量出现幅值的波段位置。

低次谐波包含了波谱的主要能量特征,高次谐波一般夹杂着噪声信息,而振幅和相位分别承载着地物波段的能量和辐射峰值位置信息,反映着波谱的局部特征信息,因此谐波分解既能抑制或消除地物背景噪声,又能以低次谐波分量突出地物的光谱特征,达到数据压缩的效果。

3.5 PCA

PCA是一种在数据特征提取、压缩和降维等方面广泛应用的统计分析方法,利用PCA方法对上述提取的数据特征:FD数据、WPT-FD数据、WPT-FD-HA数据进行变换,在尽可能多地保留原始数据信息的前提下,消除数据间的相关性,得到反演因子。一般主成分反演中选取PCA后累计方差贡献率大于90%的主成分分量作为反演的参量^[27]。

3.6 基于BP神经网络的土壤含水量反演流程

对上述预处理后获得的反演因子分别构建BP土壤含水量反演模型,其具体流程如图4所示,图中矩形框里面代表实验流程走向,无框的描述代表实验操作,阴影箭头代表向左/右参与的操作。

图 4. 土壤含水量反演流程图

Fig. 4. Inversion flow chart of soil water content

下载图片 查看所有图片

4 结果与分析

4.1 特征波段提取

将利用OS数据得到的FD光谱、WPT数据及WPT-FD数据与土壤含水量进行相关性分析,分析结果如图5所示。

图 5. 不同数据与土壤含水量的相关系数

Fig. 5. Coefficients between different data and soil water contents

下载图片 查看所有图片

由图5可得,OS数据和WPT数据分别与土壤含水量整体呈负相关,FD数据、WPT-FD数据与土壤含水量的相关性都较高且正负交替,相关性较高的波段主要集中在近红外波段。其中光谱数据经过WPT重构后,WPT-FD数据与土壤含水量的相关性明显提高。在特征波段挑选过程中,为了减少进行土壤含水量反演的波段总数,达到数据压缩的效果,通过分别计算FD数据和WPT-FD数据各个波段与土壤含水量的相关系数,最终确定特征波段挑选原则为:1) 相关系数R大于0.5;2) 相关系数依次从大到小选择;3) 波段数量适中,两种数据波段总数相同,起到控制变量作用。基于以上挑选原则最后在FD数据和WPT-FD数据中分别挑选了同波段区间内的相同的180个特征波段。

4.2 谐波分解

对提取的WPT-FD数据利用(4)~(9)式进行谐波分解,获得谐波能量谱特征分量A₀/2、A_h、B_h、C_h和φ_h,计算其与土壤含水量的相关性系数,分解中波段总数为180,考虑到正(余)弦函数的周期性,因此确定谐波分解次数为180。图6为不同谐波能量谱特征分量与土壤含水量的相关系数图。

图 6. 不同能量谱特征分量与土壤含水量的相关系数

Fig. 6. Correlation coefficient between different energy spectrum characteristic parameters and soil water contents

下载图片 查看所有图片

由图6可得,前后10次谐波分量与土壤含水量相关性较高,并且谐波能量谱特征分量A_h、B_h、C_h和φ_h呈轴对称,因此,最终选取了11个相关系数大于0.6的谐波分量,分别为:A₀/2、A₁、A₂、A₄、B₁、B₂、B₃、C₁、C₂、C₃和φ₁。

4.3 BP反演因子构建

本研究在利用BP神经网络进行反演时,为了减少输入层神经网络个数,改善网络拓扑结构与模型反演效果,分别选取了PCA结果的前5个主成分分量作为土壤含水量BP模型输入因子,表1~表3为PCA结果。

比较表1和表2中累计方差贡献率,发现同样的波段区间WPT-FD数据的PCA结果要优于FD数据,WPT-FD数据前5个主成分累计方差贡献率达到96.085%,完全达到输入因子的构建要求。

表3为谐波能量谱特征分量的PCA结果,可以看到第一个主成分的方差贡献率就达到93.900%,而前5个主成分方差贡献率的累积方差贡献率为99.013%,基本上完全包含了原有的谐波能量谱特征分量,数据降维方面效果更好。基于以上主成分分析结果,共构建3种BP模型输入因子,分别为:FD-PCA、WPT-FD-PCA和WPT-FD-HA-PCA。

4.4 土壤含水量BP模型的建立与反演结果分析

4.4.1 土壤含水量BP模型构建

本研究采用的BP模型拓扑结构为5-5-1,即输入层节点数为5,隐含层个数为5,一个输出结果层,设置网络训练迭代次数2000,学习速率0.01及附加动量因子0.9,学习误差0.001。模型反演结果的均方根误差(RMSE)为

VRMS=∑i=1n(yi-y˙i)2/n。(10)

式中y_i为实测值, y˙i为反演值,n是样本数。

4.4.2 土壤含水量BP反演结果与分析

从84组样本中随机抽取45组作为训练样本的学习目标,其余39组样本数据作为测试数据,构成训练集和测试集数据,基于FD-PCA、WPT-FD-PCA和WPT-FD-HA-PCA三种反演因子构建土壤含水量BP反演模型,最终共建立三个BP反演模型分别为:FD-PCA-BP、WPT-FD-PCA-BP和WPT-FD-HA-PCA-BP,表4分别是这三个模型的反演精度。

根据表4中决定性系数R²,比较三个模型的反演精度可得,利用WPT处理后的数据构建反演因子WPT-FD-PCA,其反演结果比FD-PCA-BP模型的结果有了较大的提高,WPT-FD-PCA-BP模型决定性系数R²达到0.9025,V_RMS为2.3648。三个模型中反演效果最好的为WPT-FD-HA-PCA-BP,决定性系数R²=0.9599,V_RMS为1.6670,明显优于其他两个模型的反演结果。与大量学者通过不同的方法来反演土壤含水量相比,其反演精度明显优于已有大量文献中较好的结果。

图7为土壤含水量实测值与三个模型反演值的拟合效果比较,从图7中可以看出,WPT-FD-HA-PCA反演因子构建的BP反演模型预测结果最好,反演值最接近于实测值。图8为模型反演值与实测值之间的相对误差。在图8中,样本号24对应的FD-PCA-BP和WPT-FD-PCA-BP模型土壤含水量反演误差极大,说明两种反演方法的效果都不如WPT-FD-HA-PCA-BP模型。样本号26对应的FD-PCA-BP模型相对于其他两种模型来说,相对误差最大,说明FD-PCA-BP模型相比其他两种模型反演效果较差。可以得出,WPT-FD-HA-PCA-BP模型的反演结果更接近于实测值,其相对误差最小;而FD-PCA-BP反演模型预测效果最差。

图 7. 土壤含水量实测值与反演值的比较

Fig. 7. Comparison of measured and inversed values of soil water content

下载图片 查看所有图片

图 8. 土壤含水量实测值与反演值的相对误差

Fig. 8. Relative error of measured and inversed values of soil water content

下载图片 查看所有图片

三个BP反演模型中,当原始数据经过WPT变换后,去掉了一部分数据噪声达到了一定的平滑效果,因此,WPT-FD-PCA-BP模型精度要高于FD-PCA-BP。而进一步通过HA不仅去掉了高频部分的噪声,并且完成了数据压缩,减小了数据冗余,所建立的WPT-FD-HA-PCA-BP模型精度达到最高,更有利于土壤含水量的BP反演。同时也研究了直接基于HA的FD-PCA-HA-BP反演模型,其模型精度低于WPT-FD-HA-PCA-BP模型,因此WPT在光谱去噪方面是比较重要的。比较FD-PCA-BP、WPT-FD-PCA-BP和WPT-FD-HA-PCA-BP这三个模型的反演结果,可以得出结论:WPT和HA能有效地去除光谱中隐含的噪声,并能够较好地压缩OS数据,构建的BP反演因子更有利于BP模型反演土壤含水量。

5 结论

在对光谱数据进行五点加权、FD和WPT光谱数据预处理的基础上,引进HA对处理后光谱数据进行谐波分解,分别基于PCA变换后的主成分构建反演因子,最终获得三种反演因子构建方法,通过比较三种反演因子利用BP模型进行的土壤含水量反演的精度,得出以下结论。

1) WPT重构有效提高了土壤光谱与含水量之间的相关性,而HA进一步去除FD光谱中的噪声,并起到了平滑与压缩数据的效果。经HA后其谐波分量与土壤含水量的相关性有了明显的提高,说明此方法在高光谱数据反演中是有效的。

2) 基于PCA,建立了三个BP反演模型:FD-PCA-BP、WPT-FD-PCA-BP和WPT-FD-HA-PCA-BP,比较三个模型的反演精度,用同样的BP模型进行反演,发现WPT-FD-HA-PCA-BP模型的反演效果最佳,其模型决定性系数R²=0.9599,V_RMSE=1.6670。此研究结果优于大多数关于土壤含水量的BP反演精度。因此,基于HA的BP反演因子的构建能有效地提高土壤含水量的反演精度。

本研究选择的研究区包含了大多数种类的土壤,且包含多组土壤含水量反射率光谱数据,避免了实验结果的局限性和单一性,其运用的HA和WPT对数据具有较好的降维效果,因此实验结果具有较好的适用性。该研究结果对确定反演土壤含水量的模型最佳反演因子具有较高的参考价值,此方法对利用高光谱遥感数据反演其他信息时反演因子的构建与选取具有借鉴意义。