本文已参与「新人创作礼」活动，一起开启创作之路。

本文介绍根据MATLAB完成大局多项式插值法与逆间隔加权法的空间插值的办法，并对不同插值办法成果加以比照剖析。

为方便咱们理解，本文将讲解与代码部分独立开来；首先将详细的操作流程与办法细心论述，其次将本文所用到的悉数代码完整附于本文末尾。其间，因为本文所用的数据并不是我的，因而惋惜不能将数据一并提供给咱们；但是根据本篇博客的思想与对代码的详细解释，咱们用自己手头的数据，能够将相关操作与剖析进程加以完整重现。

1 布景常识

空间数据的获取是进行空间剖析的根底与来源，咱们总是期望能够获取研究区域更多、更全面的准确空间特点数据信息。但在实际工作中，因为成本、资源等条件约束，不或许对悉数不知道区域加以测量，而更多只能得到有限数量的采样点的数据。因而，能够考虑选取合适的空间采样点，使用必定数学模型，根据已知采样点数据对研究区域一切方位的不知道特点信息加以猜测。

空间插值（Spatial Interpolation）便是一种将离散点测量数据转换为连续数据曲面的常用办法，包含内插（Interpolation）和外推（Extrapolation）两种使用办法。空间插值根据地理学第必定律，即一般地，间隔越近的地物具有越高的相关性。

在办法层面，空间插值一般能够分为确定性插值办法（Deterministic Interpolation）与地质统计学办法（Geostatistics）。前者根据信息点之间类似程度或整个曲面的平滑程度创立拟合曲面，后者则根据信息点归纳统计学规律，对其空间自相关性定量化，然后创立插值面。另一方面，根据插值核算时纳入考虑的采样点散布规模，又分为全体插值法与部分插值法。前者使用整个实测采样点数据集对全区进行拟合，如大局多项式插值法（Global Polynomial Interpolation）；后者则仅仅用接近某一区域内的采样点数据猜测不知道点的数据，如反（逆）间隔加权法（Inverse Distance Weighted，IDW）。此外，根据插值成果曲面中采样点猜测值与实测值的关系，又可分为准确性插值与不准确插值。

本文借助MATLAB软件自主编程，别离使用大局多项式插值法与逆间隔加权法，对湖北省荆门市沙洋县土壤pH值、有机质含量等两种特点数据进行空间插值核算，并比照对应插值办法的拟合作用。

大局多项式插值法以悉数采样点覆盖区域为根底，经过最小二乘法等手段拟合出一个最合适的平面或曲面，使得各个采样点较为均匀地散布于这一平面或曲面的邻近，且悉数高出该面的点距之和与悉数低于该面的点距之和的绝对值应当近似。大局多项式插值法不适合描绘特点数据空间散布动摇较大区域的特征。

逆间隔加权法则使用待插值点周围必定规模的已知点数据，对中心待插点加以数据插值；某点间隔待插点越近，其对这一点的影响越大，即对应系数越大。而权重的核算往往依赖于反间隔（间隔倒数）的p次方。一般地，取p等于2。根据办法原理可知，逆间隔加权法往往会导致空间散布的多点中心现象。

2 实际操作部分

2.1 空间数据读取

本文以MATLAB核算变异函数并制造经验半方差图中的数据为范例，其包含湖北省荆门市沙洋县658个土壤采样点的空间方位（X与Y，单位为米）、pH值、有机质含量与全氮含量。这些数据相同存储于data.xls文件中，而后期操作多根据MATLAB软件加以完成。因而，首先需根据详细核算需求将源数据挑选性地导入MATLAB软件中。

2.2 反常数据除掉

同上述博客所述，咱们得到的原始采样点数据在采样记载、实验室测验等进程中，或许具有必定差错，然后出现单个反常值；而反常值的存在会对后期空间插值作用（尤其是部分插值办法）产生较大影响。因而，选用均匀值加标准差法对反常数据加以挑选、除掉。这一办法的详细意义请见MATLAB核算变异函数并制造经验半方差图。

上述博客中，别离使用“2S”与“3S”办法加以处理，发现“2S”办法处理作用相对后者较好。故本文直接挑选使用“2S”办法处理成果持续进行。

得到反常值后，将其从原有658个采样点中除掉；剩余的采样点数据持续后续操作。

2.3 验证集挑选

针对不同插值办法所得成果的精度查验，本文采纳设置练习数据（Training Data）与测验数据（Test Data）并比照的办法进行。将经过上述反常除掉后的数据随机分为两个部分，其间80%作为练习数据，剩余20%作为测验数据（即验证集）。需求注意的是，测验数据仅可用于检测插值模型的构建作用，而不可参与模型构建进程，否则或许导致成果过度拟合（即模型对本次所给数据具有较好插值、拟合作用，但对外围无数据区域的猜测作用较差）。因为验证集数据为随机选取，因而后期可屡次履行程序，以期取得平稳的精度衡量目标。

2.4 最小二乘法求解

正如大局多项式插值法称号所示，这一插值办法的原理实际上是对一个描绘平面或曲面的多项式各个系数加以求解；而这一求解进程往往采用最小二乘法完成。

考虑到多项式的复杂度，本文将多项式阶数限定为二阶与三阶，并对其插值作用加以比照。二阶多项式办法如下：

三阶多项式办法如下：

上述两公式中，W(x,y)为待插点(x,y)的插值，x与y为坐标，为常数，其它字母代表各对应项系数。

若使用AcrMap等软件进行趋势面插值，则可根据实际情况恰当提升阶数。

2.5 逆间隔加权法求解

逆间隔加权法经过各已知点自身实际数值及其关于待插点的权重完成插值，公式如下：

其间，_i为第i个已知点对待插点(x_0,y_0 )的权重，d_i0为二者间间隔，n为已知点个数，p为幂参数，W(x_0,y_0 )为待插点(x_0,y_0 )的插值，Z(x_i,y_i )为第i个已知点对应特点数值。

本文中，取初始p=2，并根据插值作用恰当对其加以调整。

2.6 插值精度查验

如前所述，本文经过随机选定的测验数据对插值成果的精度进行比较与剖析。结合实际空间数据的数值特点，终究所选用以衡量插值精度的目标包含均匀差错（Mean Error，ME）、均匀绝对差错（Mean Absolute Error，MAE）、均方根差错（Root Mean Square Error，RMSE）与相关系数（Correlation Coefficient）。

其间，均匀差错能够获悉插值成果与实测点观测值的巨细关系；均匀绝对差错表明空间插值与实测点观测值之间绝对差错的均匀值，能够更好反映插值成果差错的实际情况；均方根差错表明插值成果与实测点观测值之间差异（即残差）的样本标准差；相关系数（这里选用皮尔逊相关系数）用以评价空间插值与实测点观测值之间的线性相关度。

2.7 数据导出与专题地图制造

经过上述进程，得到对应特点数据的插值成果。首先在MATLAB制造插值成果三维图，随后将插值成果数据变量保存为ASCII数据格局文件，并导入AcrMap软件制造专题地图。

其间，在数据由MATLAB导出为.txt格局后，需求在其开头部分增加以下信息：

ncols	1275
nrows	1209
xllcorner	600800
yllcorner	3364600
cellsize	50
NODATA_value	-9999

其间，ncols与nrows对应数值别离代表插值成果列数与行数，xllcorner与yllcorner对应数值别离代表插值成果x、y坐标的最小值（即起始值），cellsize对应数值代表像元巨细。

3 成果出现与剖析

经过均匀差错、均匀绝对差错、均方根差错、相关系数等四个精度衡量目标，以及对应趋势面公式与三维插值成果图，将不同空间插值办法所得成果比照、剖析如下，并制造专题地图。

3.1 大局多项式插值法二阶与三阶插值比照

使用最小二乘法，别离对二阶多项式与三阶多项式进行系数求解，并得到pH值、有机质含量等两种空间特点数据的大局多项式插值成果与各精度目标。针对不同阶数多项式，别离履行六次插值操作，并求出均匀值。得到成果别离如表1至表6所示。

由表1至表3可知，针对pH值的大局多项式插值法，二阶、三阶多项式所得插值成果的均匀差错均为负数，即两种办法均趋向于取得较之观测值高的插值成果；而后者所得均匀差错的数值较小于前者（即后者这一目标绝对值较小）。二阶、三阶多项式插值成果对应均匀绝对差错、均方根差错相差不大，但后者上述两种目标数值相同小于前者。三阶多项式相关系数相同略大于二阶多项式。综上所述，面向pH值的大局多项式插值法，其运用三阶多项式的插值作用较优于运用二阶多项式的插值作用，且这一成果在均匀差错、均匀绝对差错、均方根差错与相关系数等四个精度衡量目标中均有所表现。

由表4至表6可知，针对有机质含量的大局多项式插值法，二阶、三阶多项式所得插值成果的均匀差错均为正数，即两种办法均趋向于取得较之观测值低的插值成果；而后者所得均匀差错的数值较大于前者。二阶、三阶多项式插值成果对应均匀绝对差错相差不大，但后者上述目标数值相同略大于前者。二阶多项式的均方根差错较之三阶多项式低，且二者这一目标有着较为显着的间隔。三阶多项式相关系数相同略小于二阶多项式。综上所述，面向有机质含量的大局多项式插值法，其运用二阶多项式的插值作用与运用三阶多项式的插值作用全体区分度不大，二阶多项式插值成果在均匀差错、均匀绝对差错与相关系数等三个精度衡量目标中略优于三阶多项式，而三阶多项式插值成果则在均方根差错这一目标中表现出色。考虑到不同测验数据的选取具有随机性，因而以为上述较为接近的成果并不能特别表现出二者中更优的插值办法挑选。

归纳表1至表6可知，对于大局多项式插值法，其多项式阶数的挑选并不必定是越高越好，而是需求根据实际数据情况加以尝试、挑选。较高的阶数除了会增大核算量、延长程序履行时间外，其插值成果精度亦有或许并无增加，乃至在某些目标中出现下降趋势。

3.2 大局多项式插值法函数及其三维成果图

综上所述，别离使用二阶多项式与三阶多项式获取大局多项式插值法对应趋势面函数。

针对pH值的二阶、三阶大局多项式插值趋势面函数如下：

其间，上述二阶与三阶趋势面函数别离对应各精度衡量目标情况如表7所示。

上述二阶与三阶趋势面函数别离对应三维插值成果图如下。

针对有机质含量的二阶、三阶大局多项式插值趋势面函数如下：

其间，上述二阶与三阶趋势面函数别离对应各精度衡量目标情况如表8所示。

上述二阶与三阶趋势面函数别离对应三维插值成果图如下。

3.3 大局多项式插值法专题地图制造

经过本文前述部分的相关办法，将MATLAB插值数据成果文件导入ArcMap，经过取舍后制造湖北省荆门市沙洋县土壤pH值、有机质含量大局多项式插值专题地图。其间，专题地图的详细制造办法这里就不再赘述了，咱们如果有需求能够搜索对应的教学内容。

得到专题地图如下所示。

由上述四幅专题地图可知，三阶多项式插值所得成果较之二阶多项式遍及具有更广的数值规模，如二阶pH值插值地图中最大成果为8.99，而其三阶插值则可到达10.34；相同，二阶有机质含量插值地图中最小成果为8.91，而其三阶插值则可低至5.53。针对这一现象，个人以为或许是因为三阶多项式插值曲面曲折程度较之二阶多项式插值更大，平滑效应较之二阶多项式不显着，因而其所能到达的数值规模亦较大。

另一方面，由专题地图这一视点观之，能够发现pH值与有机质含量散布具有必定空间关系：pH值较高区域，有机质含量往往较低（首要集中于沙洋县东部地区）；反之，有机质含量则往往较高（首要集中于沙洋县中、西部地区）。而这一成果也契合“土壤有机质与pH值之间具有较高负相关关系”的定论。

3.4 逆间隔加权法插值成果及其三维成果图

根据本文前述办法，取初始p=2，并根据插值作用恰当调整，屡次重复履行逆间隔加权法，得到pH值、有机质含量等两种空间特点数据的插值成果与各精度目标。得到成果如表9所示。

由表9可知，面向pH值与有机质含量的逆间隔加权法插值作用在上述四个精度衡量目标中表现并不是特别抱负，较之前述大局多项式插值法（包含二阶、三阶）成果差错略大。尤其是有机质含量逆间隔加权法成果的均方根差错，其均匀数值已达5.10左右，说明各次IDW办法的有机质含量插值成果与实测点观测值之间差异（即残差）的样本标准差较大。此外，逆间隔加权法有机质含量插值成果的均匀相关系数未高于0.50。而与有机质含量相比，pH值的插值成果各项精度衡量目标全体稍优；由此能够看出，pH值的逆间隔加权法插值作用全体精度较好于有机质含量；而这或许与两种特点数据各自内部的空间相关性、数值取值、散布特征等有关。由全体精度衡量目标来看，编程完成的逆间隔加权法作用较之大局多项式插值法略差。

另一方面，在屡次履行逆间隔加权法插值进程中发现，每次履行作用精度差异改变较大；尤其是有机质含量，有时上述各项精度目标能够到达一个很好的水平，而有时则较差，如均方根差错乃至有时可达8左右。由此抑或能够看出，逆间隔加权法作为一种部分插值办法，其履行作用所遭到练习数据与测验数据的选取情况影响较大。

pH值与有机质含量的逆间隔加权法插值三维成果别离如下两幅图所示。

3.5 逆间隔加权法专题地图制造

经过本文前述部分的相关办法，将MATLAB插值数据成果文件导入ArcMap，经过取舍后制造湖北省荆门市沙洋县土壤pH值、有机质含量逆间隔加权法插值专题地图。

由上述两幅专题地图可知，逆间隔加权法插值所得成果较之大局多项式插值法，生成的表面崎岖改变数量更多、程度更大，而崎岖所影响的规模则较小；且如前所述，逆间隔加权法得到的插值成果具有较多小规模的中心散布区域，即插值专题图中可见包含若干气泡状小点散布的特点数据特征。此外，pH值与有机质含量的散布特征及其二者间的空间相互关系依然同大局多项式插值法成果，即二者数据巨细间出现出相反状况，沙洋县中、西部地区pH值相对较低而有机质含量较高，东部地区则pH值相对较高而有机质含量较低。

一起，正如本文榜首部分所述，因为逆间隔加权法是一种部分插值法，每一待插点的插值成果均很大程度上遭到其接近点数值的影响；因而上述空间散布特征亦仅仅其成果的全体趋势，其间也会有部分特例。例如，结合上述两幅专题地图能够看到，逆间隔加权法所得pH值插值成果在沙洋县西部地区亦有部分部分极大值点，而这些极大值点的数值乃至与东部地区持平；一起，其有机质含量插值成果在沙洋县中、西部地区亦存在零散散布的部分极小值点。

4 完整代码

4.1 大局多项式插值法MATLAB代码

%% 文件信息读入
clc;clear;
info=xlsread('data.xls');
opoX=info(:,1);
opoY=info(:,2);
oPH=info(:,3);
% oOM=info(:,4);
% oTN=info(:,5);
%% 2S法核算反常值
mPH=mean(oPH);
sPH=std(oPH);
num2=find(oPH>(mPH+2*sPH)|oPH<(mPH-2*sPH));
%% 反常值除掉
PH=oPH;
for i=1:length(num2)
    n=num2(i,1);
    PH(n,:)=[0];
end
PH(all(PH==0,2),:)=[];
poX=opoX;
for i=1:length(num2)
    n=num2(i,1);
    poX(n,:)=[0];
end
poX(all(poX==0,2),:)=[];
poY=opoY;
for i=1:length(num2)
    n=num2(i,1);
    poY(n,:)=[0];
end
poY(all(poY==0,2),:)=[];
%% 验证集挑选
very=[randperm(length(PH),floor(length(PH)*0.2))]';
%% 验证集除掉
cPH=PH;
vPH=zeros(length(very),1);
vpoX=vPH;
vpoY=vPH;
for i=1:length(very)
    m=very(i,1);
    vPH(i,:)=cPH(m,:);
    cPH(m,:)=[0];
end
cPH(all(cPH==0,2),:)=[];
cpoX=poX;
for i=1:length(very)
    m=very(i,1);
    vpoX(i,:)=cpoX(m,:);
    cpoX(m,:)=[0];
end
cpoX(all(cpoX==0,2),:)=[];
cpoY=poY;
for i=1:length(very)
    m=very(i,1);
    vpoY(i,:)=cpoY(m,:);
    cpoY(m,:)=[0];
end
cpoY(all(cpoY==0,2),:)=[];
%% 最小二乘法求解预处理
inva2=[ones(size(cpoX)),cpoX.^2,cpoY.^2,cpoX.*cpoY,cpoX,cpoY];
inva3=[ones(size(cpoX)),cpoX.^3,cpoY.^3,(cpoX.^2).*cpoY,cpoX.*(cpoY.^2),cpoX.^2,cpoY.^2,cpoX.*cpoY,cpoX,cpoY];
%% 最小二乘法求解
[coef2,bint2,r2,rint2,stats2]=regress(cPH,inva2);
[coef3,bint3,r3,rint3,stats3]=regress(cPH,inva3);
%% 趋势面法作用图制造预备
step=50;
npoX=600800:step:664500;
npoY=3364600:step:3425000;
[mnpX,mnpY]=meshgrid(npoX,npoY);
%% 趋势面法插值
pPH2=coef2(1,:)+coef2(2,:)*mnpX.^2+coef2(3,:)*mnpY.^2+coef2(4,:)*mnpX.*mnpY+coef2(5,:)*mnpX+coef2(6,:)*mnpY;
pPH3=coef3(1,:)+coef3(2,:)*mnpX.^3+coef3(3,:)*mnpY.^3+coef3(4,:)*(mnpX.^2).*mnpY+coef3(5,:)*mnpX.*(mnpY.^2)+coef3(6,:)*mnpX.^2+coef3(7,:)*mnpY.^2+coef3(8,:)*mnpX.*mnpY+coef3(9,:)*mnpX+coef3(10,:)*mnpY;
%% 趋势面法作用图制造
scatter3(cpoX,cpoY,cPH);
hold on;
mesh(mnpX,mnpY,pPH2);
title('Global Polynomial Interpolation Results of Quadratic of Organic Matter');
figure();
scatter3(cpoX,cpoY,cPH);
hold on;
mesh(mnpX,mnpY,pPH3);
title('Global Polynomial Interpolation Results of Cubic of Organic Matter');
%% 趋势面法精度比照
vpPH2=coef2(1,:)+coef2(2,:)*vpoX.^2+coef2(3,:)*vpoY.^2+coef2(4,:)*vpoX.*vpoY+coef2(5,:)*vpoX+coef2(6,:)*vpoY;
MEERan2=mean(vPH-vpPH2);
MEERan21=mean(abs(vpPH2-vPH));
RMSEan2=sqrt(sum(vpPH2-vPH).^2/length(vPH));
COCOan2=corrcoef(vpPH2,vPH);
vpPH3=coef3(1,:)+coef3(2,:)*vpoX.^3+coef3(3,:)*vpoY.^3+coef3(4,:)*(vpoX.^2).*vpoY+coef3(5,:)*vpoX.*(vpoY.^2)+coef3(6,:)*vpoX.^2+coef3(7,:)*vpoY.^2+coef3(8,:)*vpoX.*vpoY+coef3(9,:)*vpoX+coef3(10,:)*vpoY;
MEERan3=mean(vPH-vpPH3);
MEERan31=mean(abs(vpPH3-vPH));
RMSEan3=sqrt(sum(vpPH3-vPH).^2/length(vPH));
COCOan3=corrcoef(vpPH3,vPH);
%% 趋势面法导出ASCII
save 3.txt pPH2 -ASCII;
save 4.txt pPH3 -ASCII;

4.2 逆间隔加权法MATLAB代码

%% 文件信息读入
clc;clear;
info=xlsread('data.xls');
opoX=info(:,1);
opoY=info(:,2);
oPH=info(:,3);
power=2;
% oOM=info(:,4);
% oTN=info(:,5);
%% 2S法核算反常值
mPH=mean(oPH);
sPH=std(oPH);
num2=find(oPH>(mPH+2*sPH)|oPH<(mPH-2*sPH));
%% 反常值除掉
PH=oPH;
for i=1:length(num2)
    n=num2(i,1);
    PH(n,:)=[0];
end
PH(all(PH==0,2),:)=[];
poX=opoX;
for i=1:length(num2)
    n=num2(i,1);
    poX(n,:)=[0];
end
poX(all(poX==0,2),:)=[];
poY=opoY;
for i=1:length(num2)
    n=num2(i,1);
    poY(n,:)=[0];
end
poY(all(poY==0,2),:)=[];
%% 验证集挑选
very=[randperm(length(PH),floor(length(PH)*0.2))]';
%% 验证集除掉
cPH=PH;
vPH=zeros(length(very),1);
vpoX=vPH;
vpoY=vPH;
for i=1:length(very)
    m=very(i,1);
    vPH(i,:)=cPH(m,:);
    cPH(m,:)=[0];
end
cPH(all(cPH==0,2),:)=[];
cpoX=poX;
for i=1:length(very)
    m=very(i,1);
    vpoX(i,:)=cpoX(m,:);
    cpoX(m,:)=[0];
end
cpoX(all(cpoX==0,2),:)=[];
cpoY=poY;
for i=1:length(very)
    m=very(i,1);
    vpoY(i,:)=cpoY(m,:);
    cpoY(m,:)=[0];
end
cpoY(all(cpoY==0,2),:)=[];
%% IDW作用图制造预备
step=50;
npoX=600800:step:664500;
npoY=3364600:step:3425000;
[mnpX,mnpY]=meshgrid(npoX,npoY);
%% IDW分母核算
temdeno=zeros(1,length(cPH));
deno=zeros(length(npoY),length(npoX));
for p=1:length(npoY)
    for q=1:length(npoX)
        for i=1:length(cPH)
            temdeno(1,i)=(sqrt((npoY(1,p)-cpoY(i,1))^2+(npoX(1,q)-cpoX(i,1))^2))^(-power);
        end
        deno(p,q)=sum(temdeno(:));
    end
end
%% IDW求解
temPH=zeros(1,length(cPH));
pPH4=zeros(length(npoY),length(npoX));
for p=1:length(npoY)
    for q=1:length(npoX)
        for i=1:length(cPH)
            temPH(1,i)=cPH(i,1).*((sqrt((npoY(1,p)-cpoY(i,1))^2+(npoX(1,q)-cpoX(i,1))^2))^(-power))./deno(p,q);
        end
        pPH4(p,q)=sum(temPH(:));
    end
end
%% IDW作用图制造
scatter3(cpoX,cpoY,cPH);
hold on;
mesh(mnpX,mnpY,pPH4);
title('IDW Results of Organic Matter');
%% IDW验证1
temdeno4=zeros(1,length(cPH));
deno4=zeros(length(vpoY),length(vpoX));
for p=1:length(vpoY)
    for q=1:length(vpoX)
        for i=1:length(cPH)
            temdeno4(1,i)=(sqrt((vpoY(p,1)-cpoY(i,1))^2+(vpoX(q,1)-cpoX(i,1))^2))^(-power);
        end
        deno4(p,q)=sum(temdeno4(:));
    end
end
%% IDW验证2
temPH4=zeros(1,length(cPH));
vpPH4=zeros(length(vpoY),1);
for p=1:length(vpoY)
    for i=1:length(cPH)
        temPH4(1,i)=cPH(i,1).*((sqrt((vpoY(p,1)-cpoY(i,1))^2+(vpoX(q,1)-cpoX(i,1))^2))^(-power))./deno4(p,q);
    end
    vpPH4(p,1)=sum(temPH4(:));
end
%% 精度核算
MEERan4=mean(vPH-vpPH4);
MEERan41=mean(abs(vpPH4-vPH));
RMSEan4=sqrt(sum(vpPH4-vPH).^2/length(vPH));
COCOan4=corrcoef(vpPH4,(vPH)');
%% 文件转存
save 5ph.txt pPH4 -ASCII;

声明：本站所有文章，如无特殊说明或标注，均为本站原创发布。任何个人或组织，在未征得本站同意时，禁止复制、盗用、采集、发布本站内容到任何网站、书籍等各类媒体平台。如若本站内容侵犯了原著者的合法权益，可联系我们进行处理。

MATLAB实现全局多项式插值（趋势面法）与逆距离加权IDW插值