matlab离散点生成dem

用Matlab 实现移动曲面拟合法生成DEM

标签：? HYPERLINK "/Matlab/" \t "_blank" Matlab?? HYPERLINK "/%D2%C6%B6%AF%C7%FA%C3%E6%C4%E2%BA%CF%B7%A8/" \t "_blank" 移动曲面拟合法?? HYPERLINK "/DEM/" \t "_blank" DEM?分类：? HYPERLINK "/entry/2065064/" [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

(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;? %网格纵坐标

??? 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','