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

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

杜玉军

(武汉大学测绘工程0408班 31610007 武汉 430079)

摘要：移动曲面拟合法是DEM格网点内插常用的一种方法，利用Matlab可以轻松实现该方法生成DEM。

关键字：移动曲面拟合法 DEM Matlab

1．概述

为了获取规则格网DEM，内插是必不可少的过程。内插的方法很多，其中移动曲面拟合法由于其方法灵活、计算简便、精度较高、占用内存较少等诸多优点而经常被使用。

2．实现原理

移动曲面拟合法是一种以待定点为中心的逐点内插法，它以每个待定点为中心，定义一个局部函数去拟合周围的数据点。其过程为：

(1) 对每个格网点，从数据点中检索出邻近的n个(至少6个)数据点。

以待定点(X,Y)为圆心，以选定长R为半径作圆，凡落入圆内的数据点都被采用。

Xpi=Xi-X；Ypi=Yi-Y

di2= Xpi2+Ypi2

di

(2) 列立误差方程式。

选择二次曲面Z=Ax2+Bxy+Cy2+Dx+Ey+F为拟合面，则数据点pi对应的误差方程式为

vi=Xpi2A+XpiYpiB+Ypi2C+XpiD+YpiE+F-Zi

n个数据点列出的误差方程可写为：

v=MX-Z X=[A B C…F]T

(3) 计算每一数据点的权。

本文选取Pi=1/di2定权。

(4) 求解待定点高程。

根据平差理论解出二次方程的系数阵：

X=(MTPM)-1MTPZ

则系数F就是待定点内插高程Zp。

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

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)

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个数据点，32×32格网，间距为0.5m；例2对peaks函数随机取了100个数据点，100×100格网，间距1m。

程序运行效果见图1~3

图3 DEM生成效果图(例1)

图4 用该方法拟合生成的peaks模型(例2)

4．结果分析

(1)关于移动曲面拟合法：

该方法基于平差理论，采用局部函数加权拟合，较为严密，但计算精度与参与平差的数据量n(n 6)有关：若n=6，即没有多余数据，误差较大，故本人认为至少要有8~10个点，这样拟合的曲面才较为光滑。

拟合结果还与权的选择有关，这要根据实际地形决定，本例采用的是距离平方倒数定权的方式。

此外，移动曲面拟合法也有一些缺点，比如计算速度较慢。

(2)关于本程序：

在搜索每个点上邻近数据点时，经验做法是以数据点距该格网点最小距离的3倍为初始搜索半径；本程序采用的是预估计范围的方法，即：先计算出单个数据点平均占用的面积S0，继而得到n个点的大致范围n×S0，将其视为正方形，以其中心到正方形顶点的距离(对角线长一半)为初始搜索半径(R=(n×S0)× /2)。

可以看到，生成的DEM与原貌拟合的较好，程序运行效果较为满意。但本程序算法上还存在一些缺陷，如没有将数据进行分块，对于每个格网点都要进行全范围的搜索，计算所有数据点到该点的距离，使计算时间非常长。本例共15个数据点、32×32的格网，大约要1秒钟；若有100个数据点、200×200的格网，则计算时间达到2.5分钟；如果数据再多，则不会在短时间内完成。故还无法处理大量数据。

参考文献：

(1) 王佩军,徐亚明.摄影测量学(测绘工程专业).武汉大学出版社..4

(2) 张祖勋,张剑清.数字摄影测量学.武汉大学出版社..1

如果觉得《matlab曲面拟合的算法 用Matlab 实现移动曲面拟合法生成DEM》对你有帮助，请点赞、收藏，并留下你的观点哦！