matlab 离散点生成dem.docx

用matlab 实现移动曲面拟合法生成dem标签：matlab移动曲面拟合法dem分类：e【 matlab 】2006-12-22 10:02用matlab实现移动曲面拟合法生成dem杜玉军（武汉大学测绘工程0408班 200431610007 武汉 430079）摘要：移动曲面拟合法是dem格网点内插常用的一种方法，利用matlab可以轻松实现该方法生成dem。关键字：移动曲面拟合法 dem matlab1概述为了获取规则格网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部分数据说明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); %取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; %网格长宽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,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); break;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; 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 nin8 %若点数不足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.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输出结果说明：运行程序后会弹出“pick a file for read”打开文件对话框，选择datas.txt数据文件（格式为：点号 x坐标 y坐标 z坐标）并打开；接着会弹出消息提示窗“请在命令窗口输入步长！”，在command window 输入格网间距回车便可以绘制dem及等高线。例1选用了15个数据点，32x32格网，间距为0.5m；例2对peaks函数随机取了100个数据点，100x100格网，间距1m。程序运行效果见图13图3 dem生成效果图（例1）图4 用该方法拟合生成的peaks模型（例2）4结果分析（1）关于移动曲面拟合法：该方法基于平差理论，采用局部函数加权拟合，较为严密，但计算精度与参与平差的数据量（）有关：若6，即没有多余数据，误差较大，故本人认为至少要有810个点，这样拟合的曲面才较为光滑。拟合结果还与权的选择有关，这要根据实际地形决定，本例采用的是距离平方倒数定权的方式。此外，移动曲面拟合法也有一些缺点，比如计算速度较慢。（2）关于本程序：在搜索每个点上邻近数据点时，经验做法是以数据点距该格网点最小距离的3倍为初始搜索半径；本程序采用的是预估计范围的方法，即：先计算出单个数据点平均占用的面积s0，继而得到个点的大致范围s0，将其视为正方形，以其中心到正方形顶点的距离（对角线长一半）为初始搜索半径（r=(ns0) /2）。可以看到，生成的dem与原貌拟合的较好，程序运行效果较为满意。但本程序算法上还存在一些缺陷，如没有将数据进行分块，对于每个格网点都要进行全范围的搜索，计算所有数据点到该点的距离，使计算时间非常长。本例共15个数据点、3232的格网，大约要1秒钟；若有100个数据点、200200的格网，则计算时间达到2.5分钟；如果数据再多，则不会在短时间内完成。故还无法处理大量数据。参考文献：（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=delau