多尺度土壤质地与光谱空间非平稳性关系探究

郑嵛珍，陈奕云，陈 敏，吴子豪，蒋江俊男

（1.武汉大学资源与环境科学学院/自然资源部数字制图与国土信息应用重点实验室/地理信息系统教育部重点实验室，武汉 430079；2.广州市城市规划勘测设计研究院，广州 510060）

土壤与人类的生存和发展有着密切联系，在解决全球现实问题，尤其是粮食安全、水安全、气候变化、生物多样性等方面，土壤居于中心地位［1］。在全球可持续发展目标（Sustainable development goals，SDGs）中，有多个目标与土壤密切相关，人类社会已经认识到，科学合理地利用和管理土壤信息对人类自身的生存和可持续发展至关重要［2］。土壤质地是土壤重要的物理属性，也是陆面过程模拟的关键参数，对土壤结构、孔隙状况、保肥性、保水性、耕性等均有显著影响［3］。有效地获取、存储、表达、传输与分析土壤质地信息，理解和表征土壤质地的空间变化信息，能为粮食安全、生态文明建设、乡村振兴、精准扶贫等国计民生提供重要决策支持［4，5］。

传统土壤质地信息获取具有周期长、成本高、过程复杂等显著性缺点，难以进行大范围、高覆盖度的重复调查。基于高光谱技术，可以快速获取土壤属性性质及其空间分布特征。高光谱在土壤研究中的应用历史可以追溯到19世纪20—30年代［6］。随后，人们开始利用反射光谱研究土壤含水量、土壤组分、粒径及有机质含量。自20世纪60年代起，研究者基于可见光-近红外（Vis-NIR）光谱较为成功地反演了一系列土壤理化属性。如今高光谱技术被广泛地用于反演土壤中有机质含量、含水量、重金属及土壤质地等方面的研究［7-13］。

目前，高光谱技术在土壤质地领域研究热点多集中在光谱预测模型方面。学者们大多采用多元线性逐步回归法、偏最小二乘回归法、随机森林算法、BP神经网络、支持向量机等一系列线性和非线性方法进行土壤质地高光谱预测。张娜等［14］利用一元线性回归、多元逐步回归及BP神经网络3种方法建立土壤光谱反射率与土壤沙粒、粉粒之间的预测模型并对比了预测精度。王德彩等［15］应用Vis-NIR光谱结合BP神经网络预测土壤质地组分颗粒黏粒和沙粒，结果显示BP人工神经网络建立的模型精度更优。Hobley等［16］结合偏最小二乘回归模型和随机森林模型预测土壤质地，发现偏最小二乘回归（PLSR）模型对黏粒和沙粒预测效果最好，而RF模型对粉粒组分预测效果较好。然而，这些方法忽略了土壤质地与高光谱关系的空间异质性，使得残差存在空间结构。

地理加权回归（GWR）是一种常用于估计目标变量的空间分布和探究解释变量和响应变量之间非平稳关系的空间分析方法。由于GWR考虑了关系的空间异质性，能反映局部情况且具有更高准确性，已被广泛用于土壤属性空间预测［17］。在此基础上，Fotheringham等［18］于2017年提出了多尺度地理加权回归（MGWR）方法。MGWR不仅继承了GWR反映空间对象局部效应的优点，还考虑了不同影响因子的差异化作用尺度，能更灵活地识别不同空间尺度下各自变量影响因变量的地理过程［19］。沈体雁等［20］应用MGWR识别出不同尺度下北京市2011—2017年二手住宅交易价格的空间分异特征。MGWR模型还能应用于探究土壤属性和光谱关系空间非平稳性的问题。乔磊等［21］应用MGWR方法，很好地实现土壤有机质的空间预测，得到各环境协变量对有机质影响效应的空间变化，这种影响效应在不同环境协变量中有着不一样的空间尺度。因此，使用MGWR可以更好地探究不同空间尺度下光谱与土壤质地之间的非平稳性关系。

云南省杞麓湖是典型的高原湖泊，湖滨平原大部分被开垦为农田，农田分布细碎化特征导致农民土地管理方式多样化，进而使土壤质地分布可能具有异质性特征［22］。针对受人类活动影响较大的农耕区，土壤质地和土壤光谱关系是否具有空间非平稳特征更值得被探索。因此，选取云南省杞麓湖流域农田为研究区，以光谱的不同潜变量作为协变量，土壤质地的3种颗粒含量作为因变量，对比PLSR、GWR、MGWR 3种方法下土壤质地拟合模型效果，并应用MGWR探究多尺度效应影响下土壤质地和高光谱关系的空间非平稳特征，有益于加深对流域土壤属性空间分布的认识和理解，以期为相似流域开展空间异质性研究提供普适性参考。

1 材料与方法

1.1 研究区概况

杞麓湖流域（102°33′48″—102°52′36″E，24°4′36″—24°14′2″N）位于云南省中部，隶属云南省玉溪市通海县，南北衔接曲江干流和江川星云湖，东西邻近华宁龙洞河与玉溪大河（曲江上游段）。该流域属亚热带季风气候，干湿季分明，多年平均降水量约899 mm，年平均气温约15.6℃，年日照时间约2 300 h。流域地形为向南突出的新月形断陷盆地，平均海拔约1 800 m。湖盆土壤类型主要有山地红壤、紫色土和水稻土等。沿湖平原是重要的农耕区，土地利用方式主要是农业用地，水稻、小麦、烤烟、油料是主要农作物，蔬菜产业尤为发达。流域农田较为细碎，分布零星，70%以上现有耕地田块面积小。由于受耕作、施肥等人类活动影响，土壤质地与光谱之间的关系可能具有非平稳性。210个土壤样本均采集自杞麓湖流域周边农田（图1）。

图1 土壤样点空间分布

1.2 土壤机械组成测定和光谱预处理

采集土壤分成2份，分别用于土壤质地室内分析和室内光谱测量。样本测定前先经自然风干、过筛、研磨等操作后，过2 mm筛孔，经H2O2-HCl-（NaPO3）6法预处理后，形成最终上机样本。土壤质地测定方法采用激光粒度仪法，试验仪器为马尔文激光粒度仪（Mastersizer 3000），量程达0.01～3 500.00μm。土壤颗粒采用美国农业部土壤质地分类标准（USDA）进行分级：黏粒（＜0.002 mm）、粉粒（0.002～0.050 mm）和沙粒（0.050～2.000 mm）。

研磨土样过20目筛孔后进行室内高光谱数据测量，试验条件为避光室内环境。使用ASDFieldSpec 3地物光谱仪测定反射光谱，获取的地物波长范围为350～2 500 nm，输出波段数可达2 151 nm。每次测量前及时进行标准白板校正，以保证测量效果稳定。光谱分析前去除边缘噪声，保留信噪比较高的400～2 350 nm波段。光谱预处理步骤包含Savitzky-Golay平滑（7点平滑数）、一阶微分变换、均值中心化等，以上操作均在MATLAB中进行。

1.3 模型建立与评价

1.3.1 偏最小二乘回归（PLSR） 偏最小二乘回归由Wold等［23］首次提出，最早应用于化学工程领域，目前常用于定量光谱分析领域［16］。PLSR集成了多元线性回归、主成分回归和典型相关分析的优点，能有效解决光谱多波段间存在的多重共线性问题。PLSR的计算公式如下。

式中，y代表目标因变量；x代表土壤光谱参数，β（ii=1，2，……n）表示xi的回归系数。使用PLSR建立多自变量（光谱）与因变量（土壤属性）之间的回归关系，PLSR因子即潜在变量（LVs）的最佳数量通常使用交叉验证法来确定。模型建立采用Matlab 2014软件和PLS_toolbox工具箱实现。

1.3.2 偏最小二乘-地理加权回归（PLS-GWR） 地理加权回归是探索潜在空间非平稳关系最为广泛应用的方法之一［24-26］。GWR模型是对普通线性回归模型的扩展，将数据的地理位置嵌入到回归参数中。与普通线性回归相比，GWR原理简单、操作性强，还可提高模型拟合度和降低残差空间自相关性。与传统全局回归模型相比，GWR不会对模型中的每个关系产生“平均”全局估计，而是允许自变量和因变量之间的关系因空间而异。PLS-GWR即偏最小二乘法与地理加权回归法的结合，是将偏最小二乘方法中获取的多个潜变量作为地理加权回归的自变量，考虑了变量的空间位置，从而建立起自变量和因变量之间的回归关系。

假设共有n个观测样点，每个观测点i∈（1，2，……n）所属的空间位置为(ui,vi)，GWR的线性回归方程如下。

式中，(ui,vi)指第i个样本的空间坐标；xij指i位置的第j个自变量；yi是i位置上的因变量；βj(ui,vi)是i位置(ui,vi)上的第j个自变量的回归系数；εi是随机误差项。

空间权重矩阵和带宽的确定是地理加权回归建模计算和分析的关键，GWR分析中最常使用高斯函数（Gaussian）和双平方函数（Bisquare）2种核函数确定空间权重矩阵。带宽的确定方法仍以交叉验证准则（CV）或修正赤池信息量准则（AICc）为主。模型各项参数选择如下：建模采用双平方自适应核函数，带宽的选取采用默认黄金分割搜索法，以修正赤池信息量准则（AICc）作为信息评价准则。经典GWR模型中预设所有建模的过程皆在相同的空间尺度中操作，即所有自变量共用相同的有效带宽和空间权重矩阵。

1.3.3 偏最小二乘-多尺度地理加权回归（PLS-MG⁃WR） 通过放宽模型中所有空间变化过程在相同空间尺度上运作的假设，Fotheringham等［18］提出多重尺度地理加权回归（MGWR），这是在经典GWR基础上优化的一种改善有效带宽选择的地理空间回归方法。不同于GWR模型中各自变量共用一个带宽，MGWR模型在拟合过程中针对每个自变量选择了差异性带宽，这意味着每个空间位置上针对各变量都有特定的空间权重矩阵。PLS-MGWR是偏最小二乘方法与多尺度地理加权回归方法的结合，是将偏最小二乘方法中获取的多个潜变量作为多尺度地理加权回归的自变量，可以看作对PLS-GWR的改进。MGWR模型的表达式如下。

式中，βbwj指第j个自变量的回归系数；bwj是计算第j个自变量回归系数时使用的带宽；εi是随机误差项。

MGWR和经典GWR选择的核函数和带宽选择准则相同，分别为自适应高斯函数和修正赤池信息量准则（AICc）。MGWR模型预设所有建模的过程皆在不同的空间尺度中操作，即针对每个自变量确定差异性带宽和空间权重矩阵。后向拟合算法（Back-fitting）是校准MGWR模型的核心算法，其基本思想是模型中每一个加性项都使用单个光滑函数来估计，而每一个加性项都可以解释因变量和自变量之间的因变关系。遵循广义加性模型（GAM）的逻辑，MGWR模型还可表示如下。

式中，MGWR中的fj是指第j个加性项fj；ε是残差项。

多尺度地理加权回归主要基于线性可加模型和后向拟合法进行参数估计迭代，其主要做法是以平滑因子或者局部变量参数的地理加权回归估计值作为初始值，结合后向拟合法和迭代过程，来实现对平滑因子或变量参数的估计［18］（图2）。以经典GWR各参数和结果作为MGWR初始值，得到初始化残差ε（y的实测值与初始估计值之差）。将初始化残差ε与第一加性项f0相加得到新值，对第一个自变量x0和新值进行经典GWR回归，得到有效带宽bw0和新残差，以此更新下一个加性项f（jβbwjx）j，随后新残差加第二个加性项f1与变量x1再进行GWR回归，得到新带宽bw1和新残差，重复以上过程，直到与最后一个变量xm相关的局部参数都已估算完成即可。

图2 MGWR中后向拟合算法原理

1.3.4 模型评价指标 根据土壤反射光谱和黏粒、粉粒、沙粒3种颗粒之间的关系，建立PLSR、PLS-GWR、PLS-MGWR模型，随后进行模型质量评价。评价采用以下指标，包含决定系数（R2），平均绝对误差（MAE），均方根误差（RMSE），全局莫兰指数（Moran′sI）。

式中，n是参与建模的土壤样本数量是模型中黏粒（粉粒或沙粒）的拟合值；yi是实验室黏粒（粉粒或沙粒）测量值是测量值均值。决定系数越大，平均绝对误差和均方根误差越小，表明模型相对越稳定，相对拟合程度越好。莫兰指数介于［-1，1］，绝对值越大则表明空间自相关性越强；反之，空间自相关性越弱。使用Moran’sI来检测模型残差的空间结构［5］，若结果显著，则说明残差具有空间结构，模型的参数估计结果不可靠，反之，则说明残差空间分布随机，模型结果可靠。

2 结果与分析

2.1 杞麓湖流域农田土壤质地空间分布状况

本研究土壤质地分类标准为美国制，其中土壤质地等边三角形三顶点分别代表沙粒、黏粒、粉粒的0或100%的含量。杞麓湖土壤样本的质地分类结果如图3所示。

图3 土壤质地三角形

由图3可知，杞麓湖流域农田土壤样本质地类型较为齐全，共有10类，其中以粉质壤土、粉沙质黏壤土、粉沙质黏土为主，少量为沙土、壤质沙土、粉沙土、壤土、沙质壤土、黏壤土、沙质黏壤土。土壤中黏粒、粉粒、沙粒含量波动范围均较大，分别在0～57.51%、12.80～91.03%、0～87.20%波动（表1）。3种颗粒含量中粉粒均值最高，达68.09%，黏粒次之，沙粒最小。这一粒径分布特征很可能与当地高强度的土地利用有关。

从空间分布（图4）来看，黏粒含量大致呈现东高西低、南高北低的空间格局。其中，黏粒含量的高值集中在流域的东南部，这与杞麓湖流域南部土壤黏性较大的实际调研情况一致。粉粒含量在流域东西方向上分布较为平缓，在南北方向上呈现中部高、两侧低的趋势。其中，粉粒的高值集中在杞麓湖南侧。沙粒含量总体偏低，大致呈现西高东低、北高南低的空间格局。其中，沙粒的最高值出现在杞麓湖南侧，根据实地调研发现由于当地农民为了改良土壤性质，在比较黏重的土壤中掺沙，以此保证土壤通透性。

表1 土壤质地的特征统计值

图4 土壤质地的空间趋势

2.2 模型比较

分别以土壤的黏粒、粉粒、沙粒为因变量，波段为400～2 350 nm的反射光谱为自变量，根据式（1）至式（7），运用PLSR、PLS-GWR、PLS-MGWR方法分别建立3种颗粒含量与反射光谱之间的关系，得到9种拟合模型。其中，PLSR模型中的最佳潜变量数（LVs）使用舍一交叉验证法确定。经拟合，黏粒、粉粒、沙粒3种颗粒对应的LVs分别为5、2和3，利用评价指标对9种模型拟合效果进行评价，如表2所示。与PLSR拟合的土壤质地模型相比，MGWR和GWR各模型精度均有所提升，这是由于这2种方法均考虑了光谱在不同空间位置对土壤质地响应的差异。从整体上来看，MGWR对土壤质地（黏粒、粉粒、沙粒）的拟合效果（R2）和稳定性指标均优于PLSR和GWR，这是因为MGWR不仅考虑了土壤光谱与质地关系的空间异质性，还针对各光谱潜变量选择了不同空间尺度，对每个自变量的回归参数都进行局部估计，从而提升了模型精度。此外，PLSR方法建立的粉粒和沙粒模型残差的莫兰指数统计显著，表明残差具有空间依赖性，即存在空间结构，违背了经典统计方法中的残差独立性假设，使得模型精度较低；而GWR和MGWR拟合的各模型残差的莫兰指数在0.005水平上均不显著，说明GWR和MGWR模型残差为空间随机，其拟合结果更为可靠。综上所述，MGWR是一种适于拟合光谱与土壤质地关系的潜力模型。

表2 多尺度地理加权回归、地理加权回归和偏最小二乘回归模型建立效果比较

2.3 MGWR标准化回归系数空间格局分析

由于杞麓湖流域田块较为分散，光谱的不同潜变量与土壤质地的关系在实际情况中可能随着地理位置的变化而变化，呈现出空间非平稳性。多模型比较结果表明，MGWR方法针对实际情况的拟合效果较好，因此，选择MGWR方法对黏粒、粉粒、沙粒分别进行建模，进一步探究每种模型中光谱各潜变量与3种土壤颗粒含量之间的关系。MGWR各系数统计描述见表3。多尺度地理加权回归系数的空间格局如图5所示。

表3 多尺度地理加权回归系数统计描述

2.3.1 黏粒模型 潜变量LV1至LV5整体均具备正向响应。

1）LV1系数反映了光谱潜变量LV1对黏粒含量变化的响应程度。LV1回归系数的取值在0.361～0.387，均值为0.373，标准差为0.004。高值主要集中在杞麓湖的北部区域，这说明在杞麓湖北部，随着黏粒含量增加，潜变量LV1增大；低值出现在杞麓湖流域的最西部区域，在此处，黏粒含量变化较小，潜变量LV1对黏粒的响应较弱。从系数的绝对值来说，LV1在5个变量中对黏粒含量的响应程度最大。

2）LV2对黏粒含量产生的是正向的响应，在空间上呈现出一定的分层结构。高值集中在杞麓湖流域西北部，低值集中在东南部，从流域西北方向至东南方向，潜变量LV2对黏粒含量响应程度逐渐降低。回归系数的取值在0.253～0.304，均值为0.265，标准差为0.013。从系数的绝对值上来看，其响应程度为中等程度。

3）对于潜变量LV3，同样对黏粒含量具有正向响应，黏粒含量越大，潜变量LV3显示出越大的趋势。在空间分布上，高值和低值区域相衔接，分别集中在流域的北部和西北部区域。LV3回归系数的取值在0.091～0.432，均值为0.213，标准差为0.070。从系数的绝对值上来看，潜变量LV3的响应强度最小。

4）潜变量LV4因黏粒含量具有正向响应，其系数取值在0.117～0.440，均值为0.277，标准差为0.108，整体呈现由南北两端向双侧逐渐降低的趋势，在流域西南部，LV4对黏粒的正相关响应程度达到最小。从系数的绝对值上来说，潜变量LV4的响应强度适中。

5）LV5回归系数的取值在-0.097～0.611，均值为0.283，标准差为0.195，呈现出南高北低的空间格局。总体上，潜变量LV5对黏粒的响应为正向，正向效应最大值集中在杞麓湖南部区域，但北部有较小区域产生了负效应，这意味着南部区域潜变量LV5随黏粒含量同向增减，北部负效应区域随黏粒变化呈现相反趋势。在5个潜变量中，潜变量LV5的响应强度略高于LV2。

2.3.2 粉粒模型

1）潜变量LV1系数数值变化较大，在-0.196～1.238波动，均值为0.426，标准差为0.335。整体上呈现出西低东高的格局，高值集中东南部，低值集中在西北部。东南部区域LV1的响应程度最强，此处粉粒的变化程度也最大。西部区域LV1表现出负效应，即LV1由于粉粒含量的提升而有所降低，反之亦然。粉粒模型中，LV1的绝对值远大于LV2，对粉粒的响应程度极高。

2）潜变量LV2正向响应粉粒含量，即粉粒含量与LV2同步增减。LV2系数取值范围在0.040～0.244，均值为0.180，标准差为0.073。在空间分布上，LV2系数在流域北、南、东3个方向数值偏高，西部地区普遍偏低。在最西部区域LV2对粉粒的响应程度最弱，南部区域响应程度最强。从系数的绝对值上来看，LV2对粉粒的响应程度偏弱。

2.3.3 沙粒模型 潜变量LV1、LV2、LV3对沙粒响应程度各不相同，总体上均为正效应。

1）LV1正向响应沙粒含量，沙粒含量越高，LV1越高。LV1回归系数的取值在0.051～1.698，均值为0.422，标准差为0.323。空间上，系数高值低值交错分布，南部区域沙粒的变化幅度达到最大，LV1对沙粒含量响应程度达到最大。从系数的绝对值上来讲，LV1对沙粒响应程度为3个响应潜变量中最大。

图5 多尺度地理加权回归系数的空间格局

2）LV2回归系数的取值在-0.245～0.713，均值为0.194，标准差为0.233。整体上高低值分布较为分散，局部为负，大半为正。负值主要集中在杞麓湖东南区域，沙粒含量随着LV2的增加而降低。从系数的绝对值上来讲，LV2的响应程度较低。

3）LV3系数空间分布上高低值区域参半，低值零散，高值集中在流域北部。LV3回归系数的取值在-0.110～1.408，均值为0.199，标准差为0.243。负值集中在远离湖周的西部，此处沙粒含量对LV3有负向影响，LV3随沙粒的降低而增加。LV3对沙粒的响应程度与LV2相似，都为较弱程度。

3 小结与讨论

本研究使用210个杞麓湖流域的实测土壤质地数据，揭示土壤质地（黏粒、粉粒、沙粒）的空间格局，并对比了偏最小二乘回归、地理加权回归和多尺度地理加权回归3种不同方法拟合土壤质地和光谱关系的效果，最后探究差异化尺度作用下两者关系的空间异质性。

土壤质地描述性统计结果显示，研究区土壤质地类型以粉质壤土、粉沙质黏壤土、粉沙质黏土3种为主要类型，3种颗粒百分含量均值由大到小排序分别是粉粒（68.09%）、黏粒（22.14%）、沙粒（9.77%）。空间趋势上，黏粒含量呈现东高西低、南高北低的格局；粉粒含量在南北方向上呈现中部高、两侧低的特点；沙粒含量呈现西高东低、北高南低的趋势。模型评估结果显示，MGWR的MAE、RMSE和R2指标均优于PLSR和GWR，且残差为空间随机，因此使用MG⁃WR拟合土壤质地与光谱间关系效果最好。MGWR模型回归结果表明，光谱各潜变量随着空间位置变化对土壤质地3种颗粒组分的响应效应显示出不同的空间格局。从系数绝对值上来讲，潜变量LV1对黏粒含量响应程度最大，系数为0.373±0.004；对粉粒含量响应程度最大的潜变量是LV1，系数为0.426±0.335；LV1对沙粒含量的响应程度最大，系数为0.422±0.323。

与GWR相比，MGWR建立了更加有效的空间模型，这主要是由于MGWR考虑了光谱空间上的不同潜变量的不同响应尺度，针对每个变量选择特定带宽作为其拟合过程中的空间尺度指标，从而避免更大误差，使得模型结果更为稳健。本研究考虑了土壤质地与高光谱关系的空间异质性，使用MGWR更好地拟合了不同空间尺度下土壤质地和光谱各潜变量之间的关系，使得空间分析结果更加真实可靠，为开展土壤质地空间异质性研究提供方法借鉴。

由于条件限制，获取的近地高光谱实测数据样本数较少且为点状，并非区域连续的面状数据，因此需要通过大量样本来验证模型的普适性。未来结合机载高光谱、高分五号等遥感影像可以更好地实现大区域尺度的土壤质地空间分析，这将更好地服务于土壤质量风险性评价和生态环境保护。