土壤侵蚀造成土地资源退化，生产力下降，是危害人类生存与发展的世界性环境问题之一。深入了解土壤侵蚀规律、获取可靠的土壤侵蚀信息是开展水土保持工作的基础。核示踪技术是土壤侵蚀研究的一种重要方法，具有成本低、高效可靠等优点，其中又以¹³⁷Cs示踪技术应用最为广泛和成功^[1]。自20世纪60年代以来，在国际原子能机构(IAEA，International Atomic Energy Agency)等国际组织的大力推动下，其技术体系逐渐完善，测定的土壤侵蚀数据的准确性得到广泛验证和认可。

在一定区域内，未受扰动、无侵蚀亦无沉积发生的土壤条件下¹³⁷Cs的面积浓度代表该区域内每m²上¹³⁷Cs大气沉降的累积量，称为该区域的¹³⁷Cs本底值(¹³⁷Cs reference inventory，Bq/m²)。¹³⁷Cs示踪法研究土壤侵蚀的原理是分析研究区某一土壤剖面的¹³⁷Cs含量与¹³⁷Cs本底值之间的差异，选择合适的计算模型，定量计算得到该处的土壤侵蚀量或沉积量，因此准确的¹³⁷Cs本底值是运用¹³⁷Cs示踪技术研究土壤侵蚀的关键因素。目前，获取¹³⁷Cs本底值数据的主要手段是野外选择符合没有人为扰动、无侵蚀和沉积发生等条件的样地采集土壤样品，然后在实验室测试得到；但由于近几十年来人类活动强度大，在许多区域难以找到完全符合标准的采样地，该方法的实现存在一定难度。有学者指出，由于该方法在采样过程中普遍存在的一系列等问题，导致获取的¹³⁷Cs本底值数据的可靠性存在疑问^[2-3]。

以实测数据为基础、基于GIS的空间插值技术已广泛用于获取土壤、环境和气象等要素的高分辨率空间分布数据。¹³⁷Cs具有与上述要素类似的空间分布特性，空间插值技术亦适用于¹³⁷Cs本底值研究；但由于现阶段我国¹³⁷Cs本底值的实测样点数量偏少(近30年来，全国¹³⁷Cs本底值研究样点仅有100余个^[4])，且空间分布不均匀，并不适合通过¹³⁷Cs本底值实测数据直接进行空间插值。采用GIS的空间插值技术获取区域¹³⁷Cs本底值资料，须另辟途径。

笔者以长江上游的四川省为研究对象，尝试采用¹³⁷Cs本底值理论计算模型、主成分分析与地理加权回归克里格插值相结合的技术方法，对¹³⁷Cs本底值空间分布进行研究，探讨通过气象站点相对丰富的实测数据、基于本底值理论计算模型和GIS空间插值技术获取区域¹³⁷Cs本底值数据的合适方法，以期获取更高精度的区域¹³⁷Cs本底值资料，为在该区开展土壤侵蚀的核示踪研究、进行水土流失治理等工作提供有效的数据支持。

1 研究区概况

四川省位于我国西南，地处长江上游，E 92°21'~108°12'，N 26°03'~34°19'，区域面积48.5万km²。本区处于第1级青藏高原和第2级长江中下游平原的过渡带，区内海拔高差悬殊，地貌复杂多样，区域以山地地貌为主，全省可分为东部四川盆地、川西北高原和川西南山地等3大部分。区内季风气候明显，气候类型多样，差异显著：东部为亚热带湿润、半湿润气候，雨热同季，水热条件好，年降水量1 000~1 200 mm，但多云雾、日照时间较短；西部为高山高原高寒气候，气候垂直变化明显，总体上以寒温带气候为主，冬寒夏凉，水热不足，年降水量500~900 mm。四川省土壤类型丰富，区域分异明显，据第2次土壤普查，全省土壤类型共有25个土类、66个亚类，以水稻土、紫色土、黄壤等为主，本区是我国紫色土分布最集中的地区。

四川省是我国土壤侵蚀较严重的省区之一，土壤侵蚀分布广，面积大，危害严重。本区土壤侵蚀面积和数量约占整个长江上游总量的70%，长江上游的泥沙主要来源于地处金沙江和嘉陵江中下游的四川盆地丘陵地区。根据第1次全国水利普查公报数据，全省水蚀面积达11.4万km²，占幅员面积的23.54%，其中中度以上土壤侵蚀面积占比达到了57.63%^[5]。

2 数据来源与研究方法 2.1 数据来源

笔者采用的四川省147个气象站点的经纬度、年均降水量(观测精度为0.1 mm)等数据来自国家气象信息中心提供的“中国基本、基准和一般地面气象观测站1981—2010年累年值年值数据集”。各气象站点的空间分布如图 1所示。虽然该数据集与本文的研究对象在时间尺度上并不完全相同; 但是，该数据集具有较长的时间序列(30年)，均为地面实测数据，可信度高，且站点数量多，能保证插值的精度，满足研究要求。在参照类似研究的技术方法^[4]后，我们采用该数据集作为本研究的降水数据。

2.2 数据处理

D.E.Walling等^[6]建立的¹³⁷Cs本底值全球空间分布模型(以下简称W&H模型)是基于核爆产物大气沉降数据、建立在认识¹³⁷Cs沉降过程的基础上的理论模型，物理意义明确，目前应用最为广泛。本文各气象站点的¹³⁷Cs本底值模型计算值(以下简称¹³⁷Cs本底计算值)通过该模型计算得到，本文涉及的¹³⁷Cs本底值数据均衰变校正至2016年1月1日。

W & H模型^[6]按照10°纬度带划分计算带，对各带分别计算¹³⁷Cs本底值数据，而四川省地处N 26°~34°之间，跨越了N 20°~30°和30°~40° 2个计算带。为提高预测结果的精度，笔者对四川省进行了细分，以N30°为界划分为2个研究区域(图 1)。其中，低于N30°区域(以下简称L30区域)包含63个站点，高于N 30°区域(以下简称H30区域)包含84个站点，对这2个研究区域分别预测其¹³⁷Cs本底计算值的空间分布。

降水、地理位置及进行大气核试验的场所等是影响我国¹³⁷Cs本底值分布的主要因素^[7]。而我国的大气核试验均在新疆罗布泊地区进行，对本研究区域的¹³⁷Cs本底值分布基本没有影响；因此，笔者选择年降水量(代表降水)、经度和纬度(代表地理位置)等3个影响因素作为辅助变量，用于回归分析和空间插值。笔者采用的技术方法主要包括主成分分析和地理加权回归克里格插值。

2.3 主成分分析

主成分分析(principal component analysis, PCA)是一种应用广泛的统计分析方法，主要用于降维和消除变量间的相关性。其通过正交变换将1组存在相关性的变量转换为1组数量较少的、线性不相关的变量(称之为主成分)，在尽可能多地保留原来变量的信息的基础上，可减少变量的数量，有效降低分析问题的复杂性^[8]。

2.4 地理加权回归克里格插值

克里格插值法是一种基于地统计学的经典空间插值方法，目前已在多个领域得到广泛应用。有学者提出，通过引入与目标变量具有较好相关性的影响因素变量、构建回归克里格模型(regression Kriging, RK)，可以有效提高克里格法的插值精度^[9-11]。RK法的基本思想是建立目标变量和影响因素之间的回归方程，分离出趋势项，然后对回归残差(即随机因素)采用传统的普通克里格插值法(ordinary Kriging, OK)进行插值，最后对回归预测的趋势项和回归残差的估计值进行求和，从而得到目标变量的预测值。RK法可用以下式^[12]表示：

$ y\left( {{p_i}} \right) = \sum\limits_{j = 1}^n {{\beta _{j, \left( {{p_i}} \right)}}} {x_{j, \left( {{p_i}} \right)}} + {\varepsilon _{\left( {{p_i}} \right)}}。$

(1)

式中:y(p_i)为位置p_i处的目标变量；β_{j, (pi)}为回归系数；x_{j, (pi)}为位置p_i处的第j个自变量；ε_(pi)为位置p_i处的回归残差；n为自变量的个数。

传统的RK法中，回归分析部分大都采用基于最小二乘法(ordinary least squares, OLS)的全局回归(global regression, GR)，各回归系数为常数，即为全局回归克里格法(global regression Kriging, GRK)。由于没有考虑变量的空间非平稳性，因此在处理土壤、降水等具有较强空间非平稳性的变量时，该方法的插值精度仍受到一定的限制。

地理加权回归克里格法(geographically weighted regression Kriging, GWRK)法对GRK法进行了改进，是将GRK法中基于OLS的全局回归替换为地理加权回归的局部加权回归(local regression, LR)，其余步骤、方法保持不变的一种空间插值方法。GWRK法的核心是地理加权回归(geographically weighted regression, GWR)。GWR法对全局回归模型进行了改进，通过构建空间权重函数，对空间非平稳性进行了量化。在GWR中，回归系数不再是常数，而是空间位置的函数，位置(u_i, v_i)处的回归系数β通过下式^[12]计算得到：

$ \hat \beta \left( {{u_i}, {v_i}} \right) = {\left[ {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{W}}\left( {{u_i}, {v_i}} \right)\mathit{\boldsymbol{X}}} \right]^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{W}}\left( {{u_i}, {v_i}} \right)\mathit{\boldsymbol{Y}}。$

(2)

式中:β为回归系数矩阵；X为自变量设计矩阵(X^T为其转置矩阵)；W(u_i，v_i)为空间权重矩阵，由空间权重函数W(i)求得，其作用是定量衡量领域内不同空间位置k(k=1，2，…，n)样点的观测值对于回归点(u_i，v_i)回归系数估计的影响程度^[12](n为领域样点数)；Y为因变量矩阵。

已有研究表明，GWRK的插值结果能揭示被空间非平稳性所掩盖的一些局部变化，更真实地反映目标变量的空间变异情况，使插值精度得到进一步的提高^[13]。基于上述分析，笔者采用GWRK法进行四川省¹³⁷Cs本底计算值的插值研究。

2.5 模型精度的评价

从147个站点数据中随机选取30个站点构建验证数据集。采用其余117个站点数据建立¹³⁷Cs本底值插值模型，并求取30个验证站点的¹³⁷Cs本底值预测值，通过比较验证站点的理论计算值与插值预测值之间的平均误差(mean error, ME)、平均绝对误差(mean absolute error，MAE)、均方根误差(root mean square error, RMSE)^[13]、平均绝对百分误差(mean absolute relative error, MARE)和均方根相对误差(root mean square relative error, RMSRE)等指标评估不同插值模型的插值精度。其中，MARE和RMSRE定义如下：

$ {\rm{MARE}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left| {\frac{{{y_{a, i}} - {y_{e, i}}}}{{{y_{a, i}}}}} \right|} \times 100\% , $

(3)

$ {\rm{RMSRE}} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {\frac{{{y_{a, i}} - {y_{e, i}}}}{{{y_{a, i}}}}} \right)}^2}} } \times 100\% 。$

(4)

式中:y_{a, i}为位置i处的¹³⁷Cs本底计算值，Bq/m²；y_{e, i}为位置i处的¹³⁷Cs本底预测值，Bq/m²；n为验证数据集的站点个数。

2.6 研究软件平台

常规统计分析和PCA采用SPSS 23完成，基于OLS的GR和GWR回归分析由GWR4完成，半方差函数计算通过GS⁺9进行，克里格插值和各空间分布图的绘制由ArcGIS 10.4完成。各变量空间分布图的分辨率均为1 km×1 km。

3 结果与分析 3.1 ¹³⁷Cs本底值描述性统计分析

对研究区域各站点的¹³⁷Cs本底计算值指标进行了描述性统计分析(表 1)。分析结果表明，在2个研究区域各站点的¹³⁷Cs本底计算值在数量上存在一定的差异，H30区域和L30区域的平均值分别为1 167.42和616.39 Bq/m²，而¹³⁷Cs本底计算值的变异系数分别为8.65%和7.48%，在分异程度上基本一致。2个区域的¹³⁷Cs本底计算值数据均通过了J-B(Jarque-Bera)检验(P＞0.05)，符合正态分布。

3.2 影响因素主成分分析

表 2为计算的研究区域¹³⁷Cs本底计算值与各影响因素间的相关系数，结果显示：¹³⁷Cs本底计算值与各影响因素均存在较高的相关，与经度、降水量呈极显著正相关，与纬度呈极显著负相关。另一方面，各影响因素间也存在着不同程度的相互作用，如降水量与经度间达到了极显著正相关。因此，为了消除变量间的信息冗余性和多重共线性，更好地反映各变量间的真实关系，首先对降水量、经度和纬度等3个影响因素进行了主成分分析。

对各影响因素变量进行主成分分析后，根据方差累计贡献率≥85%的原则提取了主成分(表 3)，结果表明：对于2个研究区域，前2项主成分(F₁、F₂)的方差累计贡献率均超过90%，分别达到了95.39%和94.66%，表明前2项主成分已经表达了原始数据中绝大部分信息，采用前2项主成分即可较好地替代原数据信息；因此，后续的分析中将降水量、经度和纬度等3个原始影响因素指标转换为2个主成分变量F₁和F₂，2个主成分的空间分布如图 2所示。

3.3 回归分析

以各站点的2个主成分值为自变量，¹³⁷Cs本底计算值为因变量，分别采用基于OLS的GR和GWR 2种模型进行了回归分析。2种回归分析的结果分别如表 4和表 5所示。

方差膨胀因子(variance inflation factor，VIF)是衡量自变量间是否存在多重共线性的主要指标。由表 4的计算结果来看，对原始变量进行主成分分析后，回归分析各自变量的VIF值远小于10，这说明经过PCA处理后，各自变量间基本不存在多重共线性问题，回归分析的结果是可靠的。从GWR模型的计算结果来看，各解释变量的回归系数值不再是常数，如H30区域的F₁回归系数的取值范围介于-57.19与-34.18间，变异系数达到了33.46%，这说明¹³⁷Cs本底计算值与各解释变量间的相互关系具有显著的空间非平稳性。

对2种模型的回归效果进行了比较。表 6的分析结果表明：GWR模型的各项诊断指标均明显优于基于OLS的GR模型。如在H30区域，GWR模型的RSS与GR模型相比，大幅减少了57.33%；AICc值也从GR模型的674.61显著下降到GWR模型的639.37。这说明基于局部回归的GWR模型能够更好的探测和解释变量间的空间非平稳性，相对于GR模型，具有更好的回归效果。

3.4 回归残差的克里格插值

正态性检验的结果显示，GR模型在2个研究区域、GWR模型在L30区域的残差均通过了J-B检验，为正态分布，而GWR模型在H30区域的回归残差不符合正态分布；因此，首先对该残差数据进行5/7次方转换，使其满足进行克里格插值的条件。

在GS⁺ 9软件中，对2种模型的回归残差(包括H30区域转换后的回归残差)进行半方差函数拟合。从表 7的拟合结果可以看出，各半方差函数的R²均达到0.80，取得了较好的拟合效果。此外，各残差的半方差函数的块基比(即C₀/(C₀+C))大都小于25%，说明各残差存在较强的空间自相关性，适合使用克里格法进行空间插值。

根据上述半方差函数的拟合结果，采用OK法分别对2种模型的回归残差进行了空间插值。根据2个主成分变量及其回归系数、残差插值等数据的计算结果，分别获得了基于GRK模型和GWRK模型的¹³⁷Cs本底计算值空间分布图(图 3)。

3.5 模型预测精度评价

作为对照，同时采用OK法对各站点的¹³⁷Cs本底计算值进行了空间插值。通过构建的验证数据集分别评价了OK、GRK和GWRK等模型的预测精度。各评价指标的值越接近0，说明相应模型的预测精度越高。由表 8的分析结果可以看出，在各研究区域，模型预测精度均为GWRK＞GRK＞OK。综合考虑了降水、空间位置等多个影响因素的GRK模型和GWRK模型对¹³⁷Cs本底值的预测精度较OK法均有较大提高，而进一步探究了各影响因素空间非平稳性的GWRK模型的预测效果是这3种模型中最好的。

同时收集了已发表文献中¹³⁷Cs本底值的实测数据，与GWRK模型的预测结果进行了对比分析。从表 9的对比结果可以看出，GWRK模型的预测值普遍略低于实测值，差异幅度大致为10%。考虑到测定¹³⁷Cs活度的γ能谱仪的测试误差普遍在5%~10%之间，我们取其上限10%作为实测本底值和预测值是否存在差异的标准。由此可见，虽然预测值略低于实测值，但差异不大，实测值与本文的预测结果是基本吻合的。

3.6 四川省¹³⁷Cs本底值的空间分布

由图 3的插值结果可以看出，通过GRK和GWRK法得到的四川省¹³⁷Cs本底计算值空间分布的总趋势基本一致，均呈现北高南低、东高西低及高原低盆地高、盆周山地高盆中丘陵低的趋势，该预测结果与前人对全球¹³⁷Cs沉降规律的认识基本一致^[7]。一些局部的、细节性的¹³⁷Cs本底计算值空间分布特征在插值结果中也有所反映：在四川盆地西部边缘出现了1个局部¹³⁷Cs本底计算值峰值区，大致位于雅安市、乐山一线，¹³⁷Cs本底计算值达到了1 300~1 400 Bq/m²。出现这一现象主要是由于该区的“华西雨屏”^[19]现象造成降水量较高，而¹³⁷Cs主要来源于降水带来的湿沉降所致。此外，由于W & H计算模型的局限——模型中按照10°的纬度带划分计算区域，造成北纬30°两侧区域的¹³⁷Cs本底计算值存在“断崖式”突变。在资料精度允许的条件下，今后可以通过细分基本区域，改进模型，提高模拟结果的精度和可信度^[7]。

4 结论

以地处长江上游、土壤侵蚀特征具有典型性的四川省为研究对象，基于147个气象站点的实测数据以及W & H模型的¹³⁷Cs本底值分布模型，采用GWRK法进行空间插值，获得四川省¹³⁷Cs本底计算值的空间分布数据。初步研究结果表明：

1) 应用PCA法将经纬度、降水量等3个影响因素转换为2个主成分变量，主成分变量能表达原始数据中90%以上的信息，同时有效的消除了变量间的多重共线性，为后续的空间插值分析奠定了基础。

2) 根据PCA的分析结果，通过基于OLS的GRK和GWRK模型对四川省¹³⁷Cs本底计算值空间分布进行了预测，并评价了不同模型的预测精度。回归分析的结果表明，基于局部回归的GWRK模型比GRK模型能更好地解释¹³⁷Cs本底计算值与各影响因素间相互关系的空间非平稳性，具有更好的拟合效果；从最终的插值结果来看，GWRK模型综合考虑了降水、空间位置等影响因素，并进一步探究了各影响因素对¹³⁷Cs本底值影响的空间非平稳性，能更好地揭示影响因素和¹³⁷Cs本底值间的相互关系，其预测精度最高，且与实测结果较吻合。

3) 以实测数据为基础、基于GIS空间插值技术获取的四川省¹³⁷Cs本底计算值数据具有较高空间分辨率，并能比较准确地表达局部区域¹³⁷Cs本底值发生突变的细节信息。笔者提出的PCA+GWRK插值法可以有效的揭示¹³⁷Cs本底值的空间分布规律及不同因素对其的影响作用，对运用¹³⁷Cs示踪技术开展侵蚀泥沙研究具有一定的借鉴作用，可以为水土流失治理、土壤退化防治及水土保持工程等工作提供有效的技术支持。