matlab离散点生成dem

用Matlab实现移动曲面拟合法生成DEM标签：Matlab移动曲面拟合法DEM分类：[E]【MATLAB】2006-12-2210:02用Matlab实现移动曲面拟合法生成DEM杜玉军（武汉大学测绘工程0408班200431610007武汉430079）摘要：移动曲面拟合法是DEM格网点内插常用的一种方法，利用Matlab可以轻松实现该方法生成DEM。关键字：移动曲面拟合法DEMMatlab1．概述为了获取规则格网DEM，内插是必不可少的过程。内插的方法很多，其中移动曲面拟合法由于其方法灵活、计算简便、精度较高、占用内存较少等诸多优点而经常被使用。2．实现原理移动曲面拟合法是一种以待定点为中心的逐点内插法，它以每个待定点为中心，定义一个局部函数去拟合周围的数据点。其过程为：（1）对每个格网点，从数据点中检索出邻近的个（至少6个）数据点。以待定点（）为圆心，以选定长为半径作圆，凡落入圆内的数据点都被采用。Xpi=Xi-XYpi=Yi-Ydi2=Xpi2+Ypi2diR时（Xi,Yi）就被选用。若不足个，则扩大搜索半径。（2）列立误差方程式。选择二次曲面Z=Ax2+Bxy+Cy2+Dx+Ey+F为拟合面，则数据点pi对应的误差方程式为vi=Xpi2A+XpiYpiB+Ypi2C+XpiD+YpiE+F-Zin个数据点列出的误差方程可写为：[]T（3）计算每一数据点的权。本文选取Pi=1/di2定权。（4）求解待定点高程。根据平差理论解出二次方程的系数阵：X=(MTPM)-1MTPZ则系数就是待定点内插高程。3．实例计算3.1部分数据说明ptxptyptz数据点坐标向量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);%取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);%网格数量l_x=nx*step;l_y=ny*step;%网格长宽fori=1:(nx+1)x(i)=x0+(i-1)*step;%网格横坐标forj=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');%hiddenoff;colorbar;holdon;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);·函数文件1functionDt=GetData%读入数据[filename,pathname]=uigetfile('*.txt','Pickafileforread');%打开标准对话框fid1=fopen(strcat(pathname,filename),'rt');%以只读形式打开iffid1==-1%若没有选择文件则警告msgbox('InputFileorPathisnotcorrect','Warning','warn');break;endDt=load(filename);%获取数据fclose(fid1);·函数文件2functionzp=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;%权阵fork=1:Nxp(k)=pt(k,1)-x;yp(k)=pt(k,2)-y;d2(k)=xp(k)^2+yp(k)^2;ifd2(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;endendnin1=n-1;nin=nin1;whilenin8%若点数不足8个则继续搜索区外n0=n0+(8-nin);%扩大搜索区forl=1:N-nin1ifdout2(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;ptin.y.^2;ptin.x;ptin.y;ones(1,n-1)];X=inv(M*P*M')*M*P*ptin.z';zp=X(6);%待定点高程3.3输出结果说明：运行程序后会弹出“Pickafileforread”打开文件对话框，选择datas.txt数据文件（格式为：点号x坐标y坐标z坐标）并打开；接着会弹出消息提示窗“请在命令窗口输入步长！”，在CommandWindow输入格网间距回车便可以绘制DEM及等高线。例1选用了15个数据点，32x32格网，间距为0.5m；例2对peaks函数随机取了100个数据点，100x100格网，间距1m。程序运行效果见图1~3图3DEM生成效果图（例1）图4用该方法拟合生成的peaks模型（例2）4．结果分析（1）关于移动曲面拟合法：该方法基于平差理论，采用局部函数加权拟合，较为严密，但计算精度与参与平差的数据量（）有关：若6，即没有多余数据，误差较大，故本人认为至少要有8~10个点，这样拟合的曲面才较为光滑。拟合结果还与权的选择有关，这要根据实际地形决定，本例采用的是距离平方倒数定权的方式。此外，移动曲面拟合法也有一些缺点，比如计算速度较慢。（2）关于本程序：在搜索每个点上邻近数据点时，经验做法是以数据点距该格网点最小距离的3倍为初始搜索半径；本程序采用的是预估计范围的方法，即：先计算出单个数据点平均占用的面积S0，继而得到个点的大致范围×S0，将其视为正方形，以其中心到正方形顶点的距离（对角线长一半）为初始搜索半径（R=(n×S0)×/2）。可以看到，生成的DEM与原貌拟合的较好，程序运行效果较为满意。但本程序算法上还存在一些缺陷，如没有将数据进行分块，对于每个格网点都要进行全范围的搜索，计算所有数据点到该点的距离，使计算时间非常长。本例共15个数据点、32×32的格网，大约要1秒钟；若有100个数据点、200×200的格网，则计算时间达到2.5分钟；如果数据再多，则不会在短时间内完成。故还无法处理大量数据。参考文献：（1）王佩军,徐亚明.摄影测量学(测绘工程专业).武汉大学出版社.2006.4（2）张祖勋,张剑清.数字摄影测量学.武汉大学出版社.2005.1点坐标生成DEM（inmatlab)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=[20212223242526];%等高距1m%figure;plot(x,y,'.','markersize',10);tri=delaunay(x,y);holdon;triplot(tri,x,y);holdoff;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坐标');