基于分数阶微分优化光谱指数的土壤电导率高光谱估算

2019-11-13亚森江喀哈尔杨胜天尼格拉塔什甫拉提

生态学报 2019年19期

亚森江·喀哈尔,杨胜天,*,尼格拉·塔什甫拉提,张飞

1 新疆大学资源与环境科学学院, 乌鲁木齐 830046 2 新疆大学绿洲生态教育部重点实验室, 乌鲁木齐 830046

土壤盐渍化作为荒漠化的一种表现形式,会引起生态、环境、社会和经济等一系列问题,这在干旱、半干旱区表现的尤为明显[1],严重阻碍着区域生态文明建设的推进；另一方面,盐渍土作为重要的后备耕地资源,迫于人口剧增、粮食不足、环境恶化和生态破坏等的压力,开发利用盐渍土的局势也已经迫在眉睫[2]。要达到治理与利用大范围盐渍土的目的,必须及时获取有关盐渍土的可靠信息。因此,土壤盐渍化成因、空间分布、变化规律、监测方法、提高监测精度和盐渍化预警能力等成为地理学、生态学、农学等学者们的研究热点。

在土壤盐渍化动态变化监测中,具有尺度大、效率高和破坏小等特点的遥感技术近年来被广泛应用。其中,拥有更精细光谱波段信息的高光谱遥感对于提高土壤盐渍化监测精度提供了有效途径,所以国内外众多学者利用高光谱数据进行了大量土壤盐渍化监测理论与定量反演模型的研究。在1984年,Clark等[3]已经初步开始利用高光谱反射率研究地物特征,为后续研究地物高光谱特征及高光谱遥感反演提供了理论基础。Csillag等[4- 6]通过对不同盐分含量土壤在400—2500nm范围光谱进行测定,分析归纳了不同盐分含量土壤高光谱敏感吸收特征波段,指出大气对土壤盐渍化的高光谱诊断特征产生一定的影响,土壤盐分含量的反演应集中于波谱形状,而与高光谱吸收带参数关系较少。陈皓锐等[7]测定内蒙古河套灌区沙壕渠灌域沙壤土样品的高光谱曲线并进行平滑处理、相关分析和去包络线处理,分别采用偏最小二乘回归法和稳健估计法构建土壤电导率的估算模型,为大面积快速获取含盐土壤电导率和盐渍化特征提供了参考。姚远[8]基于EM38测得的盐渍土电导率数据和高光谱仪测得的盐渍土高光谱反射率数据,对高光谱反射率进行11种光谱变换后与电导率数据作相关分析,选取相关性最好的变换形式及响应波段计算5种盐分指数,并筛选最优高光谱指数,以此建立区域土壤盐渍化监测模型。彭杰等[9]利用土壤样品高光谱数据以及室内测定的盐分和电导率数据,分析耕作土壤盐分含量与电导率之间的关系,比较高光谱信息对二者的敏感性以及高光谱反演模型的精度,发现对于耕作土壤而言,电导率与盐分含量之间没有显著的相关性,因此不能利用电导率数据进行土壤盐渍化的高光谱遥感监测。

综上所述,以往对于土壤盐渍化高光谱定量估算的研究是通过全波段反射率及其对应的数学变换在一维层面上选取单个敏感波段或多个敏感波段,并利用多种回归方法建立预测模型,模型精度有待进一步提高,且基于两波段优化组合算法进行光谱指数的波段二维层面优化的应用研究相对较少。基于此,本研究以新疆艾比湖流域为研究靶区,进行野外土壤采样,室内土壤高光谱采集及理化分析等工作,尝试运用土壤原始高光谱反射率及对应的5种数学变换,对简化光谱指数(nitrogen planar domain index, NPDI) 进行波段优化计算,分析基于不同形式变换光谱的优化光谱指数与盐渍土电导率的相关性,通过变量重要性准则(variable importance in projection, VIP)筛选最优高光谱参数,并利用偏最小二乘回归(partial least square regression, PLSR)分析法建立土壤电导率高光谱定量估算模型,为土壤盐渍化动态监测及星载传感器等相关研究提供科学支持和应用参考。

1 研究区概况

图1 研究区位置与采样点分布Fig.1 Location of study area and distribution of sampling sites

地处新疆博尔塔拉境内的艾比湖流域,区域年平均气温7.7℃、年均降水量102 mm,潜在蒸散量达1447 mm,典型土壤为灰漠土、灰棕漠土及风沙土等,是典型的干旱半干旱地区盐渍化代表区[10]。艾比湖作为新疆最大的咸水湖,湖水主要依赖地表径流补给,近十几年来,由于地表径流被引入灌区、渗漏地下和消耗于地面蒸发与植物蒸腾的量增多,艾比湖入不敷出,湖泊面积严重萎缩,干涸湖底面积不断增加,湖滨荒漠化及周边区域土壤盐渍化程度不断加剧,加之阿拉山口全年8级以上大风达164d,常年侵蚀大面积裸露湖床及盐壳,使艾比湖日渐成为中国西北部沙尘暴、盐尘暴的主要策源地之一,严重影响天山北坡绿洲生态文明建设与可持续发展。

2 实验材料与数据采集

2.1 土壤样品采集与分析

土壤样品的采集时间为2016年10月15日至10月23日,针对艾比湖流域典型自然盐渍化土壤,设置5 m×5 m样方(图1),利用GPS记录每个样方位置,采用5点混合法进行样品采集,土壤样品采样深度为0—10 cm,共计57个土壤样品。土壤样品带回实验室并进行自然风干、研磨后,通过2 mm孔筛分为两部分,分别用于室内高光谱数据采集及土壤电导率和含盐量分析。其中,土壤样本的电导率和全盐量测定方法参照《土壤农业化学分析方法》,土壤电导率在土水比例为1∶5的土壤悬浊液中利用德国WTW(WissenschaftlichTechnischeWerkstätten)公司生产的inoLab® Cond 7310台式电导率测试仪进行测定,全盐量采用水浴烘干法进行测定。

2.2 高光谱测定及预处理

本研究通过美国ASD(Analytical Spectral Devices)公司生产的FieldSpec3型高光谱仪测定土壤室内高光谱数据,波段范围为350—2500 nm。高光谱波段在350—1000 nm区间的采样间隔为1.4 nm,在1000—2500 nm区间的采样间隔为2 nm,全波段范围重采样间隔为1 nm。将处理好的土壤样品分别装入黑色盛样皿(直径12 cm,深1.8 cm)中,对装满的土壤进行表面平滑处理。测定样品的光源为50 W卤素灯,测量时距土壤样品表面为50 cm,光源的天顶角为15°,待测土壤样品表面与探头的距离为5 cm。每测定两个土壤样品后进行一次白板校正,每个土壤样品重复测定5次,取5条高光谱曲线的算术平均值作为该样品的实际高光谱反射率[11]。

将ASD FieldSpec3高光谱仪测定的土壤高光谱数据通过ViewSpec Pro软件处理与导出,为了降低噪声引起的影响,去除信噪较低的边缘波段(350—399 nm及2401—2500 nm),利用Savitzky-Golay滤波方法对57个土壤样品的高光谱数据(400—2400 nm)进行平滑去噪预处理。

2.3 研究方法

本研究利用Li等[12]提出的三波段光谱指数(NPDI)归一化比值算法,其计算公式为NPDI=(RNIR/RRED-1)/[(RNIR-RRED)/(RNIR+RRED)],为进一步发挥高光谱数据的优势,将光谱指数从三波段组合降级到两波段组合,计算公式变为(Ri+Rl)/Rl,其中i与l作为土壤高光谱反射率。基于此,本文使用自主开发的《高光谱数据两波段组合软件 V1.0》(No: 2018R11S177501),将两波段优化算法应用于全波段范围进行光谱指数波段优化。通过该算法计算57个土壤样品的原始高光谱反射率R以及对应的5种光谱变换(倒数变换、对数变换、对数倒数变换、平方根变换、分数阶微分变换)在波段400—2400 nm之间进行所有可能两波段组合的优化光谱指数(NPDIs)。分数阶微分变换处理在《高光谱数据分数阶微分计算软件 V1.0》(No: 2016SR006487)中完成。为了选取模型最佳自变量参数,本文将变量投影重要性准则应用到其中,VIP值代表自变量对模型拟合的程度,自变量对因变量的解释作用相同,则自变量的VIP均接近于1。Wold[13]建议VIP值小于1的自变量对因变量的贡献较小,可以考虑剔除。因此,本研究自变量VIP的阈值为大于等于1。基于PLSR建模方法的优点,VIP技术主要用于样本较少且几个自变量间相关性较强的情形,在一定程度上弥补了传统线性回归的不足[14]。

3 分析与讨论

3.1 土壤电导率统计分析

由表1可知,建模集和验证集对应的土壤电导率最大值分别为55.70mS/cm和48.30mS/cm,最小值分别为0.20mS/cm和0.07mS/cm,均值分别为11.27mS/cm和11.18mS/cm,变异系数分别为123.64%和117.50%；研究区所有采样点土壤电导率平均值为11.32mS/cm,变异系数为122.70%,介于建模集与验证集之间,数据离散程度属于强变异性系数(C.V>100%)。

表1 土壤电导率的统计分析

3.2 土壤电导率与含盐量的相关性分析

表2是57个土壤样品的电导率和含盐量数据,从中可以看出,含盐量的变异系数也大于100%,表现为较强空间变异性。土壤电导率与含盐量的决定系数达到0.99(图2),表现为极显著相关关系,因此在本研究中土壤电导率可以很好地指示土壤含盐量。

3.3 土壤高光谱曲线特征分析

基于上述分析,本文利用盐渍土的分级标准[16]对土壤样品高光谱反射率进行归类、求平均,绘制不同电导率的土壤高光谱曲线(图3),以此大致分析土壤高光谱反射率对土壤电导率的响应。本研究区不同盐渍土高光谱反射率曲线的波动基本一致,土壤样品的高光谱反射率与其电导率未呈现明显的正负相关关系。从图3a中可知,在可见光及近红外波段(400—1500 nm)范围内高光谱曲线呈平缓上升态势,并在1400 nm、1950 nm、2350 nm左右处有3个比较明显的水分吸收谱段。总体而言,重度盐渍土高光谱反射率偏低,非盐渍土高光谱反射率偏高；除600—700 nm波段范围,不同土壤电导率的高光谱反射率曲线较易区分,且在1200—1600 nm波段区间,4种盐渍土的高光谱曲线差异性最大(图3)。

表2 土壤电导率和含盐量的统计分析

图2 土壤电导率和含盐量的相关性Fig.2 Correlation between soil conductivity and salt content

图3 土壤高光谱特征及不同盐渍化程度的土壤高光谱反射率Fig.3 Spectral reflectance of soil and soil spectral reflectance with different salinization degreeSSC: 土壤含盐量soil salt content)

3.4 优化光谱指数与土壤电导率相关性分析

图4 不同高光谱数据变换下NPDIs与土壤电导率相关性二维等势图Fig.4 Two dimensional contour map of NPDIs and soil electrical conductivity under different hyperspectral data transformationsNPDIs: 简化光谱指数nitrogen planar domain indices

高光谱数据变换Hyperspectral data transformation显著性检验Test of significance数量Amount最大值Maximum value/(r)波段组合Band combination/(inm,lnm)原数据Original data/(R)(r≥0.870,ρ<0.01)340.876(2011,1890),(2011,1891)对数Logarithmic/(lgR)(r≥0.646,ρ<0.01)350.650(2027,1881)对数倒数Logarithmic reciprocal/(1/lgR)(r≥0.870,ρ<0.01)750.876(2016,1887),(2027,1883)倒数Reciprocal/(1/R)(r≥0.870,ρ<0.01)550.880(2009,1892),(2010,1892),(2011,1891)平方根Square root/(R)(r≥0.870,ρ<0.01)410.877(2011,1891),(2010,1892)1.6阶微分1.6 order differential/(FOD)(r≥0.870,ρ<0.01)250.888(2020,1893)

3.5 土壤电导率PLSR估算模型建立及精度分析

偏最小二乘回归法是当前应用最为广泛的高光谱建模方法,为PLSR法更深入地分析数据,在建立研究区土壤电导率PLSR估算模型之前,本文利用VIP技术进一步对自变量进行筛选,图5展示了VIP准则筛选过程与对应波段组合。由图5可知,基于原始高光谱反射率及其5种变换优化的NPDI在自变量选择的情况上基本一致,筛选出的自变量数分别为17、17、28、21、23、15个,较原有的敏感波段组合数,VIP准则筛选效果明显,剔除了对模型贡献小及相对累赘的参数。

当前，针对我国经济体制改革而言，其最大障碍即为企业活力的缺乏。究其原因，主要有如下方面，其一，企业在具体体制上比较陈旧，产权不清，政企不分，造成了企业会计出现主体错位的情况。而积极构建现代企业制度，乃是奇特体制改革的基本方向。其二，企业管理不当，经营疲软。在整个企业管理架构中，人们越发注重会计人员理财管理，且已经成为企业管理的重点，因此，积极加强会计管理，同样是企业发展的基本趋势。

图5 不同高光谱数据变换下基于变量VIP值筛选最佳模型自变量Fig.5 Filtering the best model independent variables based on variable VIP valuesunder different hyperspectral data transformations

3.6 讨论

盐渍土的形成与水盐运移有密切联系,土壤水盐运移模型可以模拟区域土壤水盐运移过程,预报土壤水盐动态变化,对于改造和利用盐渍土具有重要的作用[18]。表层土壤水盐是土壤水盐运移模型的重要边界条件参数,准确的表层土壤水盐信息可以提高水盐运移模型的模拟与预测精度。通过高光谱遥感监测土壤电导率便可以及时高效地掌握表层土壤盐分的状况,这是因为对于自然土壤而言,尽管电导率和含盐量都可以反映出土壤盐渍化的程度,但相关研究[17,19]表明土壤高光谱信息对土壤电导率的响应较含盐量敏感,以土壤电导率替代含盐量进行土壤盐渍化高光谱估算研究是一种精度更高、速度更快的方法。

表4 偏最小二乘回归建模精度分析

RPD: 相对分析误差Relative prediction deviation;AIC: 最小信息准则Akaike information criterion;FOD: 分数阶微分Fractional order differentially;VIP: 重要值

传统的高光谱处理方法对土壤电导率高光谱建模时,敏感波段通常是在一维层面上以土壤电导率与高光谱反射率的相关性分析来确定,相关性越高,波段的敏感程度越高。对高光谱反射率的预处理目的就是提高土壤电导率与反射率之间的相关性,并进一步提高预测模型的精度[20-22]。高光谱丰富的波段信息为两波段优化算法提供了更多的可能组合,海量光谱数据的两波段优化算法能充分提取与土壤电导率相关性最大的波段组合,在复杂的高光谱参数中达到快速寻优的效果,深度挖掘高光谱数据从而进一步提高土壤电导率的高光谱估算精度,减少环境因素等对建模的影响[23-28]。于是本文想探讨的主要问题就是对高光谱反射率的预处理是否同样能够有利于优化光谱指数更好地估算土壤电导率。

图6 基于1.6阶微分预处理估算模型的实测与预测插值图Fig.6 Validation interpolation diagrambased onFODpreconditioning estimation modelEC: 电导率Electrical conductivity，FOD: 分数阶微分Fractional order differentially

分数阶微分在阶数上对整数阶微分的概念进行了扩展,相比整数阶微分,具有记忆性、遗传性以及非局部性,在系统控制与诊断、数字滤波、信号与图像处理等领域有着较为广泛的应用。在光谱分析领域,近期的相关研究指出,对于高光谱这类具有海量信息的高维数据源,分数阶微分也能够很好地挖掘潜在信息,弥补整数阶微分可能造成某些信息丢失的不足,并极大地扩充光谱数据预处理的方法,为高光谱研究提供一个全新的角度[29- 32]。但是,分数阶微分在单波段高光谱预处理中存在的一个问题是：它虽然可以较好地增加敏感波段的数量,却不能有效地提高相关性。因此,本文将分数阶微分与传统光谱变换方法进行对比研究,探究分数阶微分在优化光谱指数中的效用,发现分数阶微分预处理方法可以有效地提高优化光谱指数相关性的极值。

为比较一维与二维高光谱数据处理效果的差异,绘制了高光谱反射率及其不同数学变换后在一维层面上与土壤电导率进行相关性分析得出的最大相关系数图。由图7可知,传统光谱变换的相关性提升效果不明显,平方根变换的相关性最高为0.633,对数倒数、倒数变换的相关性明显降低,最低为-0.566,而在优化光谱指数中表现较差的对数变换没有太大的降低；在分数阶微分变换中,多数都有提升效果,最高的为一阶微分变换,达到-0.673,这与优化光谱指数的结果形成了鲜明的差异,最低的为二阶微分,仅达-0.436,而在二维层面表现最好的1.6阶微分在一维中没有起到相同的作用；从相关性质方面来看,原数据、对数和平方根变换呈正相关,对数倒数、倒数变换呈负相关,这与优化光谱指数的结果相异,分数阶微分的相关性质以一阶微分为界,0.2-0.8阶为正相关,1.0—2.0阶为负相关。综上所述,不管在一维层面上还是二维层面上,合适的高光谱数据预处理方法都在一定程度上对提高相关性有所帮助,分数阶微分总体上优于传统预处理方法,而且优化光谱指数对比于传统的高光谱处理方法来说具有明显的优势[33-34]。

图7 不同预处理下单波段高光谱数据与电导率的最大相关系数Fig.7 Maximum correlation coefficient of single-band hyperspectral data and conductivity under different pretreatments

与已有研究[35-42]相比,本文的研究特色在于：将高光谱数据分析从传统的一维层面上升至二维层面,并结合分数阶微分预处理进行光谱指数波段优选,用于建立土壤电导率估算模型,以提高土壤电导率反演精度,为土壤盐渍化相关研究提供一种新的思路和方法。而且本文得到的优化光谱参数可为快速准确寻求卫星传感器中监测干旱、半干旱地区土壤电导率的最佳波段提供依据,此外,波段的优化也可以为设计特定波段的主动传感器提供理论基础,进一步减少高光谱海量数据处理的工作量,为实现土壤盐分信息的高效监测服务。最后,研究区虽然属于典型的干旱、半干旱区,但是干旱、半干旱区乃至中国具有区域异质性[43],这就会不可避免地导致本研究确定的土壤电导率反演优化光谱参数仍有一定的地域局限性,因此,比较光谱指数优化算法在不同地区的最佳参数并找出普适性高光谱参数,将是值得研究的方向。