CN106779061A - 一种地形平缓区土壤重金属镉空间分布预测方法 - Google Patents

Abstract

本发明公开了一种地形平缓区土壤重金属镉空间分布预测方法，首先采用径向基函数神经网络模型建立各影响因素与土壤镉含量间的非线性映射关系；针对土壤镉含量在距离分析区域不同范围内与各因素间相关性的变化，将分析区域划分为距离分析区域10km范围内和外两个区域，分别构建揭示土壤镉含量与土壤镉含量影响因素间非线性映射关系的神经网络模型；再以HASM模型对神经网络模型预测结果的残差进行模拟，得到对分析区土壤镉含量空间分布的预测结果。本发明对验证点预测结果的精度提高了5.56％‑17.65％；同时，该方法能较好地反映分析区土壤镉含量空间变化的细节信息。

Description

一种地形平缓区土壤重金属镉空间分布预测方法

技术领域

本发明属于土壤重金属含量测定技术领域，尤其涉及一种地形平缓区土壤重金属镉空间分布预测方法。

背景技术

土壤变化是目前全球环境变化的重要内容之一。近年来，随着社会经济的快速发展，化肥农药的大量施用以及工业排污的不断增加，我国的土壤环境面临的压力越来越大。其中，移动性差、难以被微生物降解的土壤重金属污染问题尤其严重，土壤重金属污染不仅会导致农作物减产或绝收，还会进入食物链或污染大气危害人畜的健康，引发慢性疾病。目前，我国已约有占耕地总面积的1/5的耕地正遭受不同程度重金属污染，这是导致粮食减产的原因之一，直接造成的经济损失高达200亿元。从广东、沈阳等城市的重点调查结果来看，污染区的癌症死亡率和癌症发病率明显高于未污染区。因此，开展土壤重金属污染防控和治理已是当前土壤环境保护工作的迫切需要。

准确掌握土壤重金属的空间变化特征是进行土壤重金属污染防治的关键和前提。受各种因素的综合影响，土壤重金属的空间分布存在明显的空间异质性。研究表明，不同成土母质、土壤类型以及不同土地利用方式[下土壤重金属含量存在明显差异。在地形平缓区，主要受区域母质、土壤类型、城市化进程、经济和交通发展以及土地利用方式的影响，土壤重金属的空间分布也具有较高的空间异质性。因此，在当前土壤环境问题日趋严峻的背景下，土壤资源的优化利用以及土壤环境的保护和治理迫切需要准确掌握区域土壤重金属含量的空间变化信息。

土壤重金属的空间变化信息通常可以通过土壤采样分析获取。但是，受采样条件、研究区范围和研究成本(人力、物力)等因素的影响，现实过程中大量采样并不可行。因此，使用科学的方法以较少的样点获得区域土壤性质的空间变化信息，就成为了一个关键问题。

传统的统计方法是通过最小值、最大值、均值和标准差等统计参数分析计算，结果可以反映区域土壤重金属含量特征的全貌，但不能描述其空间变化特征。仅基于采样点数据和自身相关性的空间插值方法，如反比距离加权法、克里格插值法等能在一定程度上避免传统统计方法的局限性，获得区域土壤重金属含量的空间分布趋势，但这类模型没有较充分地考虑各环境因素对土壤性质空间分布的影响，不能揭示复杂环境背景下土壤性质空间变化的细节信息。研究表明，引入土壤类型、河流、道路和居民点等因素作为辅助变量的预测方法其误差比仅基于采样点数据进行空间插值的方法有明显的降低，因而，将环境因子作为辅助变量进行土壤性质的空间分布模拟已是提高其模拟精度的必然途径。但在地形平缓区内，由于地形等常用环境因子差别不大，很难用于作为土壤性质空间分布模拟预测中的辅助变量。因此，合理选择和表达地形平缓区土壤重金属空间分布的影响因子，开展土壤重金属空间变异的预测方法研究，可为准确获取区域土壤重金属空间分布信息提供相应的方法参考，为经济社会发展迅速的平原区土壤环境保护提供科学依据。

综上所述，现有的土壤重金属镉空间分布预测方法不能合理选择和表达地形平缓区土壤重金属空间分布的影响因子，从而不能准确获取区域土壤重金属空间分布信息。

发明内容

本发明的目的在于提供一种地形平缓区土壤重金属镉空间分布预测方法，旨在解决现有的土壤重金属镉空间分布预测方法不能合理选择和表达地形平缓区土壤重金属空间分布的影响因子，以及未能准确捕捉各影响因子与土壤重金属间非线性关系，从而不能准确获取区域土壤重金属空间分布信息的问题。

本发明是这样实现的，一种地形平缓区土壤重金属镉空间分布预测方法，该地形平缓区土壤重金属镉空间分布预测方法包括：

首先采用径向基函数神经网络模型建立各影响因素与土壤镉含量间的非线性映射关系；针对土壤镉含量在距离分析区域不同范围内与各因素间相关性的变化，将分析区域划分为距离分析区域10km范围内和外两个区域，分别构建揭示土壤镉含量与土壤镉含量影响因素间非线性映射关系的神经网络模型；

再以HASM模型对神经网络模型预测结果的残差进行模拟，得到对分析区土壤镉含量空间分布的预测结果；

该地形平缓区土壤重金属镉空间分布预测方法记为RBF1_HASM，表达如下：

Z(x_i，k，y_j，k)＝RBFNN[F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，F_n(x_i，k，y_j，k)]+HASM(x_i，k，y_j，k)(21)

式中，Z(x_i，k，y_j，k)为土壤镉含量的预测值，RBFNN为径向基函数神经网络模型，F1～Fn为影响因素，HASM为高精度曲面模型。

进一步，该地形平缓区土壤重金属镉空间分布预测方法具体包括以下步骤：

1)将各影响因素图层栅格化，针对计算量，取分辨率为10m；以矩阵形式转入软件MATLAB中；

2)在距离分析区域10km以内和10km以外两区域内分别以各自区内建模点建立和训练神经网络模型，得到各自的优参数和对建模点的预测结果及残差；

3)以最优模型和各影响因素的空间分布数据，完成基于神经网络模型的土壤镉含量空间分布预测结果；

4)采用HASM模型，在软件MATLAB中完成对神经网络模型对建模点预测结果残差空间分布形态的模拟；

5)将运用径向基函数神经网络模型对区域镉含量值的预测结果与HASM模型对残差值的模拟结果相加，得到分析区域土壤镉含量的空间分布模拟结果；将该结果以文本格式导出，在ArcGIS中转换为栅格数据。

进一步，所述HASM模型根据微分几何学理论，空间曲面由第一类基本量和第二类基本量决定；空间曲面表达为Z＝u(x，y)，高精度曲面模型(HASM模型)基本理论被表示为：

其中E＝1+u_x ²，F＝u_xu_y，G＝1+u_y ²，

假设计算区域在x方向和y方向的最大长度分别用Lx和Ly表示，计算区域被表示为的矩形区域；h为插值步长(像元分辨率)，I+2和J+2代表x方向和y方向的栅格数，栅格(0.5h+(i-1)h，0.5h+(j-1)h)的中心点被表示为(x_i，y_j)，其中i＝0，1，…，I，I+1，j＝0，1，…，J，J+1。u(x+h，y)和u(x-h，y)用下面的泰勒展开式表示：

式(3)减式(4)得：

所以有，

对于足够小的步长h，u_x(x，y)和u_y(x，y)的有限差分方程可表示为：

式(7)加式(8)得：

所以有，

对于足够小的h，u_xx(x，y)和u_yy(x，y)的有限差分方程为：

假设是采样点{(x_i，y_j)}在曲面u上的采样值，(其中n≥0，0≤i≤I+1，0≤j≤J+1)是中心点为(x_i，y_j)的栅格值的第n次迭代值，其中 是基于采样点的插值结果。根据数值计算，由式(11)和(12)可以得出HASM的第n+1次迭代的有限差分基本表达式为：

其中n≥0；0＜i＜I+1；0＜j＜J+1；

和是HASM的边界条件；

(11-13)和(11-14)的矩阵表达形式为：

其中A₁和分别是式(11-13)的系数矩阵和右端常数项矩阵；A₂和分别是式(1-4)的系数矩阵和右端常数项矩阵；

假设为了确保采样点的真实值与采样点的估计值相等或相近，HASM的计算公式转换为以等式约束的最小二乘问题：

其中C和d分别为采样点系数矩阵和采样点的值，它们的表达式分别为：C(k，(i-1)·J+j)＝1，即第k个采样点(x_i，，y_j)的值为

为了求解最小二乘问题的方程组(11-17)，引入一个正的足够大的参数，λ是赋予采样点的权重，它决定于采样点对模拟曲面的贡献；对于足够大的λ，式(11-17)可被转化为无约束的最小二乘问题，

也就是求解：

进一步，HASM模型利用采样点数据通过一系列迭代模拟步骤对某一生态曲面的空间分布进行全局模拟，进而实现对未采样点值的预测；HASM的迭代模拟步骤具体包括：：

(a)利用采样点数据在计算区域内进行预插值，得到未采样点的预插值结果，即HASM模型的初始输入曲面；

(b)采用地理变量的初始输入曲面计算曲面的第一基本量、第二基本量以及HASM等式的系数。

(c)通过求解HASM等式得到接近真实空间分布格局的结果。

(d)重复上面迭代过程，直至达到满意精度。

进一步，径向基函数人工神经网络模型(Radial Basis Function NeuralNetwork，RBFNN)具有单隐层的三层前馈网络，其表达式如下：

f(x_i，k，y_j，k)＝RBFNN[F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，Fn(x_i，k，y_j，k)] (20)

f(x_i，k，y_j，k)为点位(x_i，y_j)处土壤镉含量值，由该点位所处环境因素共同决定；F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，Fn(x_i，k，y_j，k)为点位(x_i，k，y_j，k)处影响土壤镉含量的n个影响因素；土壤镉含量与影响因子间存在的非线性关系，采用径向基函数人工神经网络模型(RBFNN)来表达；

人工神经网络(Artificial Neural Networks，ANNs)能通过自动分析多源输入与输出间的非线性映射的关系的方式，有效解决非线性预测。

进一步，对地形平缓区土壤重金属镉空间分布预测的结果还需进行评价，该评价方法为：

以各方法对验证样点的模拟值与实测值间的相关系数(R)、平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)评价各模拟方法的模拟结果精度；其中，MAE、RMSE和MRE的计算公式如下：

式中，n为样本数；Z_obs(i)为第i个样点的实测值(mg/kg)；Z_pred(i)为预测方法对第i个样点的预测值(mg/kg)，MAE、MRE和RMSE越小，模拟方法的误差越小；相关系数R越大，预测方法的预测结果越好。

本发明提供的一种地形平缓区土壤重金属镉空间分布预测方法，利用4个误差评价指标对不同方法精度的验证结果如下所示。从表中可以得出，RBF1_HASM对验证点预测结果与实测值的相关系数最大，各项预测误差明显低于其他3种方法。与OK、MLR_HASM和RBF2_HASM相比，RBF1_HASM预测结果的平均绝对误差(MAE)分别降低了15.00％、10.53％和5.56％，均方根误差(RMSE)分别降低了17.65％、14.29％和6.67％，平均相对误差(MRE)分别降低了13.10％、11.02％和5.46％。从预测值与实测值间的相关性来看，RBF1_HASM对验证点的预测值与实测值间关系的更接近直线y＝x，其中对极大值的预测结果更明显优于其它方法。

不同方法对验证样点预测结果的精度评价

从4种预测方法对验证点预测结果的误差来还可以看，能综合考虑影响土壤镉含量的预测方法精度最高。本发明提出的方法综合考虑了土壤镉含量与影响因素间相关关系的非线性和空间非平稳性，精度明显高于其他3种方法。RBF2_HASM与RBF1_HASM相比，没有考虑土壤镉含量与各影响因素间相关关系的空间非平衡性，精度略低于RBF1_HASM。MLR_HASM方法考虑的因素进一步减少，其误差进一步增加。以上3种方法中，HASM模型基于样点残差空间自相关性对其空间分布进行模拟，进一步提高了模拟的精度。而OK法只利用了样点间的空间自相关性，在4种方法中误差最大。这表明，合理充分地利用土壤性质的空间分布性质及其与环境影响因子间的关系，选择合适的模型描述土壤性质与其影响因子间的关系，可实现区域土壤性质空间分布较高精度的预测结果。

与地形复杂丘陵区和山区不同，平原区很难利用地形因素作为土壤性质空间分布预测的辅助变量。本发明在充分分析平原特点的基础上，分析影响土壤镉含量的影响因素，合理地选择和表达了平原区预测土壤镉含量空间分布的影响因素，结合较成熟的人工神经网络模型(RBFNN)和高精度曲面模型(HASM)，实现了对成都平原研究区土壤镉含量空间分布格局的较高精度的预测，这为类似区域的研究提供了方法上的参考。

通过本发明的预测方法，获得的结论有：

(1)分析区土壤重金属镉含量范围在0.13-0.61mg.kg^-1之间，变异系数为27.23％，平均值为0.21mg.kg^-1，分别高出该区域背景值(0.16mg.kg^-1)31.25％和中国土壤镉背景值(0.097mg.kg^-1)1.16倍。根据国家土壤环境质量标准(GB15618-1995)，分析区土壤重金属镉含量总体上未超过二级标准。从区域分布来看，分析区西部靠近金马河一边的土壤重金属镉明显高于东部地区。

(2)半方差分析表明，分析区土壤镉含量的块金值和基台值之比为43.7％，说明区内土壤镉含量的空间变异性受结构性因素和随机性因素共同影响。回归分析进一步表明，与主要水系距离(金马河)、植被指数(可在一定程度上反映区域农地种植制度)、县级以下道路密度(二级道路，含：乡镇道路Y、村道C\Z)以及县级以上道路密度(一级道路，含：国道G、省道S和县道X)对分析区土壤镉含量变异的影响作用明显，分别能解释分析区土壤镉含量22.0％、12.8％、5.7％和1.4％的空间变异，即分析区土壤镉元素含量由自然因素和人为活动因素共同作用影响，这也表明自然因素对研究区土壤镉含量的影响作用大于人为因素。其中，河流(金马河)的影响范围为10km；县级以上和以下道路密度的影响范围分别为1.2-1.5km和150-200m。

(3)以水系距离、植被指数和不同等级的道路密度为辅助变量，利用径向基函数神经网络模型和高精度曲面模型构建分析区土壤镉含量空间分布预测方法(RBF1_HASM)。与3种对照方法相比，该方法对验证点预测结果的精度提高了5.56％-17.65％。同时，该方法能较好地反映分析区土壤镉含量空间变化的细节信息。

附图说明

图1是本发明实施例提供的地形平缓区土壤重金属镉空间分布预测方法流程图；

图2是本发明实施例提供的土壤重金属各数据项的半方差函数图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图对本发明的应用原理作详细描述。

本发明实施例提供的地形平缓区土壤重金属镉空间分布预测方法，该地形平缓区土壤重金属镉空间分布预测方法包括：

首先采用径向基函数神经网络模型建立各影响因素与土壤镉含量间的非线性映射关系；针对土壤镉含量在距离分析区域不同范围内与各因素间相关性的变化，将分析区域划分为距离分析区域10km范围内和外两个区域，分别构建揭示土壤镉含量与土壤镉含量影响因素间非线性映射关系的神经网络模型；

再以HASM模型对神经网络模型预测结果的残差进行模拟，得到对分析区土壤镉含量空间分布的预测结果；

该地形平缓区土壤重金属镉空间分布预测方法记为RBF1_HASM，表达如下：

Z(x_i，k，y_j，k)＝RBFNN[F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，F_n(x_i，k，y_j，k)]+HASM(x_i，k，y_j，k) (21)

式中，Z(x_i，k，y_j，k)为土壤镉含量的预测值，RBFNN为径向基函数神经网络模型，F1～Fn为影响因素，HASM为高精度曲面模型。

如图1所示，本发明实施例提供的地形平缓区土壤重金属镉空间分布预测方法具体包括以下步骤：

S101：将各影响因素图层栅格化，针对计算量，取分辨率为10m；以矩阵形式转入软件MATLAB中；

S102：在距离分析区域10km以内和10km以外两区域内分别以各自区内建模点建立和训练神经网络模型，得到各自的优参数和对建模点的预测结果及残差；

S103：以最优模型和各影响因素的空间分布数据，完成基于神经网络模型的土壤镉含量空间分布预测结果；

S104：采用HASM模型，在软件MATLAB中完成对神经网络模型对建模点预测结果残差空间分布形态的模拟；

S105：将运用径向基函数神经网络模型对区域镉含量值的预测结果与HASM模型对残差值的模拟结果相加，得到分析区域土壤镉含量的空间分布模拟结果；将该结果以文本格式导出，在ArcGIS中转换为栅格数据。

下面结合具体实施例对本发明的应用原理作进一步描述。

1材料与方法

1.1分析区概况

分析区域位于成都平原西部，东经103°66′-103°97′，北纬30°69′-30°94′之间，该区域共涉及3个县区，分别是都江堰、温江、郫县，总面积约480km²。该区域属于亚热带季风气候区，热量丰富，雨热同期，无霜期长达300d，年降水量900mm-1300mm，年平均温度15.2℃-16.7℃，年日照时数1100h-1400h。区域内成土母质主要是为第四系灰色冲积物和紫灰色冲积物。土壤类型主要为灰色及紫灰色冲积物基础上发育而成的水稻土，主要包括含4个土属。研究区属于岷江水系的自流灌溉区，处于著名水利工程都江堰的调控范围，境内最大的河流金马河和清白江分别从分析区西部和东部边缘过境。区内主要农业产品有水稻、小麦、玉米、油菜等。受成都市经济发展的影响，该区域近年来经济发展迅速，工业主要有能源、机械、医药、蓄电池、化工和食品等。

1.2土壤样点布设、采集与测定

结合分析区土壤类型、水系分布、典型种植制度以及交通道路等因素，于2013年1月至4月在兼顾均匀性的基础上以1.5km为间隔进行布点采样，采样所涉及的土地利用类型以农地为主。实际采样过程中根据情况对布设样点作一定调整，以避开水域和建设用地。以混合取样的方式在半径约为5m的范围内采集土样，取样深度为0-20cm。取样的同时详细记录每个采样点的坐标信息以及土地利用类型、海拔、土壤类型等地表环境信息。全区共采集样点339个。土样带回实验室经自然风干后，去除杂物，用木板木棒研磨，分别过2mm和0.15mm尼龙筛。土壤样品先用氢氟酸(HF)-高氯酸(HCIO₄)-硝酸(HNO₃)进行消解，再用石墨炉原子吸收法测定样品中镉含量，土壤pH值的测定则采用电位法测定。为保证数据的准确性，分析测试过程中所使用的试剂药品均为优级纯，实验所用的器皿均用10％的硝酸浸泡过夜，刷洗并用超纯水冲洗；测试过程同时以重复样和国家标准物质进行控制。

利用独立验证样点来对对各模拟方法对分析区土壤镉含量的空间分布预测结果进行精度评价。采用ArcGIS软件中的地统计模块从339个土壤样点随机抽取20％作为验证样点，其余80％作为建模样点，即建模点273个，验证点66个。建模样点和验证样点的空间分布位置。

1.3数据处理

数据的一般统计特征主要借助经典统计法完成分析区土壤重金属数据的统计特征分析，方差分析以及回归分析均在SPSS软件中完成，并采用K-S法对土壤重金属进行正态分布检验。其中，方差分析用于揭示不同定性因素对土壤重金属的影响是否存在显著的差异；回归分析用于定量表达各影响因素对土壤重金属含量空间变异的独立解释能力进行赋值。

1.4建模方法

本发明中分析区土壤重金属镉含量的空间分布预测方法的构建主要利用径向基函数人工神经网络模型(Radial Basis Function Neural Network，RBFNN模型)和高精度曲面模型(High accuracy surface modeling，HASM模型)来完成。

1.4.1 HASM模型

根据微分几何学理论，空间曲面由第一类基本量和第二类基本量决定。如果曲面可表达为Z＝u(x，y)，高精度曲面模型(HASM模型)基本理论可被表示为：

其中E＝1+u_x ²，F＝u_xu_y，G＝1+u_y ²，

假设计算区域在x方向和y方向的最大长度分别用Lx和Ly表示，计算区域可被表示为的矩形区域。h为插值步长(像元分辨率)，I+2和J+2代表x方向和y方向的栅格数，栅格(0.5h+(i-1)h，0.5h+(j-1)h)的中心点被表示为(x_i，y_j)，其中i＝0，1，…，I，I+1，j＝0，1，…，J，J+1。u(x+h，y)和u(x-h，y)可用下面的泰勒展开式表示：

式(3)减式(4)得：

所以有，

对于足够小的步长h，u_x(x，y)和u_y(x，y)的有限差分方程可表示为：

式(7)加式(8)得：

所以有，

对于足够小的h，u_xx(x，y)和u_yy(x，y)的有限差分方程为：

假设是采样点{(x_i，y_j)}在曲面u上的采样值，(其中n≥0，0≤i≤I+1，0≤j≤J+1)是中心点为(x_i，y_j)的栅格值的第n次迭代值，其中 是基于采样点的插值结果。根据数值计算，由式(11)和(12)可以得出HASM的第n+1次迭代的有限差分基本表达式为：

其中n≥0；0＜i＜I+1；0＜j＜J+1；

和是HASM的边界条件。

(11-13)和(11-14)的矩阵表达形式为：

其中A₁和分别是式(11-13)的系数矩阵和右端常数项矩阵；A₂和分别是式(1-4)的系数矩阵和右端常数项矩阵。

假设为了确保采样点的真实值与采样点的估计值相等或相近，HASM的计算公式可转换为以等式约束的最小二乘问题：

其中C和d分别为采样点系数矩阵和采样点的值，它们的表达式分别为：C(k，(i-1)·J+j)＝1，即第k个采样点(x_i，y_i)的值为

为了求解最小二乘问题的方程组(11-17)，引入一个正的足够大的参数，λ是赋予采样点的权重，它决定于采样点对模拟曲面的贡献^[31]。对于足够大的λ，式(11-17)可被转化为无约束的最小二乘问题，

也就是求解：

HASM模型可利用采样点数据通过一系列迭代模拟步骤对某一生态曲面的空间分布进行全局模拟，进而实现对未采样点值的预测。HASM的迭代模拟步骤总结如下：

(1)利用采样点数据在计算区域内进行预插值，得到未采样点的预插值结果，即HASM模型的初始输入曲面。

(2)采用地理变量的初始输入曲面计算曲面的第一基本量、第二基本量以及HASM等式的系数。

(3)通过求解HASM等式得到接近真实空间分布格局的结果。

(4)重复上面迭代过程，直至达到满意精度。

1.4.2人工神经网络模型

人工神经网络(Artificial Neural Networks，ANNs)能通过自动分析多源输入与输出间的非线性映射的关系的方式，有效解决非线性预测问题。径向基函数人工神经网络模型(Radial Basis Function Neural Network，RBFNN)是一种具有单隐层的三层前馈网络，在参数设置、学习速度等方面拥有诸多优点，在土壤性质空间分布预测分析中得到了广泛应用，取得了较好的预测效果。其表达式如下：

f(x_i，k，y_j，k)＝RBFNN[F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，Fn(x_i，k，y_j，k)] (20)

f(x_i，k，y_j，k)为点位(x_i，y_j)处土壤镉含量值，由该点位所处环境因素共同决定。F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，Fn(x_i，k，y_j，k)为点位(x_i，k，y_j，k)处影响土壤镉含量的n个影响因素。土壤镉含量与影响因子间可能存在的非线性关系，采用径向基函数人工神经网络模型(RBFNN)来表达。

1.4.3平原区重土壤镉含量空间分布预测模型构建

首先采用径向基函数神经网络模型建立各影响因素与土壤镉含量间的非线性映射关系；考虑到土壤镉含量在距离金马河不同范围内与各因素间相关性的变化，将分析区域划分为距离金马河10km范围内和外两个区域，分别构建揭示土壤镉含量与其影响因素间非线性映射关系的神经网络模型；再以HASM模型对神经网络模型预测结果的残差进行模拟，得到对分析区土壤镉含量空间分布的预测结果。该方法记为RBF1_HASM，其表达如下：

Z(x_i，k，y_j，k)＝RBFNN[F1(x_i，k，y_j，k)，F2(x_i，k，y_j，k)，…，F_n(x_i，k，y_j，k)]+HASM(x_i，k，y_j，k) (21)

式中，Z(x_i，k，y_j，k)为土壤镉含量的预测值，RBFNN为径向基函数神经网络模型，F1～Fn为影响因素，HASM为高精度曲面模型。

该方法的计算过程如下：

1)将各影响因素图层栅格化，考虑到计算量，取分辨率为10m；以矩阵形式转入软件MATLAB中。

2)在距离金马河10km以内和10km以外两区域内分别以各自区内建模点建立和训练神经网络模型，得到各自的优参数和对建模点的预测结果及残差。

3)以最优模型和各影响因素的空间分布数据，完成基于神经网络模型的土壤镉含量空间分布预测结果。

4)采用HASM模型，在软件MATLAB中完成对神经网络模型对建模点预测结果残差空间分布形态的模拟。

5)将运用径向基函数神经网络模型对区域镉含量值的预测结果与HASM模型对残差值的模拟结果相加，即可得到研究区土壤镉含量的空间分布模拟结果。将该结果以文本格式导出，在ArcGIS中转换为栅格数据。

1.4.4对照方法

1)RBF2_HASM

该方法与RBF1_HASM不同之处在于，在全区只用一个神经网络模型来建立研究区土壤镉含量与其影响因素间非线性映射关系，即没有考虑土壤镉含量与各影响因素间相关关系的空间非平衡性。其计算过程与RBF1_HASM相同。该方法记为RBF2_HASM。

2)MLR_HASM

该方法采用多元线性回归模型来建立分析区全域内土壤镉含量与其影响因素间的映射关系。回归模型是目前土壤性质空间分布预测研究中融合环境辅助变量最常用的方法之一。该方法与RBF1_HASM不同之处在于，没有考虑土壤镉含量与各影响因素间相关关系的非线性和空间非平衡性。方法计算过程中，回归模型预测部在ArcGIS中完成，HASM模型模拟部分仍然在软件MATLAB中进行。方法记为MLR_HASM。

3)OK

普通克里格法(OrdinaryKriging，OK)是一种最优线性无偏估计插值方法。它是以区域化变量为基础基础，应用变异函数理论，根据已知采样点的属性值在其邻近区域内对未采样区域的属性值进行无偏最优估计，具体的理论与方法介绍在许多文献中均有详细叙述。该方法是目前土壤性质空间分布研究中最常用的方法。本研究过程中，以建模点在GS⁺中完成半方差分析，获得普通克里格法插值所需要的最优参数，再在ArcMap10.1中完成空间插值分析。与RBF1_HASM相比，该方法仅依据采样点间的空间自相关性进行内插，没有考虑其他与土壤性质相关的影响因素。

1.4.5预测结果评价方法

以各方法对验证样点的模拟值与实测值间的相关系数(R)、平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)评价各模拟方法的模拟结果精度。其中，MAE、RMSE和MRE的计算公式如下：

式中，n为样本数；Z_obs(i)为第i个样点的实测值(mg/kg)；Z_pred(i)为预测方法对第i个样点的预测值(mg/kg)，MAE、MRE和RMSE越小，模拟方法的误差越小。而相关系数R越大，预测方法的预测结果越好。

2结果与分析

2.1土壤重金属镉含量统计特征

统计结果表明(表1)，分析区土壤pH值在4.18-8.42之间，平均为6.41，总体上表现为弱酸性土壤。基于339个土样数据统计显示，研究区土壤重金属镉含量范围在0.13-0.61mg.kg^-1之间，变异系数为27.23％，平均值为0.21mg.kg^-1，根据国家土壤环境质量标准(GB15618-1995)，分析区土壤重金属镉含量总体上未超过二级标准，但高提出的该区域背景值(0.16mg.kg^-1)31.25％和中国土壤镉背景值(0.097mg.kg^-1)1.16倍。该结果也明显低于成都平原其他区域。如李等人的分析表明金马河西侧崇州地区土壤重金属镉含量均为0.33mg.kg^-1；王等人2010年在成都平原西部6个县区得到的研究结果为0.32mg.kg^-1。

从数据分布类型来看(表1)，分析区土壤镉呈明显的偏态分布(K-S检验，P＜0.01)。由于本发明中方差分析和回归分析要求变量为正态分布，同时也为避免数据出现负值，后续分析过程中将土壤镉含量值扩大10倍再取对数值，对数转换后的值符合正态分布(K-S检验，P＝0.192)，因此将该对数转换值用于本发明后续的方差分析和回归分析。

表1 研究区农地土壤pH值和镉含量统计特征

2.2土壤重金属镉含量空间结构

统计分析只能反映出分析区土壤镉含量变化的总体特征，不能较好揭示区域土壤镉含量空间变化的结构性，需进一步利用地统计学方法中的半方差模型对分析区土壤镉含量空间分布的结构特征进行分析。从半方差拟合结果来看，分析区土壤镉含量符合高斯模型，拟合决定系数为0.968(表2、图2)，拟合程度较高，可以用来反映分析区土壤镉元素的空间变异结构特征。半方差分析中，土壤性质的块金值与基台值之比(块金系数)＜25％时，空间变量具有强烈的空间自相关性，该变量主要受结构性因素的影响；若该比值在25％-75％之间，表明该变量具有中等程度的空间自相关性，受结构性因素和随机性因素共同作用；若比值＞75％，说明该变量的空间自相关性很弱，主要受随机因素控制。其中，结构性因素主要包括土壤在发育形成过程中的母质、地形因素、及成土后所属的土壤类型等，随机性因素包括人类的施肥和耕作活动等因素。结构因素使得土壤属性具有空间自相关性，而随机因素则减弱土壤属性的空间自相关性，增大其异质性。

从半方差函数拟合结果的参数来看，分析区土壤镉元素的块金值和基台值之比为43.7％，说明区内土壤镉含量的空间分布同时受结构因素和随机因素共同影响。从土壤镉含量块金系数值的大小来看，受结构性因素的影响程度高于随机因素的影响，与邻近区域的研究结果一致。统计分析表明，成都平原及周边区域的土壤镉含量在不同母质、土壤类型、地形及耕作施肥背景下差异显著。但各影响因素对本发明区土壤镉含量空间变异的影响程度大小并不明确，需进一步进行探讨和分析。

表2 研究区土壤重金属镉含量的半方差函数及其拟合参数

2.3土壤重金属镉的影响因素分析

根据已有数据的可获取性，以分析区水系分布和土壤类型作为结构因素，交通道路以及能在一定程度上反映农地种植制度的植被指数作为随机因素，分别探讨不同因素影响对分析区土壤重金属镉含量的影响。

2.3.1水系分布

分别以左侧和右侧的金马河和清白江分别为一级河流，在ArcGIS软件中生成两条河流的缓冲距离，提取各样点分别与金马河和清白江的距离；同时以过境的走马河、江安河、柏条河和徐堰河为二级河流，生成二级河流的缓冲距离，并提取各样点与二级河流的距离，分析不同的河流对分析区土壤重金属镉的影响。一级河流对土壤镉含量有较大的影响，二级河流对土壤镉含量的影响不大。在靠近金马河区域内，土壤镉含量与距金马河距离成反比，而土壤镉含量与距清白江距离成正比。其中，在金马河区域内10km范围内土壤镉含量下降趋势尤其显著。距离水系10km范围内土壤镉含量下降趋势显著，其中在5km范围内土壤镉含量明显较高，样点土壤镉含量超过二级标准的比例远大于5km范围内之外。

成都平原西部土壤镉含量的背景平均值达到0.258mg.kg^-1，分析区西北部山区土壤镉含量背景平均值为0.376mg.kg^-1，其中各出露的地层中二叠系(P)、石炭系(C)、泥盆系(D)和震旦系(Z)地层出露区域镉的背景平均含量更高达0.659mg.kg^-1，这些因素导致了由岷江水系形成的都江堰至崇州市元通镇的冲积扇群及湔江冲积扇群上土壤中镉出现明显的富集特征。长江流域范围内四川段各河流中，岷江水系的一级阶地和水系沉积物中镉含量值分别达到了0.27和0.53mg.kg^-1，显著高于其他支流。

本发明的分析区主要处于岷江水系的金马河流域区，该区域母质直接来源于北部的二叠系(P2)和泥盆系(D2g)地层，土壤镉的较高含量值主要受到了来自金马河上游流经区域出露地层的影响，而冲积平原区距离河流的远近可反映案例内成土过程中母质和水文特征的差异性，进而对土壤性质产生影响。因此，与金马河的距离对区内土壤镉含量有较大的影响，影响范围为10km左右。这一影响范围大于金马河西侧的4km，这可能是因为金马河距西侧同等级的西河约11km，而西河所带来的冲积母质同样来自成都西北部山区镉背景值较高的地层，金马河和西河中间地带区域距两条河的距离在4-5km左右。

土壤镉含量在靠近青白江一侧明显较低，这可能与青白江的流域范围以及青白江和柏条河流经的区域地质背景有关。在青白江和柏条河流经的天马镇、丽春镇和桂花出现低山，出露地层为白垩系灌口组(K2g)和夹关组砂岩(K2j)。已有研究表明，成都平原区紫色冲积物发育而来的土壤镉含量在各母质中最低^[25]。

表3 距金马河不同距离段土壤重金属镉含量统计特征

2.4.2土壤类型

从不同土壤类型镉含量的统计结果来看(表4)，不同土壤类型镉含量差异不大(方差分析，P＞0.05)。其中，脱潜潮田的镉平均含量相对较高，为0.25mg.kg^-1，该土壤类型所属的样点镉含量超过二级标准的比例也远高于其他土壤类型；潜育潮田、渗育灰潮田和潴育灰潮田发育的土壤镉平均含量相差不大，二级超标率分别为8.7％、7.53％和2.74％。这可能与区域的土壤类型分布有关，区内主要为渗育灰潮田，其余类型分布较少。

表4 研究区不同土属重金属镉含量特征

2.4.3植被系数

分析区土壤镉含量与植被指数MODIS NDVI的关系，两者间的相关系数为-0.362(P＜0.01)，表现为极显著的负相关关系。不同种植方式下农地上肥料农药施用等农作管理不同。区内农地利用主要以水稻-油菜轮作和水稻-小麦轮作为主。不同的种植制度中水稻与油菜轮作下土壤镉含量低于水稻-小麦轮作。成都地区油菜体含镉量大于小麦，而小麦生物量明显小于油菜的生物量。分析区内农地在收获油菜时，油菜杆一般会收走；而收获小麦时一般在收割的同时将小麦杆粉碎还田。本发明对植被指数的处理可以在一定程度上反映区域农地的种植制度。高NDVI主要反映的是水稻-油菜轮作。因此，分析区土壤镉含量与植被指数表现出极显著的负相关关系。

2.4.4交通道路

以分析区县级以上交通道路为一级道路(G、S、X)、县级以下交通道路二级道路(Y、C、Z)，生成分析区两种不同等级交通公路的缓冲距离，分析交通公路对土壤镉含量的影响。在距离一级道路1.2-1.5km范围内，土壤镉含量值随着与交通道路距离的增大呈减少趋势；在距离二级道路150-200m范围内，土壤镉含量随着与交通道路距离的增大而略有下降。这表明，县级以上交通道路对土壤镉的影响范围在1.2-1.5km左右，县级以下交通道路的影响为150-200m。

根据土壤镉含量与道路缓冲距离的分析结果，在ArcGIS中分别计算分析区500×500m格网和100×100m格网中一级道路和二级道路的长度，并转换为每平方公里范围内一级道路和二级道路的长度，得到分析区1km²范围内不同等级道路密度，进而分析研究区土壤镉含量与道路密度间的关系。土壤镉含量随着交通道路密度的增加而增加的趋势。相关分析结果表明，分析区土壤镉含量与一级道路(G、S、X)密度的相关系数为0.138，表现为显著的正相关关系(P＜0.05)，与二级道路(Y、C、Z)密度的相关系数为0.241，表现为极显著的正相关关系(P＜0.01)。与二级道路密度的相关性高于一级道路，这是因为尽管一级道路影响范围大和交通流量大，但密度远小于二级道路。

2.4.5各因素综合影响分析

以对分析区土壤镉含量有显著影响的水系距离(与金马河距离)、植被指数和交通密度作为自变量进行回归分析，以定量揭示不同因素对研究区内土壤镉含量空间分布的影响作用。结果表明，各回归方程都达到了显著以上水平(表5)，说明这些因素对区域内土壤镉含量空间变异有着显著的影响。

从回归分析结果可以看出(表5)，各影响因素对分析区土壤镉含量变异的独立解释能力大小为：与水系的距离＞NDVI＞县级以下道路密度(二级道路)＞县级以上道路密度(一级道路)。其中，土壤镉含量与可体现自然地质过程因素的与水系距离这一因子呈对数关系，为非线性关系；

对分析区土壤镉含量空间变异的影响程度最大，达到22.0％，远大于人为活动因素中的道路密度和可体现农地种植方式的NDVI(表5)。这主要是因为分析区处于冲积平原区，水系在土壤的形成和发过程中具有较重要的作用，与水系的距离可反映了成土母质在区域内的分选过程和水文变化。同时这一区域的冲积母质主要来源于西北部山区镉背景含量较高的地层。

可反映农地种植活动的植被指数NDVI对土壤镉含量变异的独立解释能力为12.8％，高于道路密度的影响，这是因为农地种植活动是在整个区域内进行，其施肥耕作活动是区域性，而交通公路则主要对其道路两侧的土壤产生影响，影响范围相对较小。由于分析区县级以下道路密度较大，交通道路也能解释一定的土壤镉含量空间变异。

上述分析结果表明，分析区土壤镉含量受自然因素和人为活动因素共同影响；而影响因素作用程度的定量分析结果表明，结构性因素自然地质过程(主要是水系分布)对分析区土壤镉含量的影响作用总体上高于土地利用方式和交通等人为活动因素，与半方差分析结果一致。

表5 不同影响因素下土壤镉含量的回归分析结果

注：县级以上道路含(高速G、省道S、县道X)，县级以下道路(乡镇道路Y和村道C)

考虑到水系的影响范围主要在距离金马河10km范围内，本发明采用相关分析分别计算分析区内距离金马河10km以内和10km以外土壤镉含量与各因素间的相关性。结果表明(表6)，距离金马河10km以内，土壤镉含量与水系距离、植被指数和县级以下道路密度呈现极显著的相关关系；距离金马河10km之外，与县级以上道路密度和植被指数呈现极显著和显著的相关关系。这说明，分析区土壤镉含量与各因素间的相关性存在明显的空间非平衡性，构建该区域土壤镉含量空间预测模型时需考虑土壤镉与其影响因素间关系的空间非平衡性。

表6 距金马河不同距离段内各影响因素对土壤镉含量的相关系数

注：^**和^*分别表示显著水平为0.01和0.05。

2.4土壤重金属镉含量空间分布预测结果及精度验证

2.4.1不同方法预测结果

4种预测方法预测结果。4种方法得到的成都平原研究区土壤镉含量的空间分布趋势基本一致，低值区主要分布于中部和东部地区，而高值区则分布于研究区西部地区。4种方法得到的高值区和低值区所分布的位置基本一致。

但不同模拟方法对分析区土壤镉含量空间分布格局的预测结果也存在着明显的差异。4种方法中，普通克里法(OK)具有明显的平滑效应，其预测结果最为光滑，高值区和低值区呈明显的块状分布，仅能反映土壤镉的空间分布趋势，如西南部的高值出呈明显的块状分布，可以看该结果受个别高值样点的影响较大。融合环境影响因素的3种预测方法MLR_HASM、RBF1_HASM和RBF2_HASM明显在反映土壤镉空间变异更多的细节信息，如交通等因素的影响。

2.6.2精度评价

利用4个误差评价指标对不同方法精度的验证结果如表7所示。从表中可以得出，RBF1_HASM对验证点预测结果与实测值的相关系数最大，各项预测误差明显低于其他3种方法。与OK、MLR_HASM和RBF2_HASM相比，RBF1_HASM预测结果的平均绝对误差(MAE)分别降低了15.00％、10.53％和5.56％，均方根误差(RMSE)分别降低了17.65％、14.29％和6.67％，平均相对误差(MRE)分别降低了13.10％、11.02％和5.46％。从预测值与实测值间的相关性来看，RBF1_HASM对验证点的预测值与实测值间关系的更接近直线y＝x，其中对极大值的预测结果更明显优于其它方法。

表7 不同方法对验证样点预测结果的精度评价

从4种预测方法对验证点预测结果的误差来还可以看，能综合考虑影响土壤镉含量的预测方法精度最高。本发明提出的方法综合考虑了土壤镉含量与影响因素间相关关系的非线性和空间非平稳性，精度明显高于其他3种方法。RBF2_HASM与RBF1_HASM相比，没有考虑土壤镉含量与各影响因素间相关关系的空间非平衡性，精度略低于RBF1_HASM。MLR_HASM方法考虑的因素进一步减少，其误差进一步增加。以上3种方法中，HASM模型基于样点残差空间自相关性对其空间分布进行模拟，进一步提高了模拟的精度。而OK法只利用了样点间的空间自相关性，在4种方法中误差最大。这表明，合理充分地利用土壤性质的空间分布性质及其与环境影响因子间的关系，选择合适的模型描述土壤性质与其影响因子间的关系，可实现区域土壤性质空间分布较高精度的预测结果。

与地形复杂丘陵区和山区不同，平原区很难利用地形因素作为土壤性质空间分布预测的辅助变量。本发明在充分分析平原特点的基础上，分析影响土壤镉含量的影响因素，合理地选择和表达了平原区预测土壤镉含量空间分布的影响因素，结合较成熟的人工神经网络模型(RBFNN)和高精度曲面模型(HASM)，实现了对成都平原研究区土壤镉含量空间分布格局的较高精度的预测，这为类似区域的研究提供了方法上的参考。

3结论

(1)分析区土壤重金属镉含量范围在0.13-0.61mg.kg^-1之间，变异系数为27.23％，平均值为0.21mg.kg^-1，分别高出该区域背景值(0.16mg.kg^-1)31.25％和中国土壤镉背景值(0.097mg.kg^-1)1.16倍。根据国家土壤环境质量标准(GB15618-1995)，分析区土壤重金属镉含量总体上未超过二级标准。从区域分布来看，分析区西部靠近金马河一边的土壤重金属镉明显高于东部地区。

(2)半方差分析表明，分析区土壤镉含量的块金值和基台值之比为43.7％，说明区内土壤镉含量的空间变异性受结构性因素和随机性因素共同影响。回归分析进一步表明，与主要水系距离(金马河)、植被指数(可在一定程度上反映区域农地种植制度)、县级以下道路密度(二级道路，含：乡镇道路Y、村道C\Z)以及县级以上道路密度(一级道路，含：国道G、省道S和县道X)对分析区土壤镉含量变异的影响作用明显，分别能解释研究区土壤镉含量22.0％、12.8％、5.7％和1.4％的空间变异，即分析区土壤镉元素含量由自然因素和人为活动因素共同作用影响，这也表明自然因素对研究区土壤镉含量的影响作用大于人为因素。其中，河流(金马河)的影响范围为10km；县级以上和以下道路密度的影响范围分别为1.2-1.5km和150-200m。

(3)以水系距离、植被指数和不同等级的道路密度为辅助变量，利用径向基函数神经网络模型和高精度曲面模型构建分析区土壤镉含量空间分布预测方法(RBF1_HASM)。与3种对照方法相比，该方法对验证点预测结果的精度提高了5.56％-17.65％。同时，该方法能较好地反映分析区土壤镉含量空间变化的细节信息。

本发明在充分分析平原区土壤镉含量影响因素的基础上，合理选择和表达了平原区预测土壤镉含量空间分布的影响因素，结合人工神经网络模型和高精度曲面模型，实现了对平原研究区土壤镉含量空间分布格局的较高精度的预测，为类似区域的研究提供了方法上的参考。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种地形平缓区土壤重金属镉空间分布预测方法，其特征在于，该地形平缓区土壤重金属镉空间分布预测方法包括：

首先采用径向基函数神经网络模型建立各影响因素与土壤镉含量间的非线性映射关系；针对土壤镉含量在距离分析区域不同范围内与各因素间相关性的变化，将分析区域划分为距离分析区域10km范围内和外两个区域，分别构建揭示土壤镉含量与土壤镉含量影响因素间非线性映射关系的径向基函数人工神经网络模型；

再以HASM模型对神经网络模型预测结果的残差进行模拟，得到对分析区土壤镉含量空间分布的预测结果；

该地形平缓区土壤重金属镉空间分布预测方法记为RBF1_HASM，表达如下：

Z(x_i,k,y_j,k)＝RBFNN[F1(x_i,k,y_j,k),F2(x_i,k,y_j,k),...,F_n(x_i,k,y_j,k)]+HASM(x_i,k,y_j,k)；

式中，Z(x_i,k,y_j,k)为土壤镉含量的预测值，RBFNN为径向基函数神经网络模型，F1～Fn为影响因素，HASM为高精度曲面模型。

2.如权利要求1所述的地形平缓区土壤重金属镉空间分布预测方法，其特征在于，该地形平缓区土壤重金属镉空间分布预测方法具体包括以下步骤：

1)将各影响因素图层栅格化，针对计算量，取分辨率为10m；以矩阵形式转入软件MATLAB中；

2)在距离分析区域10km以内和10km以外两区域内分别以各自区内建模点建立和训练神经网络模型，得到各自的优参数和对建模点的预测结果及残差；

3)以最优模型和各影响因素的空间分布数据，完成基于神经网络模型的土壤镉含量空间分布预测结果；

4)采用HASM模型，在软件MATLAB中完成对神经网络模型对建模点预测结果残差空间分布形态的模拟；

5)将运用径向基函数神经网络模型对区域镉含量值的预测结果与HASM模型对残差值的模拟结果相加，得到分析区域土壤镉含量的空间分布模拟结果；将该结果以文本格式导出，在ArcGIS中转换为栅格数据。

3.如权利要求1～2任意一项所述的地形平缓区土壤重金属镉空间分布预测方法，其特征在于，所述HASM模型根据微分几何学理论，空间曲面由第一类基本量和第二类基本量决定；若空间曲面表达为Z＝u(x,y)，则其第一类基本量表达为，

$\left\{\begin{array}{c}E=1+u_x^2\\ F=u_x{u}_y\\ G=1+u_y^2\end{array}\right.---\left(1\right)$

第二类基本量表达为，

$\left\{\begin{array}{c}L=\frac{u_{xx}}{\sqrt{1+u_x^2+u_y^2}}\\ M=\frac{u_{xy}}{\sqrt{1+u_x^2+u_y^2}}\\ N=\frac{u_{yy}}{\sqrt{1+u_x^2+u_y^2}}\end{array}\right.---\left(2\right)$

高精度曲面模型基本理论被表示为：

$\left\{\begin{array}{c}{u}_{xx}={\mathrm{\Γ}}_{11}^1{u}_x+{\mathrm{\Γ}}_{11}^2{u}_y+\frac{L}{\sqrt{EG-F^2}}\\ u_{xy}={\mathrm{\Γ}}_{12}^1{u}_x+{\mathrm{\Γ}}_{12}^2{u}_y+\frac{M}{\sqrt{EG-F^2}}\\ u_{yy}={\mathrm{\Γ}}_{22}^1{u}_x+{\mathrm{\Γ}}_{22}^2{u}_y+\frac{N}{\sqrt{EG-F^2}}\end{array}\right.---\left(3\right)$

其中

假设计算区域在x方向和y方向的最大长度分别用Lx和Ly表示，计算区域被表示为的矩形区域；h为插值步长，I+2和J+2代表x方向和y方向的栅格数，栅格(0.5h+(i-1)h,0.5h+(j-1)h)的中心点被表示为(x_i,y_j)，其中i＝0,1,…,I,I+1，j＝0,1,…,J,J+1；u(x+h,y)和u(x-h,y)用下面的泰勒展开式表示：

$u(x+h,y)=u(x,y)+h\frac{\∂u(x,y)}{\∂x}+\frac{h^2}{2!}\frac{{\∂}^{2}u(x,y)}{\∂x^2}+\frac{h^3}{3!}\frac{{\∂}^{3}u(x,y)}{\∂x^3}+O\left(h^4\right)---\left(4\right)$

$u(x-h,y)=u(x,y)-h\frac{\∂u(x,y)}{\∂x}+\frac{h^2}{2!}\frac{{\∂}^{2}u(x,y)}{\∂x^2}-\frac{h^3}{3!}\frac{{\∂}^{3}u(x,y)}{\∂x^3}+O\left(h^4\right)---\left(5\right)$

式(3)减式(4)得：

$u(x+h,y)-u(x-h,y)=2h\frac{\∂u(x,y)}{\∂x}+\frac{2h^3}{3!}\frac{{\∂}^{3}u(x,y)}{\∂x^3}+O\left(h^5\right)---\left(6\right)$

所以有，

$u_x(x,y)=\frac{\∂u(x,y)}{\∂x}=\frac{u(x+h,y)-u(x-h,y)}{2h}+\frac{h^2}{3!}\frac{{\∂}^{3}u(x,y)}{\∂x^3}+O{\left(h\right)}^4---\left(7\right)$

对于足够小的步长h，u_x(x,y)和u_y(x,y)的有限差分方程可表示为：

$u_x(x,y)\≈\frac{u(x+h,y)-u(x-h,y)}{2h}---\left(8\right)$

$u_y(x,y)\≈\frac{u(x,y+h)-u(x,y-h)}{2h}---\left(9\right)$

式(7)加式(8)得：

$u(x+h,y)+u(x-h,y)=2u(x,y)+\frac{2h^2}{2!}\frac{{\∂}^{2}u(x,y)}{\∂x^2}+O\left(h^4\right)---\left(10\right)$

所以有，

$u_{xx}(x,y)=\frac{{\∂}^{2}u(x,y)}{\∂x^2}=\frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2}+O\left(h^2\right)---\left(11\right)$

对于足够小的h，u_xx(x,y)和u_yy(x,y)的有限差分方程为：

$u_{xx}(x,y)\≈\frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2}---\left(12\right)$

$u_{yy}(x,y)\≈\frac{u(x,y+h)-2u(x,y)+u(x,y-h)}{h^2}---\left(13\right)$

假设是采样点{(x_i，y_j)}在曲面u上的采样值，(其中n≥0，0≤i≤I+1，0≤j≤J+1)是中心点为(x_i,y_j)的栅格值的第n次迭代值，其中 是基于采样点的插值结果；根据数值计算，由式(12)和(13)得出HASM的第n+1次迭代的有限差分基本表达式为：

$\frac{u_{i+1,j}^{n+1}-2u_{i,j}^{n+1}+u_{i-1,j}^{n+1}}{h^2}={\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}^n\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2h}+{\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}^n\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}+\frac{L_{i,j}^n}{\sqrt{E_{i,j}^n+G_{i,j}^n-1}}---\left(14\right)$

$\frac{u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1}}{h^2}={\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}^n\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2h}+{\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}^n\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}+\frac{N_{i,j}^n}{\sqrt{E_{i,j}^n+G_{i,j}^n-1}}---\left(15\right)$

其中n≥0；0＜i＜I+1；0＜j＜J+1；

$F_{i,j}^n=\left(\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2h}\right)\left(\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}\right);G_{i,j}^n=1+{\left(\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}\right)}^2;$

$L_{i,j}^n=\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{h^2\sqrt{1+{\left(\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2h}\right)}^2+{\left(\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}\right)}^2}};$

$N_{i,j}^n=\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{h^2\sqrt{1+{\left(\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2h}\right)}^2+{\left(\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}\right)}^2}};$

${\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}^n=\frac{G_{i,j}^n(E_{i+1,j}^n-E_{i,j}^n)-2F_{i,j}^n(F_{i+1,j}^n-F_{i-1,j}^n)+F_{i,j}^n(E_{i,j+1}^n-E_{i,j-1}^n)}{4(E_{i,j}^n{G}_{i,j}^n-{\left(F_{i,j}^n\right)}^2)h};$

${\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}^n=\frac{2G_{i,j}^n(F_{i,j+1}^n-F_{i,j-1}^n)-G_{i,j}^n(G_{i+1,j}^n-G_{i-1,j}^n)-F_{i,j}^n(G_{i,j+1}^n-G_{i,j-1}^n)}{4(E_{i,j}^n{G}_{i,j}^n-{\left(F_{i,j}^n\right)}^2)h};$

${\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}^n=\frac{2E_{i,j}^n(F_{i+1,j}^n-F_{i-1,j}^n)-E_{i,j}^n(E_{i,j+1}^n-E_{i,j-1}^n)-F_{i,j}^n(E_{i+1,j}^n-E_{i-1,j}^n)}{4(E_{i,j}^n{G}_{i,j}-{\left(F_{i,j}^n\right)}^2)h};$

${\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}^n=\frac{E_{i,j}^n(G_{i,j+1}^n-G_{i,j-1}^n)-2F_{i,j}^n(F_{i,j+1}^n-F_{i,j-1}^n)+F_{i,j}^n(G_{i+1,j}^n-G_{i-1,j}^n)}{4(E_{i,j}^n{G}_{i,j}^n-{\left(F_{i,j}^n\right)}^2)h};u_{0,j}^{n+1}=u_{0,j}^0,(0\≤j\≤J+1);$

和是HASM的边界条件；

(11-13)和(11-14)的矩阵表达形式为：

$A_1{U}^{n+1}=b_1^n---\left(16\right)$

$A_2{U}^{n+1}=b_2^n---\left(17\right)$

其中A₁和分别是式(11-13)的系数矩阵和右端常数项矩阵；A₂和分别是式(1-4)的系数矩阵和右端常数项矩阵；

假设为了确保采样点的真实值与采样点的估计值相等或相近，HASM的计算公式转换为以等式约束的最小二乘问题：

$\left\{\begin{array}{c}min\left|\right|{\mathrm{ZU}}^{n+1}-q^n|{|}_2\\ s.t.{\mathrm{CU}}^{n+1}=d\end{array}\right.---\left(18\right)$

其中C和d分别为采样点系数矩阵和采样点的值，它们的表达式分别为：C(k,(i-1)·J+j)＝1，即第k个采样点(x_i,y_j)的值为

为了求解最小二乘问题的方程组(11-17)，引入一个正的足够大的参数，λ是赋予采样点的权重，它决定于采样点对模拟曲面的贡献；对于足够大的λ，式(11-17)可被转化为无约束的最小二乘问题：

$\underset{f}{min}\left|\right|\left[\begin{array}{c}Z\\ \mathrm{\λ}C\end{array}\right]U^{n+1}-\left[\begin{array}{c}{q}^n\\ \mathrm{\λ}d\end{array}\right]|{|}_2---\left(19\right)$

也就是求解：

$\[\begin{array}{cc}{Z}^T& {\mathrm{\λC}}^T\end{array}\]\left[\begin{array}{c}Z\\ \mathrm{\λ}C\end{array}\right]U^{n+1}=\[\begin{array}{cc}{Z}^T& {\mathrm{\λC}}^T\end{array}\]\left[\begin{array}{c}{q}^n\\ \mathrm{\λ}d\end{array}\right]---\left(20\right).$

4.如权利要求1所述的地形平缓区土壤重金属镉空间分布预测方法，其特征在于，HASM模型利用采样点数据通过一系列迭代模拟步骤对某一生态曲面的空间分布进行全局模拟，进而实现对未采样点值的预测；HASM的迭代模拟步骤具体包括：

(a)利用采样点数据在计算区域内进行预插值，得到未采样点的预插值结果，即HASM模型的初始输入曲面；

(b)采用地理变量的初始输入曲面计算曲面的第一基本量、第二基本量以及HASM等式的系数；

(c)通过求解HASM等式得到接近真实空间分布格局的结果；

(d)重复上面迭代过程，直至达到满意精度。

5.如权利要求1所述的地形平缓区土壤重金属镉空间分布预测方法，其特征在于，径向基函数人工神经网络模型具有单隐层的三层前馈网络，其表达式如下：

f(x_i,k,y_j,k)＝RBFNN[F1(x_i,k,y_j,k),F2(x_i,k,y_j,k),...,F_n(x_i,k,y_j,k)] (21)

f(x_i,k,y_j,k)为点位(x_i,y_j)处土壤镉含量值，由该点位所处环境因素共同决定；F1(x_i,k,y_j,k),F2(x_i,k,y_j,k),...,Fn(x_i,k,y_j,k)为点位(x_i,k,y_j,k)处影响土壤镉含量的n个影响因素；土壤镉含量与影响因子间存在的非线性关系，采用径向基函数人工神经网络模型表达；

人工神经网络能通过自动分析多源输入与输出间的非线性映射的关系的方式，有效解决非线性预测。

6.如权利要求1所述的地形平缓区土壤重金属镉空间分布预测方法，其特征在于，对地形平缓区土壤重金属镉空间分布预测的结果还需进行评价，该评价方法为：

以RBF1_HASM方法对验证样点的模拟值与实测值间的相关系数R、平均绝对误差MAE、平均相对误差MRE和均方根误差RMSE评价各模拟方法的模拟结果精度；其中，MAE、RMSE和MRE的计算公式如下：

$MAE=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}|Z_{obs\left(i\right)}-Z_{pred\left(i\right)}|---\left(22\right)$

$RMSE=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{(Z_{obs\left(i\right)}-Z_{pred\left(i\right)})}^2}---\left(23\right)$

$MRE=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}\frac{|Z_{obs\left(i\right)}-Z_{pred\left(i\right)}|}{{\mathrm{Cd}}_{obs\left(i\right)}}\×100---\left(24\right)$

式中，n为样本数；Z_obs(i)为第i个样点的实测值mg/kg；Z_pred(i)为预测方法对第i个样点的预测值mg/kg，MAE、MRE和RMSE越小，模拟方法的误差越小；相关系数R越大，预测方法的预测结果越好。