1、用 Matlab 实现移动曲面拟合法生成 DEM标签： Matlab 移动曲面拟合法 DEM 分类： E【 MATLAB 】2006-12-22 10:02用 Matlab 实现移动曲面拟合法生成 DEM杜玉军（武汉大学测绘工程 0408 班 200431610007 武汉 430079）摘要：移动曲面拟合法是 DEM 格网点内插常用的一种方法，利用 Matlab 可以轻松实现该方法生成 DEM。关键字：移动曲面拟合法 DEM Matlab1 概述为了获取规则格网 DEM，内插是必不可少的过程。内插的方法很多，其中移动曲面拟合法由于其方法灵活、计算简便、精度较高、占用内存较少等诸多优点而经常被
2、使用。2 实现原理移动曲面拟合法是一种以待定点为中心的逐点内插法，它以每个待定点为中心，定义一个局部函数去拟合周围的数据点。其过程为：（1） 对每个格网点，从数据点中检索出邻近的个（至少 6 个）数据点。以待定点（）为圆心，以选定长为半径作圆，凡落入圆内的数据点都被采用。Xpi=Xi-XYpi=Yi-Ydi2= Xpi2+Ypi2diR 时（X i,Yi）就被选用。若不足个，则扩大搜索半径。（2） 列立误差方程式。选择二次曲面 Z=Ax2+Bxy+Cy2+Dx+Ey+F 为拟合面，则数据点 pi 对应的误差方程式为vi=Xpi2A+XpiYpiB+Ypi2C+XpiD+YpiE+F-Zin 个
3、数据点列出的误差方程可写为：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%移动曲面拟合法生成 DEMclear;clc;%*读入数据*%Pt=GetData; %调用 GetData 函数，读入数据ptn=num2str(Pt(:,1); %取点号ptx=Pt(:,2);
4、 %取 xpty=Pt(:,3); %取 yptz=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); %网格
5、数量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; %网格纵坐标s=l_x*l_y/length(Pt); %单点大致占用面积z(i,j)=GridZ(Pt(:,2:4),x(i),y(j),s); %调用 GridZ 函数，内插网格点的高程endend%*图像输出*%mesh(x,y,z);%hidden off;colorbar;hold on;plot3(ptx,pty,ptz,r.);text(ptx,pty,ptz+0.3,
6、ptn);title(移动曲面拟合法生成 DEM 模型);xlabel(x);ylabel(y);zlabel(z);contour(x,y,z,V);函数文件 1function 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);br
7、eak;endDt=load(filename); %获取数据fclose(fid1);函数文件 2function 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;
8、 ptout.z=0; dout2=0; %搜索区外的数据点坐标和距离P=0; %权阵for k=1:Nxp(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(
9、k);m=m+1;endendnin1=n-1;nin=nin1;while nin8 %若点数不足 8 个则继续搜索区外n0=n0+(8-nin); %扩大搜索区for l=1:N-nin1if dout2(l)s0*n0/2ptin.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;endendnin=n-1;endP=diag(1./din2); %权阵M=ptin.x.2;ptin.x.*ptin.y;pt
10、in.y.2;ptin.x;ptin.y;ones(1,n-1);X=inv(M*P*M)*M*P*ptin.z;zp=X(6); %待定点高程3.3 输出结果说明：运行程序后会弹出“Pick a file for read”打开文件对话框，选择 datas.txt 数据文件（格式为：点号 x 坐标 y 坐标 z 坐标）并打开；接着会弹出消息提示窗“请在命令窗口输入步长！”，在 Command Window 输入格网间距回车便可以绘制 DEM 及等高线。例 1 选用了 15 个数据点， 32x32 格网，间距为 0.5m；例 2 对 peaks 函数随机取了 100 个数据点，100x100
11、格网，间距 1m。程序运行效果见图 13图 3 DEM 生成效果图（例 1）图 4 用该方法拟合生成的 peaks 模型（例 2）4 结果分析（1）关于移动曲面拟合法：该方法基于平差理论，采用局部函数加权拟合，较为严密，但计算精度与参与平差的数据量（）有关：若 6，即没有多余数据，误差较大，故本人认为至少要有 810 个点，这样拟合的曲面才较为光滑。拟合结果还与权的选择有关，这要根据实际地形决定，本例采用的是距离平方倒数定权的方式。此外，移动曲面拟合法也有一些缺点，比如计算速度较慢。（2）关于本程序：在搜索每个点上邻近数据点时，经验做法是以数据点距该格网点最小距离的 3 倍为初始搜索半径；本程
12、序采用的是预估计范围的方法，即：先计算出单个数据点平均占用的面积 S0，继而得到个点的大致范围S 0，将其视为正方形，以其中心到正方形顶点的距离（对角线长一半）为初始搜索半径（R=(nS 0) /2）。可以看到，生成的 DEM 与原貌拟合的较好，程序运行效果较为满意。但本程序算法上还存在一些缺陷，如没有将数据进行分块，对于每个格网点都要进行全范围的搜索，计算所有数据点到该点的距离，使计算时间非常长。本例共 15 个数据点、3232 的格网，大约要 1 秒钟；若有 100 个数据点、200200 的格网，则计算时间达到 2.5 分钟；如果数据再多，则不会在短时间内完成。故还无法处理大量数据。参考
13、文献：（1） 王佩军, 徐亚明.摄影测量学(测绘工程专业).武汉大学出版社.2006.4（2） 张祖勋, 张剑清.数字摄影测量学.武汉大学出版社.2005.1点坐标生成 DEM（in matlab)clear;x=1; 2; 3; 4; 5; 6; 7; 8; 9; 10;y=21; 12; 13; 14; 15; 16; 17; 18; 11; 19;z=21.2; 22.8; 21.7; 22.9; 21.2; 22.5; 21.3; 22.5; 21.5; 22.7;step=20 21 22 23 24 25 26; % 等高距 1m %figure;plot(x,y,.,markersize,10);tri=delaunay(x,y);hold on;triplot(tri,x,y);hold off;figure;trimesh(tri,x,y,z);xi,yi=meshgrid(0:0.1:11,10:0.1:22);zi=griddata(x,y,z,xi,yi,v4);figure;c,h=contour(xi,yi,zi,step);clabel(c,h);xlabel(x 坐标);ylabel(y 坐标);
