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

技术特征：

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越大，预测方法的预测结果越好。