matlab离散点生成demWord文档下载推荐.docx
《matlab离散点生成demWord文档下载推荐.docx》由会员分享，可在线阅读，更多相关《matlab离散点生成demWord文档下载推荐.docx（9页珍藏版）》请在冰点文库上搜索。
di2=Xpi2+Ypi2
di<
R时（Xi,Yi）就被选用。
若不足个，则扩大搜索半径。
（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;
%网格长宽
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函数，内插网格点的高程
end
%****图像输出****%
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）;
函数文件1
functionDt=GetData
%读入数据
[filename,pathname]=uigetfile（'
*.txt'
'
Pickafileforread'
%打开标准对话框
fid1=fopen（strcat（pathname,filename）,'
rt'
%以只读形式打开
iffid1==-1
%若没有选择文件则警告
InputFileorPathisnotcorrect'
Warning'
warn'
break;
Dt=load（filename）;
%获取数据
fclose（fid1）;
函数文件2
functionzp=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:
N
xp（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;
nin1=n-1;
nin=nin1;
whilenin<
8
%若点数不足8个则继续搜索区外
n0=n0+（8-nin）;
%扩大搜索区
forl=1:
N-nin1
ifdout2（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;
%将本次采用的点的距离置无穷以避免重复搜索
nin=n-1;
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
输出结果
说明：
运行程序后会弹出“Pickafileforread”打开文件对话框，选择datas.txt数据文件（格式为：
点号
x坐标
y坐标
z坐标）并打开；
接着会弹出消息提示窗“请在命令窗口输入步长！
”，在CommandWindow输入格网间距回车便可以绘制DEM及等高线。
例1选用了15个数据点，32x32格网，间距为0.5m；
例2对peaks函数随机取了100个数据点，100x100格网，间距1m。
程序运行效果见图1~3
图3
DEM生成效果图（例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分钟；
如果数据再多，则不会在短时间内完成。
故还无法处理大量数据。
参考文献：
王佩军,徐亚明.摄影测量学（测绘工程专业）.武汉大学出版社.2006.4
张祖勋,张剑清.数字摄影测量学.武汉大学出版社.2005.1
点坐标生成DEM（inmatlab）
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;
21.5;
22.7];
step=[20212223242526];
%等高距1m%
figure;
plot（x,y,'
.'
markersize'
10）;
tri=delaunay（x,y）;
triplot（tri,x,y）;
holdoff;
trimesh（tri,x,y,z）;
[xi,yi]=meshgrid（0:
0.1:
11,10:
22）;
zi=griddata（x,y,z,xi,yi,'
v4'
[c,h]=contour（xi,yi,zi,step）;
clabel（c,h）;
x坐标'
y坐标'
升级会员 