基于MGWRK的土壤有机质制图及驱动因素研究

奥鹏网院作业

基于MGWRK的土壤有机质制图及驱动因素研究
乔磊1，张吴平2，黄明镜3，王国芳1，任健1

（1山西农业大学资源环境学院，山西太谷 030800；2山西农业大学软件学院，山西太谷 030800；3山西省农业科学院旱地农业研究中心，太原 030031）

摘要：【目的】空间预测是一种获得有机质空间局部细节的重要方法，其准确性对于农田合理管理有着重要意义。本研究通过对比不同的土壤有机质空间制图方法以获得更优的预测精度，在预测的同时揭示环境协变量的空间非平稳性特征及不同环境协变量关系的空间尺度。【方法】选取晋中盆地的7个乡镇作为研究区，对比普通克里金（OK）、回归克里金（RK）、地理加权回归克里金（GWRK）和多重尺度地理加权回归克里金（MGWRK）4种不同方法对土壤有机质的预测能力和效果，并探究影响因子在空间分布中对有机质的影响效应变化和这种影响效应的空间尺度。MGWRK是一种多重尺度地理加权回归（MGWR）与普通克里金方法相结合的方法。【结果】选取坡向、坡度、年均降水量、年平均温度、海拔、植被年总初级生产力、年蒸散量、地形湿度指数、平面曲率、汇流动力指数、地形指数、地形粗糙指数、年平均NDVI为环境协变量参与建模，在多元线性回归建模中，模型统计学意义显著，这表明模型具备统计学意义。从Radius指数来看，各模型模拟效果由好到差依次为RK、OK、GWRK、MGWRK；从制图效果来看，MGWRK与GWRK制图效果相当，从有机质的空间预测图可以看出，土壤有机质在研究区呈现东西两侧偏低、中部偏高的空间格局，其中汾河以东、昌源河流经区域土壤有机质普遍偏高。坡向、年均降水量、年平均温度、海拔、地形指数、年平均NDVI对研究区东部有机质的影响强于西部，而坡度、植被年总初级生产力、年蒸散量、平面曲率、汇流动力指数、地形粗糙指数则表现出截然相反的空间非平稳性特征，地形湿度指数对有机质的影响则体现为北部强南部弱。【结论】MGWRK方法的空间预测精度分别达到了RK方法的69%、OK方法的71.74%、GWRK方法的71.17%。MGWRK在空间非平稳性特征的解释能力和空间可视化表现良好，为有机质的预测和空间非平稳性特征的描述提供方法借鉴。

关键词：土壤有机质；土壤制图；多重尺度地理加权回归克里金；地理加权回归克里金；回归克里金；空间非平稳性特征

0 引言
【研究意义】土壤有机质是土壤的重要组成部分[1]，土壤有机质对于作物生长[2]、土壤肥力[3]、环境保护[4]等方面有着重要的作用，出于有机质的重要性，土壤有机质制图就尤为重要，土壤有机质制图对于揭示土壤有机质空间分布、农业管理以及生态环境等诸多方面都有着巨大的推动作用。对土壤有机质的空间分布及预测能够为农业管理和生态建设提供有力依据，探究土壤有机质与环境协变量之间的空间非平稳特征可以更加明确不同影响因子对土壤有机质的影响程度及空间变化。【前人研究进展】近年来利用地理加权回归方法分析土壤有机质的空间分布及空间非平稳性特征已有报道[5-11]，吴春生等[5]利用高程、年均气温、年均降雨量以及归一化植被指数4个协变量完成地理加权回归建模，提高了有机质空间预测精度；陈琳等[6]对比了不同空间回归技术的预测精度，证实了地理加权回归克里金（GWRK）的预测精度优于普通克里金（OK）、逐步回归克里金（SWRK）；Zeng[7]等对比了混合地理加权回归方法、GWR等多种空间预测方法的模拟效果，展示了混合地理加权回归的良好应用前景。【本研究切入点】地理加权回归的应用非常广泛，涉及林业[12-13]、环境[14-16]、城镇化[17]、交通[18]等众多领域，但是利用多重尺度地理加权回归克里金方法分析有机质分布及空间非平稳性特征却鲜有报道。【拟解决的关键问题】本研究利用多重尺度地理加权回归克里金（multiscale geographically weighted regression kriging，MGWRK）方法选取坡向、坡度、年均降水量、年平均温度、海拔、植被年总初级生产力、年蒸散量、地形湿度指数、平面曲率、汇流动力指数、地形指数、地形粗糙指数、年平均NDVI共13个环境协变量对土壤有机质进行建模，探究不同协变量的空间非平稳性特征和尺度效应，并与普通克里金（ordinary kriging，OK）、回归克里金（regression kriging，RK）、地理加权回归克里金（geographically weighted regression kriging，GWRK）对比分析各自的预测精度与制图效果。旨在获得更高精度的有机质空间分布结果以及不同影响因素的空间非平稳性特征，为土壤管理和空间分析提供依据和参考。

1 材料与方法
1.1 研究区概况
本研究选取晋中平原中平遥县的7个乡镇作为研究区域。研究区位于平遥县西北部，是山西省晋中盆地腹地位置，境内多为平川。研究区属温带大陆性半干旱季风气候，平均气温12℃，年平均相对湿度58%，年均降水量502 mm，主要分布在7—9月份。主导风向春秋冬季为西北风，夏季多为东南风，年平均风速2.1 m·s-1。种植作物以玉米为主，土壤类型以潮土为主，部分地区为盐土、褐土。

1.2 土壤采样与数据来源
样品采集于2010年春季作物种植前，依据《耕地地力调查与质量评价技术规程》（NY/T 1634—2008），用不锈钢土钻等工具，在每个田块中采用“S”法布点采样，以GPS定位点作为中心，向四周辐射多个分样点，每个混合土样取15个以上分样点，每个分样点的采土部位、深度、数量保持一致，采样深度为0—20 cm，共采集样点2 616个，采样点分布见图1。土壤有机质（SOM）采用重铬酸钾加热法[19]测定。数字高程数据（digital elevation model，DEM）来自中国科学院计算机网络信息中心地理空间数据云网站（www.gscloud.cn），各地形因子均由DEM数据计算得出。气象数据来自国家气象信息中心中国气象数据网（http://data.cma. cn/）。研究中所用遥感数据均来自NASA网站（https://earthdata.nasa.gov/），2009年植被年总初级生产力（gross primary productivity，GPP）数据来自MOD17数据集，2009年蒸散（evapotranspiration，ET）数据来自MOD16数据集，2009年平均NDVI数据来自MOD13数据集。环境协变量数据列于表1。

表1 环境协变量数据列表

Table 1 List of environmental covariate data


width=431.45,height=308.4
图1 研究区概况及样点分布

Fig. 1 Study area and soil sample sites

1.3 研究方法
1.3.1 普通克里金（OK） 普通克里金（ordinary kriging，OK）基于空间自相关性和二阶平稳假设[20]，依据估算误差最小和半方差函数分析，赋予一定邻域内样点不同权重后，计算样点值加权和即为插值结果。其中半变异函数公式为：

width=154.2,height=27.95 （1）
式中，γ(h)为半变异函数值；N(h)为距离等于h时的点对总数；Z(xi)为区域变量在位置xi处的实测数值；Z(xi+h)为在xi+h处的实测数值。

1.3.2 回归克里金（RK） 回归克里金[21]是多元线性回归（multiple linear regression，MLR）与克里金插值方法的结合，通过多元线性回归方法获取预测模型，对预测产生的残差进行普通克里金估计，将趋势值与残差项相加得到最终预测结果。本研究使用的多元线性回归方法为普通最小二乘法（ordinary least square，OLS），是因变量yi与自变量xi的多元线性函数：

width=87.05,height=28.5 （2）
式中，β0为常数项即斜率，βi为每个自变量的回归系数，εi为随机误差项。

回归克里金[22]的公式如下：

Z=m(xi)+ e(xi) （3）

式中，Z为RK模型预测值，m(xi)为MLR过程的趋势值，e(xi)是利用OK法估计的MLR过程残差，后续的GWRK与MGWRK方法回归项计算均与此类似。

1.3.3 地理加权回归克里金（GWRK） GWRK[23]是地理加权回归（geographically weighted regression，GWR）模型与克里金插值方法的结合，也可视为是RK的改进，RK是全局回归模型[24]，在RK中多元线性回归假设因变量与自变量的关系在空间上不随着空间位置的变化而变化，而GWRK为局部回归模型。地理加权回归[25]是众多空间回归技术应用最为广泛的一种，建立模型的过程中考虑到了各样点的空间位置以及不同位置中自变量与因变量关系的变化，揭示在空间尺度中自变量与因变量的关系。

假定共n个观测点，观测点为i∈{1, 2, ..., n}，每个观测点的位置为（ui,vi），参与建模的自变量共m个。GWR模型[26]的表达式如下：

width=108.55,height=29.55 （4）
式中，xij是i观测点的第j个自变量，βi(ui,vi)为在(ui,vi)位置上第j个自变量的回归系数，εi为随机误差项，yi为因变量。

核函数和有效带宽[26]是GWR建模过程中最重要的两项指标。核函数可用于确定空间权重，包括高斯函数、双平方函数和指数函数，核函数在运行过程中可采用不同的类型，包括固定核函数和自适应核函数[13]，若观测点均匀分布时可采用固定核函数，否则选择自适应核函数，本研究使用自适应高斯函数进行空间权重确定。由于GWR是移动窗口回归方法的一种，而带宽则是核函数的重要性质之一，带宽与核函数原理详见文献[26]，有效带宽的确定方法包括修正赤池信息量准则（corrected akaike information criterion，AICc）、赤池信息量准则（akaike information criterion，AIC）、贝叶斯信息准则（bayesian information criterion，BIC）等，本研究采用修正赤池信息量准则（AICc）确定模型带宽，原则是AICc最小化，以此作为确定带宽的标准，其中AICc的计算公式如下：

width=125.2,height=40.3 （5）
式中，σ为样本标准差，n为样本含量，tr(S)为矩阵S的迹，矩阵S是模型的帽子矩阵。

在GWR模型中，确定后的最优有效带宽是所有自变量共用的，即在建模中所有自变量的有效带宽都是一样的，空间权重矩阵也是一致的。在GWR模型中未考虑不同环境协变量与因变量间关系不同的的空间尺度，而是预设所有环境协变量的空间非平稳性特征的空间尺度一致。

1.3.4 多重尺度地理加权回归克里金（MGWRK） MGWRK是多重尺度地理加权回归（multiscale geographically weighted regression，MGWR）方法与克里金方法的结合，也可看做是对于GWRK的改进。MGWR[27]是Fotheringham等提出的一种基于有效带宽优化从而使每个自变量都具有不同的空间带宽的地理加权回归方法。与GWR不同的是，GWR确定的模型带宽是所有自变量共用的，而实际情况中不同的自变量的空间尺度效应和对因变量的影响程度是不一样的，所以MGWR方法考虑到了这一点，MGWR方法使每个自变量都有各自的有效带宽，即为多重尺度。MGWR模型如下：

width=112.85,height=29.55 （6）
式中，βbwj是经过第j个自变量有效带宽校正后的回归系数。

在本研究中，MGWR与GWR模型的核函数和核类型均保持一致，核函数均采用高斯函数，核函数采用自适应核函数。两种模型最大的区别就是确定有效带宽的方法，有效带宽的本质即是空间尺度，也正是确定带宽方法的不同决定了MGWR多重尺度的特性。有效带宽的选择对于模型计算结果影响巨大，我们可以将有效带宽看作一个具有平滑性质的参数，当有效带宽太大有可能造成整个研究区域的参数相似，掩盖研究变量的局部特征；而有效带宽太小则会造成过多的局部变异，难以识别研究空间变量之间的规律[28]。选取一个最优带宽避免以上两种极端情况，提供良好的空间关系描述效果就显得尤为重要。

在MGWR模型中，在不同的自变量中有着不同的带宽，这也意味着在建模中不同自变量在同一位置将具有完全不同的空间权重矩阵，在GWR模型中使用的具有平滑作用的参数（有效带宽）将不再适用于MGWR模型。MGWR模型中采用向后拟合（Back-fitting）算法[27]进行模型校正。首先使用GWR模型计算的各位置自变量参数以及模拟结果（y）作为初始化，使用这些初始值（βbwjxj）以及求得的残差（ε=y-Sm j=0βbwjxij），残差加当前的βbwjxj（βbwjxj+ε），利用GWR模型将该值与x0（第一个自变量）进行回归，此时得到因变量与x0之间的参数βbw0的局部估计，即得到了因变量与x0之间的最优带宽bw0。接下来使用更新后的βbwjxj与下一个自变量x1进行同样的过程，得到直到所有的自变量（x0,x1,…,…, xm）都完成该过程（图2）。

width=319.15,height=212.25
图2 MGWR模型的向后拟合（Back-fitting）算法[27]

Fig. 2 Back-fitting algorithm for multiscale geographically weighted regression model (MGWR)

1.4 模型评价
本研究使用平均误差（ME）、平均绝对误差（MAE）、均方根误差（RMSE）评价模型精度[29]，这三项指数的公式如下：

width=112.3,height=28.5 （7）
width=118.2,height=27.95 （8）

width=128.4,height=30.65 （9）

式中，n是用于验证的观测点数量，Z(xi)是实测值，Z*(xi)是模型预测值。理论上MAE和RMSE越接近于0模型预测效果越好。在一些研究中3个指标很难同时达到最佳水平，所以除上述指标外，引进Yang等[29]提出的Radius指数（式10），以便于结合3个指数确定最优模型。

width=171.95,height=19.35 （10）
式中，Radius即为Radius指数，MES、MAES、RMSES是使用Z得分法进行标准化后的ME、MAE、RMSE，该指数越小则说明模型效果越好。

2 结果
2.1 环境协变量选择及MLR建模
土壤有机质的影响因素众多，在多元线性回归建模选取变量的过程当中需要考虑如下因素[30-31]：（1）各环境协变量对于有机质的分布有显著影响，影响极其微弱的不予考虑；（2）各环境协变量在空间分布中具有一定的空间异质性；（3）各环境协变量之间不存在多重共线性，这也是多元线性回归建模的必要条件之一。

本研究选取坡向、坡度、年均降水量、年平均温度、海拔、植被年总初级生产力、年蒸散量、地形湿度指数、平面曲率、汇流动力指数、地形指数、地形粗糙指数、年平均NDVI共13项指标（图3）作为建模的环境协变量，以有机质作为因变量。其中环境协变量之间的多重共线性诊断以容差和方差膨胀因子（VIF）作为依据，容差越小，表明该自变量作为因变量进行回归分析时被其他变量解释的程度越高，越可能存在严重的共线性，容差合理的范围是（0.1，+∞）；方差膨胀因子是容差的倒数，若其值≥10，说明自变量间可能存在严重的共线性现象[32]，模型的共线性诊断结果见表2。其中VIF均未超过10，故认为在该模型中不存在严重的多重共线性现象。

对土壤有机质以及各环境协变量进行描述性统计（表3），可以看出，不同环境协变量的变异系数是不同的，其中地形指数的变异系数是最大的，年均降水量的变异系数是最小的。

表2 模型多重共线性检验

Table 2 Model multicollinearity test


多元线性回归（MLR）的最终拟合模型为：

SOM=-131.09-0.0016×Aspects-0.01×Slope +0.46×Pre-3.74×Tem-0.04×Height-0.009×GPP+ 0.01×ET-0.08×TWI-92.15×PC+0.000007×SPI+0.06×TPI+0.03×TRI-1.51×NDVI

式中，SOM为土壤有机质，Aspects为坡向，Slope为坡度，Pre为年均降水量，Tem为年平均温度，Height为海拔，GPP为植被年总初级生产力，ET为年蒸散量，TWI为地形湿度指数，PC为平面曲率，SPI为汇流动力指数，TPI为地形指数，TRI为地形粗糙指数，NDVI为年平均归一化植被指数。对模型进行方差检验，F统计量[33]为43.39，模型具有统计学显著意义，即回归方程有效。

同样，在GWRK与MGWRK模型中也采用同样的环境协变量作为模型自变量。

2.2 各模型预测结果空间变异
对不同模型的残差结果空间半变异函数拟合（图4），模型参数及最优模型见表4。3种模型残差值半变异函数的最优拟合模型从指数模型、高斯模型和球状模型[34]中选取。GWR残差的半变异函数最优拟合模型为高斯模型，其他最优拟合模型均为指数模型。从块金系数[35]可以看出，各模型残差最优模型的块金系数均小于50%，这意味着模型预测过程的残差具有空间结构性，使用各自的半变异函数进行空间插值，结果见图5。

2.3 各模型评价
各模型评价指标见表5，从Radius指数来看，各模型模拟效果由好到差依次为：RK、OK、GWRK、MGWRK。MGWRK模型的空间预测精度相对于其他模型欠佳，而RK模型的空间预测精度则表现最优。

表3 土壤有机质及环境协变量描述性统计

Table 3 Descriptive statistics of soil organic matter and environmental covariates


表4 各模型结果空间变异

Table 4 Semi-variogram of the results of each model


研究从指数模型、高斯模型、球状模型中选取最佳模型

Study select optimization model from Exponential model(E), Gaussian Model(M) and spherical model(S)

表5 各模型评价指标

Table 5 Evaluation index of each model


width=425.55,height=417.5
坡向以正北为0°，顺时针计数，-1为平地 Aspect direction is defined as zero degrees due north, clockwise count, -1 is flat

图3 环境协变量空间分布

Fig. 3 Spatial distribution of environmental covariates

width=415.35,height=119.3
图4 半变异函数图

Fig. 4 Semivariogram

width=420.7,height=139.7
图5 各模型残差图

Fig. 5 Map of different Model residuals (MLR residuals, GWR residuals, MGWR residuals)

2.4 各模型土壤制图结果
根据以上分析及模型参数，依照不同模型绘制土壤有机质分布图（图6）。可以看出4种不同的土壤制图方法，得到空间基本趋势是类似的。MGWRK与GWRK方法的效果相对其他两种方法更为平滑，这可能是由于两种方法考虑了各环境协变量的空间非平稳性特征导致的，OK和RK两种方法高低值交错明显，“牛眼”现象明显，与有机质空间实际分布情况不符。从制图效果来看，GWRK方法与MGWRK方法的效果相当。

土壤有机质在研究区呈现东西两侧偏低，中部偏高的空间格局。其中汾河以东、昌源河流经区域土壤有机质普遍偏高，沿河土壤的土壤水分受到河流影响，而土壤水分影响着土壤有机质的分解与转化过程，土壤微生物的活动与土壤水分密切相关，土壤水分过高时土壤中形成厌氧环境，从而改变土壤有机质的分解过程与产物，进而影响土壤有机质含量；土壤水分过低时则有可能影响土壤微生物的生存状况，降低土壤微生物的活性[36]。

width=395.45,height=292.3
该图制图的空间分辨率为30 m The mapping spatial resolution is 30 m

图6 各模型土壤制图结果

Fig. 6 Soil mapping results of each model

2.5 有机质影响因子回归系数空间分布
结合土壤有机质及回归参数的统计特征、多元线性回归方程可知，在全局尺度上土壤有机质与全氮、年均降水量呈现正相关关系，土壤有机质与其余指标均呈现负相关关系。但在实际情况中随着地理位置的变化因变量与自变量之间的关系也是变化的，即空间非平稳性，MGWR模型针对不同环境协变量使用不同的空间尺度，在进行空间预测的同时也对自变量在每个空间位置的回归参数局部估计[37]。为了避免建模过程当中不同变量的数量级对回归系数的相对大小造成影响，在对比影响因素随着位置变化影响程度的不同时使用经过z分数方法标准化后的局部回归系数进行空间非平稳性特征描述，即标准化回归系数。

图7-a中，坡向的标准化回归系数范围在-2.45—2.26之间，整体呈现东高西低的空间格局，其中正值区域占绝大部分。在研究区东部区域为正相关影响，即坡向角度越大有机质增长幅度越大，在南政乡回归系数达到了最小值，在此处呈现负相关影响，即坡向角度越大有机质减少幅度越大。

图7-b中，坡度的标准化回归系数范围在-1.82—1.39之间，整体呈现东低西高的空间格局，正负值区域参半，正负值预取之间分界明显，大致表现为汾河以东为负相关影响，汾河以西为正相关影响。

图7-c中，年均降水量的标准化回归系数范围在-1.41—2.37之间，整体呈现东高西低的空间格局，在东部区域降水量偏少（低于500 mm），此处表现为正相关影响，即在该部分区域内降水量越大有机质随之增加；而在西部区域降水量相对较大，表现为负相关影响，即在此区域内降水量越大有机质越小。

图7-d中，年平均温度的标准化回归系数范围在-1.38—1.30之间，整体呈现东高西低的空间格局，在年平均温度小于12℃的区域表现为正相关影响，在南政乡这种影响效应达到了最强，即随着年平均温度的增加有机质随之增加且增长幅度最大；年平均温度大于14℃的区域负相关影响达到最强，该区域内随着温度的增加有机质减少幅度达到最大。

图7-e中，海拔的标准化回归系数范围在-1.64—2.06之间，整体呈现东高西低的空间格局，正负值区域参半。在西南部区域，海拔是整个研究区最低的，海拔对于有机质的负相关影响效应达到最大，即海拔降低有机质随之减少的幅度达到最大。

图7-f中，植被年总初级生产力的标准化回归系数范围在-2.07—1.26之间，整体呈现东低西高的空间格局，正值区域占大部分。在东部区域GPP水平达到区域最高，而此处回归系数则达到最低水平，即随着GPP的增加有机质减少的幅度达到最大；而南部GPP水平为区域最低，回归系数却达到最高，即此处GPP越小有机质增长幅度越大。

图7-g中，年蒸散量的标准化回归系数范围在-2.05—0.93之间，整体呈现东低西高的空间格局，正值区域占大部分。研究区西部蒸散量达到最大，与此同时年蒸散量对有机质的正相关影响效应达到最大，东部蒸散量最小，相关系数也达到最小，负相关影响程度达到最强。

图7-h中，地形湿度指数的标准化回归系数范围在-3.09—2.19之间，整体呈现北高南低的空间格局，正负值区域参半。TIW在宁固乡北部影响程度最为微弱，此时TWI的范围大致为10.81—10.91，TWI反映土壤水分的干湿状况[38]，从而间接地说明在一定土壤水分范围内对有机质的影响程度达到最弱。

图7-i中，平面曲率的标准化回归系数范围在-1.42—1.31之间，整体呈现东低西高的空间格局，正负值区域参半。结合图3，研究区的平面曲率在宁固乡最接近0，即地形更为平整[6]，此处有机质也是增长幅度最大的，表现为正相关影响。东部地形变化较为剧烈的区域均表现为负相关影响。

图7-j中，汇流动力指数的标准化回归系数范围在-2.07—1.21之间，整体呈现东低西高的空间格局，正负值区域参半。与SPI的空间分布格局恰恰相反，SPI越大对有机质的负相关影响越强，反之SPI越小则正相关影响越强。

图7-k中，地形指数的标准化回归系数范围在-1.65—1.36之间，整体呈现东高西低的空间格局，正值区域占大部分。TPI在研究区东部达到最大，TPI的标准化回归系数也达到了最大水平，即在地形变化剧烈[38]的区域有机质变化也是剧烈的。

图7-l中，地形粗糙指数的标准化回归系数范围在-1.48—1.71之间，整体呈现东低西高的空间格局，正负值区域参半。TRI在研究区呈现东高西低的空间格局，恰恰与标准化回归系数相反，TRI＜1.98时标准化回归系数达到最小，即地形异质性[39]越小有机质的增长幅度越大。当TRI＞8.59时标准化回归系数达到最大，即地形异质性越大有机质的减少幅度越大。

图7-m中，年平均NDVI的标准化回归系数范围在-2.83—1.42之间，整体呈现北高南低的空间格局，正负值区域参半。结合图3，NDVI在研究区南部达到最小，北部达到最大，在北部有机质随着年平均NDVI的增加而增加，越往北增加幅度越大，反之越往南随着NDVI的减少有机质随之降低且降低幅度变大。

2.6 MGWR与GWR模型带宽对比
GWR模型及MGWR各协变量的最优带宽见图8，可以清楚地看出，两种模型在带宽方面有着截然不同的特性，在MGWR模型中各协变量均有着不同的带宽，而GWR模型仅有一种带宽，不同环境协变量的空间非平稳性特征有着不同的空间尺度，GWR模型中忽视了这种空间尺度效应，而MGWR模型对这种特性加以描述并参与建模，达到了可观的效果。

width=440.05,height=416.95
图7 土壤有机质影响因子标准化回归系数分布

Fig. 7 Standardized regression coefficient distribution of soil organic matter influencing factors

width=212.8,height=165.5
a：坡向；b：坡度；c：年均降水量；d：年平均温度；e：海拔；f：植被年总初级生产力；g：年蒸散量；h：TPI；i：平面曲率；j：SPI；k：TPI；l：TRI；m：NDVI；n：截距项；o：GWR模型

a: Aspect of MGWR; b: Slope of MGWR; c: Annual mean precipitation of MGWR; d: Annual mean temperature of MGWR; e: Height of MGWR; f: Gross primary productivity of MGWR; g: Evapotranspiration of MGWR; h: TPI of MGWR; i: Plan curvature of MGWR; j: SPI of MGWR; k: TPI of MGWR; l: TRI of MGWR; m: NDVI; n: Intercept of MGWR; o: Geographically Weighted Regression (GWR)

图8 MGMR与GWR模型带宽对比

Fig. 8 Bandwidth comparison of MGMR and GWR

3 讨论
为了更加全面地探究MGWRK方法精度比GWRK方法精度低的原因，对MLR、GWR和MGWR 3种方法的预测精度也进行评价（表6），并结合文章2.3部分的模型评价绘制Radius指标的弧度图（图9）。由图表可以看出，MGWR模型的Radius指标比GWR模型降低了141.8024，GWRK的Radius指标相对于GWR降低了218.94，MGWRK相对于MGWR仅降低了38。MGWR模型的表现优于GWR模型，但MGWRK模型的表现却逊于GWRK模型。结合图5，对比不同模型的残差空间分布图，GWR的残差高低交错结构明显，这也是GWR模型中带宽固定的结果；而MGWR的残差空间分布图高低交错并没有GWR的残差明显，这是因为MGWR模型中不同自变量的带宽不同对残差的空间分布造成影响而导致MGWRK预测精度相对于GWRK模型有所降低。

GWR模型模拟结果的信息准则指数分别为：AICc（6705.059）、AIC（6665.335），MGWR模型计算结果信息准则指数分别为AICc（6123.467）、AIC（6090.237），MLR模型计算结果信息准则指数分别为AICc（6940.712）、AIC（6938.527）。信息准则指数越小模型效果越好，从信息准则指数看，MGWR、GWR模型的模拟效果均优于MLR模型。MLR、GWR、MGWR的拟合度分别为0.178、0.366、0.484，各自的调整后R2分别为0.174、0.309、0.442，从拟合度来看，与信息准则和Radius指标表现的结果一致。

由于MGWR模型多重尺度的特性，使得各种环境协变量对于因变量影响程度的空间尺度效应得以很好的描述，这种多重尺度特性提高了MGWR模型的预测能力，但是由于多重尺度对于MGWR残差的影响，这种特性反而降低了MGWRK的预测精度。MGWR模型针对不同环境协变量的空间影响能力给予了不同的空间带宽，切实地解释了不同环境协变量不同的尺度特性，而GWR模型由于统一带宽使得环境协变量的空间多重尺度效应被掩盖。

width=196.65,height=181.6
图9 不同模型的Radius指数弧形图

Fig. 9 Radius indices shown in an arc diagram for different models

表6 其他模型评价

Table 6 Evaluation index of others model


4 结论
研究区土壤有机质呈现东西两侧偏低、中部偏高的空间格局，其中汾河以东、昌源河流经区域土壤有机质普遍偏高。多重尺度地理加权回归克里金（MGWRK）能够良好的实现土壤有机质的空间预测，并且能够得到各环境协变量随着位置变化对有机质影响效应的变化，这种影响效应在不同环境协变量中有着不一样的空间尺度。MGWRK方法的空间预测精度分别达到了RK方法的69.52%、OK方法的71.74%、GWRK方法的71.17%。

土壤有机质的影响因素众多，包括温度、水分状况、通气情况、植物残体特性、土壤特性等，本文仅选取了部分协变量参与建模，在后续研究中可以引入更加详尽综合的环境变量以提高预测精度，并且多重尺度地理加权回归克里金方法的特性与表现有待进一步研究。

References

[1] 黄昌勇, 徐建明. 土壤学. 北京: 中国农业出版社, 2010.

HUANG C Y, XU J M. Soil Science. Beijing: China Agriculture Press, 2010. (in Chinese)

[2] 银敏华, 李援农, 李昊, Yang Yang, 徐袁博, 张天乐. 垄覆黑膜沟覆秸秆促进夏玉米生长及养分吸收. 农业工程学报, 2015, 31(22): 122-130.

YIN M H, LI Y N, LI H, YANG Y, XU Y B, ZHANG T L. Ridge-furrow planting with black film mulching over ridge and corn straw mulching over furrow enhancing summer maize’s growth and nutrient absorption. Transactions of the Chinese Society of Agricultural Engineering, 2015, 31(22): 122-130. (in Chinese)

[3] 周伟, 王文杰, 何兴元, 张波, 肖路, 王琼, 吕海亮, 魏晨辉. 哈尔滨城市绿地土壤肥力及其空间特征. 林业科学, 2018, 54(9): 9-17.

ZHOU W, WANG W, HE X Y, ZHANG B, XIAO L, WANG Q, LV H L, WEI C H. Soil fertility and spatial variability of urban green land in Harbin. Scientia Silvae Sinicae, 2018, 54(9): 9-17. DOI: 10.11707/ j.1001-7488.20180902. (in Chinese)

[4] 庞蕊, 刘敏, 李美玲, 徐兴良. 土壤碳排放组分区分的研究进展. 生态学杂志, 2017, 36(8): 2327-2335.

PANG R, LIU M, LI M L，XU X L. Advances in research on the separation of soil carbon emissions. Chinese Journal of Ecology, 2017, 36(8): 2327-2335. DOI: 10.13292/j.1000-4890.201708.007. (in Chinese)

[5] 吴春生, 刘高焕, 刘庆生, 黄翀, 张韵婕, 管续栋. 蒙古高原中北部土壤有机质空间分布研究. 资源科学, 2016, 38(5): 994-1002.

WU C S, LIU G H, LIU Q S, HUANG C, ZHANG Y J, GUAN X D. The spatial distribution of soil organic matter on the north-central Mongolian Plateau. Resources Science, 2016, 38(5): 994-1002. DOI: 10.18402/resci.2016.05.18. (in Chinese)

[6] 陈琳, 任春颖, 王宗明, 张柏. 基于克里金插值的耕地表层土壤有机质空间预测. 干旱区研究, 2017, 34(4): 798-805.

CHEN L, REN C Y, WANG Z M, ZHANG B. Prediction of spatial distribution of topsoil organic matter content in cultivated land using Kriging methods. Arid Zone Research, 2017, 34(4): 798-805. DOI: 10.13866/j.azr.2017.04.12. (in Chinese)

[7] ZENG C Y, YANG L, ZHU A X, DAVID G R, LIU J, LIU J Z, QIN C Z, WANG D S. Mapping soil organic matter concentration at different scales using a mixed geographically weighted regression method. Geoderma, 2016, 281: 69-82.

[8] 王合玲, 张辉国, 吕光辉. 艾比湖流域土壤有机质与土壤因子响应关系的空间非平稳性分析. 水土保持研究, 2016, 23(2): 224-228, 357, 358.

WANG H L, ZHANG H G, LV G H. Analysis of spatial non- stationary relationships between soil organic matter and soil factors in Ebinur Lake Basin. Researcher of Soil and Water Conservation, 2016, 23(2): 224-228, 357, 358. DOI: 10.13869/j.cnki. rswc.2016.02.042. (in Chinese)

[9] 马泉来, 高凤杰, 张志民, 单培明, 韩文文, 周军, 曲茉莉. 我国东北黑土丘陵区小流域土壤有机质空间分布模拟. 环境科学研究, 2016, 29(3): 382-390.

MA Q L, GAO F J, ZHANG Z M, SHAN P M, HAN W W, ZHOU J, QU M L. Simulation of spatial distribution of soil organic matter in a mollisol watershed in Northeastern China. Research of Environmental Sciences, 2016, 29(3): 382-390. DOI: 10.13198/j.issn.1001-6929.2016. 03.09. (in Chinese)

[10] YE H C, HUANG W J, HUANG S Y, HUANG Y F, ZHANG S W, DONG Y Y, CHEN P F. Effects of different sampling densities on geographically weighted regression kriging for predicting soil organic carbon. Spatial Statistics, 2017, 20: 79-61.

[11] 石淑芹, 曹祺文, 李正国, 许恒周. 气候与社会经济因素对土壤有机质影响的空间异质性分析——以黑龙江省中部地区为例. 中国生态农业学报, 2014, 22(9): 1102-1112.

SHI S Q, CAO Q W, LI Z G, XU H Z. Influence of spatial heterogeneity of climatic and socio-economic factors on soil organic matter —A case study of the central Heilongjiang Province, China. Chinese Journal of Eco-Agriculture, 2014, 22(9): 1102-1112. DOI: 10.13930/j.cnki.cjea.140097. (in Chinese)

[12] 杜一尘, 李明泽, 范文义, 王斌. 基于地理加权回归模型与林火遥感数据估算森林年龄. 林业科学, 2019, 55(6): 184-194.

DU Y C, LI M Z, FAN W Y, WANG B. Estimation of forest stand age based on gwr model and forest fire remote sensing data. Scientia Silvae Sinicae, 2019, 55(6): 184-194. DOI: 10.11707/j.1001-7488. 20190622. (in Chinese)

[13] 王海宾, 侯瑞萍, 郑冬梅, 高秀会, 夏朝宗, 彭道黎. 基于地理加权回归模型的亚热带地区乔木林生物量估算. 农业机械学报, 2018, 49(6): 184-190.

WANG H B, HOU R P, ZHENG D M, GAO X H, XIA C Z, PENG D L. Biomass estimation of arbor forest in subtropical region based on geographically weighted regression model. Transactions of the Chinese Society of Agricultural Machinery, 2018, 49(6): 184-190. DOI: 10.6041/j.issn.1000-1298.201.06.021. (in Chinese)

[14] 亚森江·喀哈尔, 茹克亚·萨吾提, 尼加提·卡斯木, 尼格拉·塔什甫拉提, 张飞, 阿不都艾尼·阿不里, 师庆东. 优化光谱指数的露天煤矿区土壤重金属含量估算. 光谱学与光谱分析, 2019, 39(8): 2486-2494.

YASENJIANG K, RUKEYA S, NIJAT K, NIGARA T, ZHANG F, ABDUGHENI A, SHI Q D. Estimation of heavy metal contents in soil around open pit coal mine area based on optimized spectral index. Spectroscopy and Spectral Analysis, 2019, 39(8): 2486-2494. DOI: 10.3964/j.issn.1000-0593(2019)08-2486-09. (in Chinese)

[15] 许珊, 邹滨, 宫俊霞. 2001～2015年中国城镇化与PM_(2.5)浓度时空关联特征. 中国环境科学, 2019, 39(2): 469-477.

XU S, ZOU B, GONG J X. Analysis of the spatial-temporal association between urbanization and PM2.5 concentration during 2001~2015 period in Mainland China. China Environmental Science, 2019, 39(2): 469-477. DOI: 10.19674/j.cnki.issn1000-6923.2019.0057. (in Chinese)

[16] XU C C, ZHAO J Y, LIU P. A geographically weighted regression approach to investigate the effects of traffic conditions and road characteristics on air pollutant emissions. Journal of Cleaner Production, 2019, 239: 118084.

[17] 王海军, 张彬, 刘耀林, 刘艳芳, 徐姗, 邓羽, 赵雲泰, 陈宇琛, 洪松. 基于重心-GTWR模型的京津冀城市群城镇扩展格局与驱动力多维解析. 地理学报, 2018, 73(6): 1076-1092.

WANG H J, ZHANG B, LIU Y L, LIU Y F, XU S, DENG Y, ZHAO Y T, CHEN Y C, HONG S. Multi-dimensional analysis of urban expansion patterns and their driving forces based on the center of gravity-GTWR model: A case study of the Beijing-Tianjin-Hebei urban agglomeration. Acta Geographica Sinica, 2018, 73(6): 1076-1092. DOI: 10.11821/dlxb201806007. (in Chinese)

[18] 杨晴青, 刘倩, 尹莎, 张戬, 杨新军, 高岩辉. 秦巴山区乡村交通环境脆弱性及影响因素——以陕西省洛南县为例. 地理学报, 2019, 74(6): 1236-1251.

YANG Q Q, LIU Q, YIN S, ZHANG J, YANG X J, GAO Y H. Vulnerability and influencing factors of rural transportation environment in Qinling-Daba mountainous areas: A case study of Luonan county in Shaanxi province. Acta Geographica Sinica, 2019, 74(6): 1236-1251. DOI: 10.11821/dlxb201906012. (in Chinese)

[19] 鲍士旦. 土壤农化分析. 3版. 北京: 中国农业出版社, 2015.

Bao S D. Soil and Agrochemical Analysis. 3rd ed. Beijing: Chinese Agriculture Press, 2015. (in Chinese)

[20] 刘爱利, 王培法, 丁园圆. 地统计学概论. 北京: 科学出版社, 2012.

LIU A L, WANG P F, DING Y Y. Introduction to Geostatistics. Beijing: Science Press, 2012. (in Chinese)

[21] WU Z H, WANG B Z, HUANG J L, AN Z H, JIANG P, CHEN Y Y, LIU Y F. Estimating soil organic carbon density in plains using landscape metric-based regression Kriging model. Soil & Tillage Research, 2019, 195: 104381.

[22] HENGL T, HEUVELINK G B M, ROSSITER D G. About regression- kriging: From equations to case studies. Computers and Geosciences, 2007, 33(10): 1301-1315.

[23] SHEN Q S, WANG Y, WANG X R, LIU X, ZHANG X Y, ZHANG S L. Comparing interpolation methods to predict soil total phosphorus in the Mollisol area of Northeast China. Catena, 2019, 174: 59-72.

[24] 杨顺华, 张海涛, 郭龙, 任艳. 基于回归和地理加权回归Kriging的土壤有机质空间插值. 应用生态学报, 2015, 26(6): 1649-1656.

YANG S H, ZHANG H T, GUO L, REN Y. Spatial interpolation of soil organic matter using regression Kriging and geographically weighted regression Kriging. Chinese Journal of Applied Ecology, 2015, 26(6): 1649-1656. DOI: 10.13287/j.1001-9332.20150331.023. (in Chinese)

[25] 聂磊, 舒红, 刘艳. 复杂地形地区月平均气温(混合)地理加权回归克里格插值. 武汉大学学报(信息科学版), 2018, 43(10): 1553-1559.

NIE L, SHU H, LIU Y. Interpolation of monthly average temperature by using (mixed) Geographically Weighted Regression Kriging in the complex terrain region. Geomatics and Information Science of Wuhan University, 2018, 43(10): 1553-1559. DOI: 10.13203/j.whugis20160433. (in Chinese)

[26] A. STEWART F, CHRIS B, MARTIN C. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley, 2002.

[27] A. STEWART F, YANG W B, KANG W. Multiscale Geographically Weighted Regression (MGWR). Annals of the American Association of Geographers, 2017, 107(6): 1247-1265.

[28] LI Z Q, TAYLOR O, STEWART F, KANG W, LEVI W, YU H C, LUO W. MGWR 1. 0 User Manual. Spatial Analysis Research Center (SPARC). [2019-11-1] https: //sgsup.asu.edu/sparc/mgwr.

[29] YANG S H, LIU F, SONG X D, LU Y Y, LI D C, ZHAO Y G, ZHANG G L. Mapping topsoil electrical conductivity by a mixed geographically Weighted Regression Kriging: A case study in the Heihe River Basin, northwest China. Ecological Indicators, 2019, 102: 252-264.

[30] 邱孟龙, 曹小曙, 周建, 冯小龙, 高兴川. 基于GWR模型的渭北黄土旱塬粮食单产空间分异及其影响因子分析——以陕西彬县为例. 中国农业科学, 2019, 52(2): 273-284.

QIU M L, CAO X S, ZHOU J, FENG X L, GAO X C. Spatial differentiation and impact factors of grain yield per hectare in Weibei Plateau based on GWR Model: A case study of Binxian County, Shannxi. Scientia Agricultura Sinica, 2019, 52(2): 273-284. DOI: 10.3864/j.issn.0578-1752.2019.02.007. (in Chinese)

[31] 邹润彦, 周宏冀, 郭熙, 但承龙, 吕添贵, 李洪义. 环鄱阳湖区农田土壤有机碳影响因素空间分布格局分析及制图研究. 长江流域资源与环境, 2019, 28(5): 1121-1131.

ZOU R Y, ZHOU H J, GUO X, DAN C L, LV T G, LI H Y. Spatial distribution and mapping of influencing factors of farmland soil organic carbon in the Poyang Lake Region. Resources and Environment in the Yangtze Basin, 2019, 28(5): 1121-1131. DOI: 10.11870/ cjlyzyyhj201905012. (in Chinese)

[32] 杨帆, 曾维忠, 张维康, 庄天慧. 林农森林碳汇项目持续参与意愿及其影响因素. 林业科学, 2016, 52(7): 138-147.

YANG F, ZENG W Z, ZHANG W K, ZHUANG T H. Foresters’ constant participation willingness and affecting factors in forest carbon sequestration project. Scientia Silvae Sinicae, 2016, 52(7): 138-147. DOI: 10.11707/j.1001-7488.20160717. (in Chinese)

[33] 张吴平, 杨坚. 食品试验设计与统计分析. 3版. 北京: 中国农业大学出版社, 2017.

ZHANG W P, YANG J. Food Experimental Design and Statistical Analysis. 3rd ed. Beijing: China Agricultural University Press, 2017. (in Chinese)

[34] 郑新奇, 吕利娜. 地统计学（现代空间统计学）. 北京: 科学出版社, 2018.

ZHENG X Q, LV L N. Geostatistics (Modern Spatial Statistics). Beijing: Science Press, 2018. (in Chinese)

[35] 李哈滨, 王政权, 王庆成. 空间异质性定量研究理论与方法. 应用生态学报, 1998(6): 93-99.

LI H B, WANG Z Q, WANG Q C. Theory and methodology of spatial heterogeneity quantification. Chinese Journal of Applied Ecology, 1998(6): 93-99. DOI: 10.13287/j.1001-9332.1998.0136. (in Chinese)

[36] 郑帅. 温带森林土壤有机质分解对土壤含水量变化的响应[D]. 长春: 东北师范大学, 2019.

ZHENG S. Response of soil organic matter decomposition to soil water content change in temperate forests[D]. Changchun: Northeast Normal University, 2019. (in Chinese)

[37] 肖璐, 郎艺超, 夏浪, 楼昭涵, 孙楠, 黄李童, George Christakos. 基于多源数据的PM(2.5)浓度时空分布预测与制图. 环境科学, 2017, 38(12): 4913-4923.

XIAO L, LANG Y C, XIA L, LOU Z H, SUN N, HUANG L T, GEORGE C. Space-time estimations and mapping of PM2.5 fine particulates based on multi-source data. Environmental Science, 2017, 38(12): 4913-4923. DOI: 10.13227/j.hjkx.201705122. (in Chinese)

[38] 田瑞云, 王玉宽, 傅斌, 刘援. 基于DEM的地形单元多样性指数及其算法. 地理科学进展, 2013, 32(1): 121-129. (in Chinese)

TIAN R Y, WANG Y K, FU B, LIU Y. DEM-based topographic unit diversity index and its algorithm. Progress in Geography, 2013, 32(1): 121-129. DOI: 10.3724/SP.J.1033.2013.00121. (in Chinese)

[39] RILEY S J, DEGLORIA S D, ELLIOT R. A terrain ruggedness index that quantifies topographic heterogeneity. Intermountain Journal of Sciences, 1999, 5: 1-4.

Mapping of Soil Organic Matter and Its Driving Factors Study Based on MGWRK
QIAO Lei1, ZHANG WuPing2, HUANG MingJing3, WANG GuoFang1, REN Jian1

(1College of Resources and Environment, Shanxi Agricultural University, Taigu 030800, Shanxi; 2College of Software, Shanxi Agricultural University, Taigu 030800, Shanxi; 3Dryland Agriculture Research Center, Shanxi Academy of Agricultural Sciences, Taiyuan 030031)

Abstract: 【Objective】Spatial prediction is an important approach to obtain location-specific values of soil organic matter (SOM), which is an important figure of soil fertility and farmland management properly. This study was performed to compare different digital soil spatial mapping methods of SOM to get better prediction accuracy and to reveal the spatial non-stationarity characteristics of environmental covariates and the spatial scale of different environmental covariates at the same time. 【Method】In this study, the digital soil spatial mapping method was used, which was a combination of Multiscale Geographically Weighted Regression Model with simple Kriging of the residuals (MGWGK) for mapping SOM in seven towns from Jinzhong Basin. The performance of MGWRK with those of Ordinary Kriging (OK), Regression Kriging (RK), and Geographically Weighted Regression Kriging (GWRK) were compared to explore the relationship between influence factors and SOM on the influence degree and space scale.【Result】Based on the stepwise regression method, 13 indexes were selected as environmental covariables in the modeling, including aspect, slope, height, annual average precipitation, annual average temperature, gross primary productivity (GPP), evapotranspiration (ET), topographic wetness index (TWI), plan curvature, stream power index (SPI), terrain position index (TPI), terrain ruggedness index (TRI), and the annual NDVI. In the multiple linear regression (MLR), the formula had statistical significance. In the Radius index, the performance of each model was in order from good to bad: RK, OK, GWRK, MGWRK. In the mapping performance, MGWRK was close to GWRK, and both of which were better than OK method and RK method. The SOM in the study area, showed a spatial pattern of higher in middle than east and west side, among which the SOM was high in the east of the Fenhe river and the Changyuan river. The influence of aspect, annual average precipitation, annual average temperature, height, TPI and the annual NDVI on SOM in the eastern of the study area was stronger than that in the western. Whiles slope, GPP, ET, plan curvature, SPI and TRI showed opposite influence in spatial. The influence of TWI on SOM was stronger in the northern than the southern. 【Conclusion】The spatial prediction accuracy of MGWRK was 69% of RK, 71.74% of OK, and 71.17% of GWRK. MGWRK performed well in the spatial non-stationary features and the spatial visualization, which provided a reference for prediction of SOM and description spatial non-stationarity characteristics.

Key words:soil organic matter (SOM); soil mapping; multiscale geographically weighted regression kriging (MGWRK); geographically weighted regression kriging (GWRK); regression kriging (RK); non-stationarity

开放科学（资源服务）标识码（OSID）：width=42.45,height=42.45

doi: 10.3864/j.issn.0578-1752.2020.09.011

收稿日期：2019-11-14；

接受日期：2020-02-13

基金项目：山西省重点研发计划重点项目（201703D211002-2）、山西省科技攻关项目（20130311008-5）

联系方式：乔磊，E-mail：qiaolei1995@126.com。通信作者张吴平，E-mail：zwping@126.com

（责任编辑 李云霞）