用Matlab 实现移动曲面拟合法生成DEM
标签： Matlab  移动曲面拟合法  DEM 分类： [E]【 MATLAB 】2006-12-22 10:02
 
用Matlab 实现移动曲面拟合法生成DEM
杜玉军
（武汉大学测绘工程0408班 200431610007  武汉  430079）
 
摘要：移动曲面拟合法是DEM格网点内插常用的一种方法，利用Matlab可以轻松实现该方法生成DEM。
关键字：移动曲面拟合法  DEM  Matlab
 
1．  概述
为了获取规则格网DEM，内插是必不可少的过程。内插的方法很多，其中移动曲面拟合法由于其方法灵活、计算简便、精度较高、占用内存较少等诸多优点而经常被使用。
 
2．  实现原理
移动曲面拟合法是一种以待定点为中心的逐点内插法，它以每个待定点为中心，定义一个局部函数去拟合周围的数据点。其过程为：
（1）          对每个格网点，从数据点中检索出邻近的个（至少6个）数据点。
以待定点（）为圆心，以选定长为半径作圆，凡落入圆内的数据点都被采用。
Xpi=Xi-XYpi=Yi-Y
di2= Xpi2+Ypi2
di<R时（Xi,Yi）就被选用。若不足个，则扩大搜索半径。
（2）          列立误差方程式。
选择二次曲面Z=Ax2+Bxy+Cy2+Dx+Ey+F为拟合面，则数据点pi对应的误差方程式为
vi=Xpi2A+XpiYpiB+Ypi2C+XpiD+YpiE+F-Zi
n个数据点列出的误差方程可写为：
[]T
（3）          计算每一数据点的权。
本文选取Pi=1/di2定权。
（4）          求解待定点高程。
根据平差理论解出二次方程的系数阵：
X=(MTPM)-1MTPZ
则系数就是待定点内插高程。
 
3．  实例计算
3.1   部分数据说明
ptx  pty  ptz  数据点坐标向量
x(i)  y(i)  z(i,j) 第i行j列的格网点坐标值
其他数据说明见相应注释。
3.2   实现代码
·脚本文件
%DEM.m
%移动曲面拟合法生成DEM
clear;
clc;
 
%****读入数据****%
Pt=GetData;  %调用GetData函数，读入数据
ptn=num2str(Pt(:,1)); %取点号
ptx=Pt(:,2); %取x
pty=Pt(:,3); %取y
ptz=Pt(:,4); %取z
 
%*****数据预处理*****%
msgbox('请在命令窗口输入步长！');
step=input('在此输入步长：\n'); %得到网格间距
x_max=max(ptx); y_max=max(pty); %计算x,y最大值
x_min=min(ptx); y_min=min(pty); %计算x,y最小值
x0=floor(x_min/step)*step;  y0=floor(y_min/step)*step; %Grid起始点
nx=ceil((x_max-x0)/step);  ny=ceil((y_max-y0)/step);  %网格数量
l_x=nx*step;  l_y=ny*step;    %网格长宽
 
for i=1:(nx+1)
    x(i)=x0+(i-1)*step;  %网格横坐标
for j=1:(ny+1)
    y(j)=y0+(j-1)*step;  %网格纵坐标
matlab等高线间隔
    s=l_x*l_y/length(Pt);  %单点大致占用面积
    z(i,j)=GridZ(Pt(:,[2:4]),x(i),y(j),s); %调用GridZ函数，内插网格点的高程
end
end
 
%****图像输出****%
mesh(x,y,z');
%hidden off;
colorbar;
hold on;
plot3(ptx,pty,ptz,'r.');
text(ptx,pty,ptz+0.3,ptn);
title('移动曲面拟合法生成DEM模型');
xlabel('x');ylabel('y');zlabel('z');
contour(x,y,z',V);
 
·函数文件1
function Dt=GetData
%读入数据
[filename,pathname]=uigetfile('*.txt','Pick a file for read'); %打开标准对话框
fid1=fopen(strcat(pathname,filename),'rt'); %以只读形式打开
if fid1==-1  %若没有选择文件则警告
    msgbox('Input File or Path is not correct','Warning','warn');
    break;
end
Dt=load(filename); %获取数据
fclose(fid1);
 
·函数文件2
function zp=GridZ(pt,x,y,s0)
%移动曲面拟合法内插(x,y)处的高程
%pt为数据点矩阵，s0为单点平均占用面积
    N=length(pt); %点数
    n0=8;        %搜索点数
    %初始化各值
    n=1; m=1;    %计数器
    xp=0; yp=0; d2=0; %原点在(x,y)时数据点坐标和与(x,y)距离的平方
    ptin.x=0; ptin.y=0; ptin.z=0; din2=0; %落入搜索范围内的数据点坐标值和与(x,y)距离的平方
    ptout.x=0; ptout.y=0; ptout.z=0; dout2=0; %搜索区外的数据点坐标和距离
    P=0; %权阵
    for k=1:N
        xp(k)=pt(k,1)-x; yp(k)=pt(k,2)-y;
        d2(k)=xp(k)^2+yp(k)^2;
        if d2(k)<s0*n0/2  %搜索半径为n0个点占用范围(当作正方形)的对角线的一半
            ptin.x(n)=xp(k);
            ptin.y(n)=yp(k);
            ptin.z(n)=pt(k,3);
            din2(n)=d2(k);
            n=n+1;
        else  %落入搜索范围外
            ptout.x(m)=xp(k);
            ptout.y(m)=yp(k);
            ptout.z(m)=pt(k,3);
            dout2(m)=d2(k);
            m=m+1;
        end
    end
    nin1=n-1;
    nin=nin1;
    while nin<8  %若点数不足8个则继续搜索区外
        n0=n0+(8-nin);  %扩大搜索区
        for l=1:N-nin1
            if dout2(l)<s0*n0/2
                ptin.x(n)=ptout.x(l);
                ptin.y(n)=ptout.y(l);
                ptin.z(n)=ptout.z(l);
                din2(n)=dout2(l);
                dout2(l)=inf;  %将本次采用的点的距离置无穷以避免重复搜索
                n=n+1;
          end
        end
        nin=n-1;
    end
    P=diag([1./din2]); %权阵
    M=[^*ptin.^2;ptin.x;ptin.y;ones(1,n-1)];