(matlab)移动曲面拟合详细实现

作业二

目的

掌握 移动曲面法 数字高程模型内插原理及其内插子程序的设计方法，了解其它逐点高程内插方法的基本原理。

内容及原理

原理：利用已知的点数据，利用最小二乘的原理，拟合出一个可以很好的符合这些已知数据的曲面，先假设这个曲面的表达式形式（可以是一次、二次、三次等）未知其中的一些系数，通过已知点，进行间接平差求得曲面表达式。
一次表达式：
数据:

点号	X	Y	Z
1	102	110	15
2	109	113	18
3	105	115	19
4	103	103	17
5	108	105	21
6	105	108	15
7	115	104	20
8	118	108	15
9	116	113	17
10	113	118	22

求待定点X=110，Y=110点的高程。
例如二次曲面的表达式为
$$\widehat{Z_i}=A{\bar{X}}_i^2+B\bar{X}{\bar{Y}}_i+C{\bar{Y}}_i^2+D{\bar{X}}_i+E{\bar{Y}}_i+F$$
$A、B、C、D、E、F$为待定系数，那么此处为未知数，虽然称其为系数，但是在间接平差中真正的系数应该是已知数，所以将系数写在前面为：
$$\widehat{Z_i}={\bar{X}}_i^{2}A+\bar{X}{\bar{Y}}_{i}B+{\bar{Y}}_i^{2}C+{\bar{X}}_{i}D+{\bar{Y}}_{i}E+F$$
误差方程 是这样子：
$$\begin{gathered} V_i=\widehat{Z_i}-Z_i \\ V=MX-Z; \\ X=\left[\begin{array}{c}A\\ B\\ C\\ D\\ E\\ F\end{array}\right] \end{gathered}$$
注意区别$Z$与$Z_i$，前者是余下的常数，后者是真的高程,但在此公式中因为不存在其他常数项所以两者相等。之后就是套公式求得$$\begin{gathered} (\hat{A}、\hat{B}、\hat{C}、\hat{D}、\hat{E}、\hat{F})=argmin{V}^{T}V \\ X=(M^{T}PM{)}^{-1}{M}^{T}PZ \end{gathered}$$

步骤及结果

• 读入数据以待定点为圆心转换坐标

   obj=[110,110];
  know=[102 110 15;
  109 113 18;
  105 115 19;
  103 103 17;
  108 105 21;
  105 108 15;
  115 104 20;
  118 108 15;
  116 113 17;
  113 118 22];
  know(:,1:2)=know(:,1:2)-obj;

• 计算所有点与待定点的距离，计算权阵(本例中采用反距离权，即距离越远对估计点的高程影响越小)

  dis=disCalculate(know,obj);
  P=diag(1./(dis.^2));

disCalculate()为距离计算函数自定义的

function [dis] = disCalculate(know,obj)
  n=size(know,1);
  dis=zeros(n,1);
  for i=1:n
  dis(i)=sqrt((know(i,1)-obj(1))^2+(know(i,2)-obj(2))^2);
  end
 end

• 计算M、Z矩阵

[M,Z]=Mcalculator(know);

Mcalculator（此处定义了三种曲面，就有三种计算稀疏的方式可以算出来之后求$RMSE$，确定哪个最好）
一次：

function [M,Z] = Mcalculator1(know)
n=size(know,1);
M=zeros(n,3);
for i=1:n
    M(i,1)=know(i,1);
    M(i,2)=know(i,2);
    M(i,3)=1;
end
Z=know(:,3);
end

二次：

function [M,Z] = Mcalculator2(know)
n=size(know,1);
M=zeros(n,6);
for i=1:n
    M(i,1)=know(i,1)^2;
    M(i,2)=know(i,2)*know(i,1);
    M(i,3)=know(i,2)^2;
    M(i,4)=know(i,1);
    M(i,5)=know(i,2);
    M(i,6)=1;
end
Z=know(:,3);
end

三次：

function [M,Z] = Mcalculator3(know)
%三次ax^3+bx^2y+cxy^2+dy^3+ex^2+fxy+gy^2+hx+iy+j

n=size(know,1);
M=zeros(n,10);
for i=1:n
    M(i,1)=know(i,1)^3;
    M(i,2)=know(i,1)^2*know(i,2);
    M(i,3)=know(i,2)^2*know(i,1);
    M(i,4)=know(i,2)^3;
    M(i,5)=know(i,1)^2;
    M(i,7)=know(i,2)^2;
    M(i,6)=know(i,2)*know(i,1);
    M(i,8)=know(i,2);
    M(i,9)=know(i,1);
    M(i,10)=1;
end
Z=know(:,3);
end

• 计算待定系数，确定待定点的高程
  由于待定点的X、Y都为零，那么带入表达式所以高程就是待定系数F。

X=(M'*P*M)\(M'*P*Z);
Z_obj=X(6);

• 误差分析：即用估计出的曲面在每个已知点得到的$\hat{Z}$与真实$Z$差的平方和的开根号。
  以二次曲面为例,matlab代码为

%% 误差估计
f=@(x,y) X(1)*x.^2+X(2).*x.*y+X(3)*y.^2+X(4).*x+X(5).*y+X(6);
rmse=0;
for i=1:10
    rmse=rmse+(f(know(i,1),know(i,2))-know(i,3))^2;
end
rmse=sqrt(rmse);
disp(rmse)

下表是用分别用一次、二次、三次曲面的插值结果和均方根误差

方法	目标点高程值	均方根误差
一次曲面	17.9748	7.3719
二次曲面	17.9756	2.6241
三次曲面	22.1482	22.1482

可以看到一次面与三次曲面拟合精度不如二次曲面，由此可知在此例中待拟合区域地形起伏不大，但也不是没有起伏，所以使用二次曲面。
以下是曲面插值结果

👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊👊