Matlab实现移动曲面拟合法

---

0 前言

  高程拟合就是根据参考点（已知点）上的高程异常值求出其他待定点上的高程异常值，在数学上属于拟合的问题。作为一种精度较高，计算方法简单，对计算机内存要求低的移动曲面拟合法（Moving Surface Fitting,MSF）常用于由离散数据点计算未知点的高程异常值，但该法计算速度相对于多项式拟合法和多面函数法来讲可能比较慢。

1 理论

  在指定半径为$R$的范围中，以已知的中心点为圆心并和周围的已知点建立一个拟合曲面，这个曲面在中心点上的内插值即为所求值，关键的是这个有限区域会随着中心点位置的移动而变化，由以上描述可以看出，移动曲面法属于局部逼近法。
  假设某测量点$A(x,y)$的平面坐标$(x,y)$和高程异常$\xi$间有如下关系：
$$\xi =f(x,y)+\epsilon$$
  式中，$f(x,y)$为高程异常的趋势面，$\epsilon$为两者的残差。引入二次曲面的数学模型：
$$\xi =a_0+a_{1}x+a_{2}y+a_3{x}^2+a_4{y}^2+a_{5}xy$$
  设测区内高程异常为的已知点有$n(n>6)$个，那么$a_0$～$a_5$这6个待定系数可运用二次曲面拟合的方法求定。由于各已知点到中心点的距离与这些点在最小二乘计算中所作贡献的大小有关，因此，对高程异常按距离进行加权，利用权函数式：
$$P(d_i)={[(R-d_i)/d_i]}^2(i=1,2,\cdots n)$$
  式中，$d_i=\sqrt{{(x_i-x)}^2+{(y_i-y)}^2}$；$R$为搜索半径。
  设中心点$A$周围的点坐标为$(x_i,y_i)$，其在移动坐标系中的坐标为：
$$\begin{array}{rl}& x_i^A=x_i-x\\ & y_i^A=y_i-y\end{array}$$
  假设参与拟合的点数为$n$，由（1）、（2）式可列误差方程：

$${\xi }_i+v_i=a_0+a_i{x}_i^A+a_2{y}_i^A+a_3{(x_i^A)}^2+a_4{(y_i^A)}^2+a_5{x}_i^A{y}_i^A$$
  $i$分别取$1,2,\cdots ,n$；各行的权为$P(d_i)$。表示成总误差方程式即为：
$$V=BX-L$$
  计算各参与拟合的控制点的权$P(i)(i=1,2,\cdots n)$并令
$$P=\left[\begin{array}{cccc}P(1)& 0& 0& 0\\ 0& P(2)& 0& 0\\ 0& 0& \ddots & 0\\ 0& 0& 0& P(n)\end{array}\right]$$
  最后，根据最小二乘原理解带权的极小值问题，求得：

$$X={(B^{T}PB)}^{-1}{B}^{T}PL$$
  代入可得${\xi }_i$。
  由以上叙述可以看出，多项式曲面拟合模型和加权平均拟合模型是移动曲面模型的基础，移动曲面模型包含了加权平均的局部适应性强特点及多项式曲面拟合的整体性较好的特点，两种模型的结合还起到抵消部分模型误差的作用，使该模型的适用范围增大、精度提高，因此这是一种有效的GNSS高程拟合方法.

2 实现

2.1 测试数据准备

  首先准备EXCEL数据，这里说明一下，最后一列BOOL中：1代表已知高程异常值点，2代表未知高程异常值点。

点号	X (/m)	Y (/m)	Z (/m)	BOOL
8	3566375.346	499179.74	21.0344	1
23	3566854.849	498567.506	21.0076	1
10	3566324.251	498659.474	21.0172	1
5	3563826.318	499348.917	21.056	1
7	3564312.203	500321.498	21.0856	1
11	3564001.762	500035.27	21.0789	1
22	3567961.396	498691.018	20.9978	1
26	3567524.909	500219.017	21.0516	1
4	3565549.066	498813.558	21.0285	1
28	3568016.26	499235.369	21.0138	1
2	3565858.08	499248	21.0432	0
3	3564029.592	499613.378	21.0654	0
6	3566091.401	499632.434	21.0554	0
9	3564827.161	500392.773	21.0826	0
21	3566814.669	499080.199	21.0273	0
1	3564231.786	499937.723	21.0744	0
20	3567660.247	499189.334	21.017	0

2.2 代码实现部分

%% 文件打开
clear;clc;close all;
[data]=xlsread('MLS.xlsx');%文件位置需要自行更改
%% 数据预处理
KPoint=[];                  %高程异常值已知的点
UPoint=[];                  %高程异常值未知的点
for i=1:size(data,1)
    if data(i,5)==1
        KPoint=[KPoint;data(i,1:4)];
    else
        UPoint=[UPoint;data(i,1:4)];
    end
end
x=KPoint(:,2);              %已知高程异常值点x坐标
y=KPoint(:,3);              %已知高程异常值点y坐标
z=KPoint(:,4);              %已知高程异常值点z坐标
X=UPoint(:,2);              %未知高程异常值点X坐标
Y=UPoint(:,3);              %未知高程异常值点Y坐标
Zz=UPoint(:,4);             %未知高程异常值点Z真值
%% 系数计算
%Z =MSF(x,y,z,X,Y);
%MSF	移动二次曲面插值函数
%x,y	已知高程异常值点坐标x,y
%z      已知高程异常值
%X,Y	未知高程异常值点坐标X,Y
[m0, n0]=size(x);
[m1, n1]=size(Y);
for i=1:m1
    a=m0*n1*(i-1)+1;
    b=m0*n1*(i-1)+m0;
    %生成距离矩阵Dis(m0*m1*n1,n0)
    Dis(a:b,:)=sqrt((X(i)-x).^2+(Y(i)-y).^2);
    if find(Dis(a:b,:)==0)
        [m2,n2]=find(Dis(a:b,:)==0);
        Z(i,1)=z(m2,n2);
    else
        D=Dis(a:b,:);
        dx=sort(D,'descend');%'ascend' 表示升序（默认值），'descend' 表示降序
        R=dx(end-5);         %确定搜索半径选择距离最近的7个点
        [m2,n2]=find(Dis(a:b,:)<=R);
        A=ones(size(m2,1),6);
        A(:,1)=1;
        A(:,2)=x(m2)-X(i);
        A(:,3)=y(m2)-Y(i);
        A(:,4)=(x(m2)-X(i)).*(y(m2)-Y(i));
        A(:,5)=(x(m2)-X(i)).^2;
        A(:,6)=(y(m2)-Y(i)).^2;
        L=z(m2);             %已知高程异常值
        %定权1反距离定权
        Pl=diag(1./D(m2).^2);
        %定权2
        %Pl=diag((((R-D(m2)))./D(m2)).^2);
        %定权3
        %Pl=diag(exp(-(D(m2)./mean(D(m2))).^2));
        %三次样条函数定权
%         for j=1:size(m2,1)
%             sb=abs(x(m2(j))-X(i))/R;
%             if sb<0.5
%                 w(j)=2/3-4*sb^(-2)+4*sb^(-3);
%             elseif sb>1
%                 w(j)=0;
%             else
%                 w(j)=4/3-4*sb+4*sb^(-2)-4/3*sb^(-3);
%             end
%         end
%         Pl=diag(w);
    end
    
    %--计算**最小二乘法**下的曲面系数
    X0=(A'*Pl*A)\(A'*Pl*L);    %同时作为整体最小二乘法的迭代初值
    %X0=lsqminnorm(A'*Pl*A,A'*Pl*L);
    %加入岭参数避免方程病态求解
%     for k=1:10000
%         X1=(A'*Pl*A+(k*6*eye(6)))\(A'*Pl*L);    %同时作为整体最小二乘法的迭代初值
%         if abs(X1(1)-Zz(i,1))<abs(X0(1)-Zz(i,1))
%             X0=X1;
%         end
%     end
    Z(i,1) = X0(1);
end

%% 评价指标计算
V=Z-Zz;          %残差计算(单位mm)
Aver=mean(abs(V));      %误差均值
IP=sqrt(sum(V.^2)/(size(Zz,1)-1));  %内符合精度inner precision
OP=sqrt(sum(V.^2)/(size(z,1)-1));   %外符合精度outer precision

%% 结果输出
Result=[UPoint Z];          %将结果保存再UPoint矩阵中，便于对比

  第5列为由平面位置拟合得到的高程异常值，第4列为该点处高程异常值验证值。

2.3 进一步获取结果

  根据2.2中的计算结果，我们可以进一步输出得到自己需要的图，这里只展示几种简单的，想要生成更多图形可参靠MATLAB官方支持.

%绘图输出
SResult=sortrows(Result,2); %按照点号排序

ff = figure('Name','结果输出','NumberTitle','off');
set(ff, 'Position',[351,96,1346,832],'color','white')

subplot(221); 
scatter(X,Y,20,[0 0 1],'+'); 
hold on%将拟合点和控制点同时输出在同一张图上
scatter(x,y,20,[1 0 0],'o');
text(x,y,num2str(KPoint(:,1))); %输出控制点点号
text(X,Y,num2str(UPoint(:,1))); %输出拟合点点号
legend('未知点','已知点');
title('坐标点平面位置图');
xlabel('X');
ylabel('Y');

subplot(222); 
plot(SResult(:,2),SResult(:,4),'-*',SResult(:,2),SResult(:,5),'-o');
%横轴为X，与坐标点平面位置图相对应
text(SResult(:,2),SResult(:,4),num2str(SResult(:,1))); %标注点号
legend('拟合前高程异常值（m）','拟合后高程异常值（m）');
title('拟合前后高程异常对比图');
xlabel('X');
ylabel('高程异常值(m)');

subplot(224); 
bar(SResult(:,2),(SResult(:,5)-SResult(:,4))*1000,'group','b');
text(SResult(:,2),(SResult(:,5)-SResult(:,4))*1000,num2str(SResult(:,1))); %标注点号
%横轴为X，与坐标点平面位置图相对应
title('残差直方图');
xlabel('X');
ylabel('残差值(mm)');

 
               赶紧点赞、收藏起来吧！不然下次就找不到了💕