1++设计地形表面积测算的原理和方法

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
我们应用最大Z容差法来选择三角形:
最大Z容差法是将重要点的选取作为一个优化问题来处理,就是利用格网点原始高程和包含该点的三角形估算的高程差来动态选取重要点。
图3 最大容差法的思想
具体的做法是这样的:
S1、连接格网DEM边界四个角点中任意对角的两个点,形成初始三角网。
S2、对每个格网点,找到包含它的三角形,内插该点在所在三角形面上的高程,求出内插高程与该点原始高程之差的绝对值(称为误差)。
S3、如果所有的格网点的误差都在最大容许的范围内则处理结束,输出TIN网;反之则进行下一步。
S4、将具有最大误差的格网点插入已存在的TIN中构成新的TIN,并返回S2。
这样我们就能构造出所有的三角形。若一个三角形的面积为 则地表总面积为:

由于原数据为 的坐标,现在改为 的坐标,则总的数据点有 个数据点。能反映出来的有效数据点总共有 个,假设每一个有效点的坐标为
z2=(p13+p14+p34)/2;
s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34));
s=s+s2;
end
end
end
s
z1=(p12+p23+p13)/2;
s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13));
s=s+s1;
end
if all([z(i,j),z(i+1,j+1),z(i,j+1)]>=12)
z2=(p13+p14+p34)/2;
s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34));
图4 有效点图
根据海伦公式:
公式里的p为半周长:

又由两点距离公式:
对于估算所给矩形区域内海拔在12米以上部分的地表面积,我们假设有四点:
若边为
假设1当三角形的三个点 全位于海拔12米以上时,通过附录2中的MATLAB程序计算得总面积 平方米。
假设2当三角形的三个点 有一个位于海拔12米以上时,通过附录3中的MATLAB程序计算总面积 平方米。
因此计算精度为:
获得地形山脚线各平面控制点的坐标数据,定出山脚线的位置。一方面确定每个高程控制点的平面位置;另一方面可按勾股定理算出每个高程控制点的两个高程测算值,取其平均值得出闭合高程 (取其整数定为 )。以一定的比例在计算纸上绘制成草图,在草图上依据设计地形变化趋势,在高程控制点和山脚线之间作出相应的 条等高线,两等高线的垂直高差为 。用透明计算纸描两幅草图。读出4个方向被等高线切割的平行线的各段(简称“割线”)长度 1,按两等高线之间的高差 计算出山的斜面长度 ,微分长方形的面积 ,用 积分得到设计地形的4个图面面积统计量 乘以相应的系数 (假如图上 代表 则 为 )后得到4个山的表面积统计量 用这4个面积统计量来做参数的区间估计。在总体方差 为未知时, 需由样本均方 估计,于是:
置信区间 。
平均值 ,
标准误差
标准误差 ,按自由度 ,查学生氏 值表,由于是为工程预算提供依据,故可取一尾 于是山的表面积真值区间为 。
附录
附录1
clc, clear
a=load('data5B.txt');
amin=min(min(a)), amax=max(max(a))
x0=0:5:4000; y0=0:5:3000;
[xx0,yy0]=meshgrid(x0,y0);
mesh(xx0,yy0,a)
v=[-49,1,12,44,75,107,171,184:80:357,357:100:614,614];
figure, c=contour(x0,y0,a,v); clabel(c)
附录2
clc, clear
a=load('data5B.txt');
s=s+s2;
end
end
end
s
附录3
clc, clear
a=load('data5B.txt');
a=a';
x0=0:5:4000; y0=0:5:3000;
pp=csape({x0,y0},a)
x=0:4000; y=0:3000;
z=fnval(pp,{x,y});
[m,n]=size(z);
1 设计地形表面积测算的原理和方法
根据附件data5B.txt中所给的数据运用附录1中的MATLAB程序做出山体实体图和等高线图:
图2 实体图
图1 等高线图
根据山的实体图雏形为基础,由于附件data5B.txt中已给出了该区域的步长为5米的网格节点对应的海拔高度值,则数据已经网格化,因此引入数字地面模型中的TIN模型:
a=a';
x0=0:5:4000; y0=0:5:3000;
pp=csape({x0,y0},a)
x=0:4000; y=0:3000;
z=fnval(pp,{x,y});
[m,n]=size(z);
s=0;
for i=1:m-1
for j=1:n-1
p1=[x(i),y(j),z(i,j)];
p2=[x(i+1),y(j),z(i+1,j)];
p14=norm(p4-p1); p34=norm(p4-p3);
ifany([z(i,j),z(i+1,j),z(i+1,j+1)]>=12)
z1=(p12+p23+p13)/2;
s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13));
s=s+s1;
end
ifany([z(i,j),z(i+1,j+1),z(i,j+1)]>=12)
s=0;
fori=1:m-1
forj=1:n-1
p1=[x(i),y(j),z(i,j)];
p2=[x(i+1),y(j),z(i+1,j)];
p3=[x(i+1),y(j+1),z(i+1,j+1)];
p4=[x(i),y(j+1),z(i,j+1)];
p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1);
p3=[x(i+1),y(j+1),z(i+1,j+1)];
p4=[x(i),y(j+1),z(i,j+1)];
p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1);
p14=norm(p4-p1); p34=norm(p4-p3);
if all([z(i,j),z(i+1,j),z(i+1,j+1)]>=1wenku.baidu.com)
相关文档
最新文档