第五章水准网程序设计

《第五章水准网程序设计》由会员分享，可在线阅读，更多相关《第五章水准网程序设计（15页珍藏版）》请在装配图网上搜索。

1、第五章第五章 水准网程序设计水准网程序设计5.15.1 概概述述水准网是为了确定地面点的高程而布设的控制网，网中的观测值是高程控制点之间的高差。为了统一全国的高程系统，我国采用黄海平均海水面作为全国高程系统的基准面，在该面上的任一点，其高程为零。水准网中的任一点高程以及点与点之间的高差应属于正常高系统。但是，水准测量中水准仪逐站测量并累加得到的原始高差并非正常高差，对于高精度的水准网，应在原始数据基础上添加尺长改正、正常水准网面不平行、球气差改正等系统改正，才能得到正常高差，受篇幅所限，本章不讨论观测值的归算问题，本章提到的观测值均假定是已经归算后的正常高差观测值。由于观测值存在误差，实际工作

2、中需要若干已知点作为平差的基准控制点。水准网平差的目的就是求解各观测值的最佳估值及评定未知量的精度。本章开始，我们将开始接触到测量数据处理程序的设计与开发。测量数据因其数据量大、繁杂等特殊性，往往需要频繁对矩阵进行处理计算。然而在MATLAB 软件中，矩阵的计算变得格外简单，这为广大测绘工作者提供了实用性强、简单方便的数据处理平台。平差程序就是将平差计算过程程序化，综合考虑，选择参数平差模型作为水准网平差的主要模型。测量程序设计一般包含程序功能设计、平差模型选择、算法选择等内容。在本章的结构的组织上，先对水准原理进行简要分析，在此基础上按不同功能模块设计函数，最后对程序进行分析验证。5.25.

3、2 水准路线处理程序设计水准路线处理程序设计5.2.1 水准网平差函数类设计1.类设计Main%函数主体程序calculateHO（）%计算近似高程ca_ATPA()%计算法方程系数 A，权阵 GP，常数项 GLca_V()%计算高差改正数Printlevledata()%输出结果2.函数说明具体各个成员变量的含义在后面程序中会有注释，在此对主函数main 做简要强调。主函数 main 中包括数据的读取、检查数据格式、提取相关观测值、各类成员函数和最小二平差。是本程序运行的主要M 文件。5.2.2 原始数据文件格式设计水准网平差所需的数据需要从txt 文本中提取，本节主要介绍原始数据文件的内容

4、与格式，需要说明的是，原始数据文件的设计并没有统一的、严格的标准，在方便程序设计的基础下，程序设计者可以自由设计数据文件的格式。我们将网的数据分为三类：网的概况信息、已知数据、观测数据。网的概况信息包括总点数、已知点点数、观测者总数、验前单位权中误差。已知数据包括已知点名、已知高程值。观测数据包括高差起点点名、高差终点点名、观测值高差值、线路长度。下面结合实例说明原始文件的具体数据格式图 5-2-1 为典型的水准网表 5-2-1 已知点数据点名高程/m A5.160 B6.016图5-1 水准网示意图表 5-2-2 观测高差No.1234起点AABB终点P1P2P1P2h/m1.3592.00

5、90.3631.012S/km1.11.72.32.7No.567起点P1P1P3终点P2P3Bh/m0.6570.238S/km2.41.4-0.5952.6利用以上数据进行水准网平差，数据文件的内容如下：7 5 2 0.001A1 B1 P1 P2 P3A1 5.016B1 6.016A1 P1 1.359 1.1A1 P2 2.009 1.7B1 P1 0.363 2.3B1 P2 1.012 2.7P1 P2 0.657 2.4P1 P3 0.238 1.4P3 B1 -0.595 2.6格式说明：（1）第一行为网的概况信息：观测值总数、总点数、已知点总数、验前单位权中误差。（2）第二

6、行是所有点点名 依次书写。（3）第三行是对应已知点点名、高程（单位为 m），依次往下写、列出所有一直点名和高程。（4）接下来是观测高差：高差起始点名、高差终点点名、高差观测子和高差线路长度。网中点名需要字母和数字结合，否则点名读不出来，点名之间空格；实际平差时，水准网规模和数据可能会不同，但是只要按照上面的格式和顺序输入相关数据即可用本章的程序。5.2.2 数据的存储1 点数据的存储平差程序用到的数组中，点名、高程值、高程改正数等数组与高程点名一一对应，这点数据称为点数据。点数据数组的总长度为总点数。高程值用数组 Height 保存，高差改正数用数据 dX 存放，近似值存放在数组 caheig

7、ht中。点名用数组 pname 数组存放。pname 数组不仅存放数组名，同时也可看做点名地址数组。点数据数组存放顺序必须保持严格的一致，例如某点点名地址存在pname数组中的第3单元，即 pname（3）中，该点的高程值也就存放在Height（3）中，高程改正数也存放在dX（3）中。当点数据存放一致时，可以更快捷的进行点数据存储，更利于程序设计。2 观测数据的存储观测数据包括高差观测值、高差起始点的点名、高差的权、高差的改正数，它们与观测值一一对应，数组长度都等于总点数。高差起点点名和终点点名在程序中的作用是缺点观测值的起始点号和终点点号，有了点号才能访问点数据数组，进行相应的计算。利用函数

8、strcmp 将点名和点号一一对应，通过点号可以得到点名，当然也能通过点名得到点号，高差的起始点点号和终点点号分别用数组startp，endp 来保存。个观测数据数组的存放顺序也必须是一致，即L（i）、P(i)、V(i)、startp（i）、endp（i）均对应于同一个观测量。5.2.3 近似高程计算1 近似高程值计算思路第一步，设置未知点标志。各点的高程值用数组Height 保存，近似值高程计算之前，由于正常的高程值不可能小于-9999.9，将 Height 数组中未知点高程赋值为-9999.9，根据高程数组中的值是否小于-9999.0 判断某点是否需要计算近似高程，当该点近似高程计算出来之

9、后就会被取代。第二步，计算高程值。近似高程计算的基本思路是，遍历观测值，找到每一个观测值起点号、终点号，根据点号从 Height 数组中提取起点和终点的高程值,当高差的一端是高程已知点，另一端是未知点时，就由已知点高程和观测高差算出未知点高程，并将计算结果赋值给高程数组，如果两端都是未知点，就跳到下一循环。在遍历观测值时，排在前面的观测值总是最先用于高程近似值计算，当观测值的排序不同时，近似高程计算用到的观测值也不同，计算结果也不同。有时候一次遍历观测值可能还有部分点的高程值仍未计算出来，还需要再次遍历观测值，直至算出所有高程值为止。2 函数文件 calculateHO（）代码function

10、 y=calculateHO(a,b,c,d,e,f,g,h)%计算近似高程%caheight=calculateHO(m_Knumber,m_Pnumber,m_Lnumber,height,startP,endP,pname,L);jj=0;for i=(a+1):b d(1,i)=-9999.9;%未知点的初始高程为-9999.9endfor i=1:(b-a)%循环的次数不超过（m_Pnumber-m_Knumber）for j=1:c%按照观测组进行循环 K1=e(j);%起点点号 K2=f(j);%终点点号 if d(1,K1)-9999.0&d(1,K2)-9999.0%起点未知

11、，终点已知 d(1,K2)=d(1,K1)+h(1,j);jj=jj+1;elseif d(1,K1)-9999.0%起点已知，终点未知 d(1,K1)=d(1,K2)-h(1,j);jj=jj+1;end end y=d;if jj=(b-a)break end if i(b-a)%生成产生错误的 txt fid=fopen(resultfid.txt,a+);fprintf(fid,nn 下列点无法计算出概略高程:n);for i=1:b if d(i)h&jh A(k,i)=-1;A(k,j)=1;end if ih A(k,j)=1;end if ih&j=h A(k,i)=-1;en

12、dendx=A(:,(h+1):a);%除去已知点的系数，使其成为满秩矩阵y=P;z=L;5.2.3 残差计算残差也叫观测值的平差改正数，在参数平差中观测值的平差值可以直接用参数平差值计算出来，所以计算残差的目的不是用来改正观测值，而是主要用来进行精度估计。参加计算由 M 文件 ca_V 完成，计算结果存在数据V 中，并返回pvv值（用于计算单位权中误差）。1 算法以观测值序号k为循环变量，观测值循环，由式（5-2-1）计算各观测值的残差vk，第k次循环中所要进行的工作包括：（1）获得高差的起点号i和终点号j；（2）获得起点和终点的高程值Hi和Hj（3）计算残差vk Hj Hihk，存在Vk中

13、；（4）计算pvv。2 ca_V 函数源代码functionx y=ca_V(a,b,c,d,e,f)%计算高差改正数%V pvv=ca_V(m_Lnumber,startP,endP,Height,L,P)pvv=0.0;for i=1:a k1=b(i);k2=c(i);V(i)=d(1,k2)-d(1,k1)-e(1,i);endx=V;y=V.*V*f;5.2.3 精度估计与平差结果输出精度估计与平差结果输出由函数 printlevledata（）完成。该函数的平差结果写入结果文件。输出内容包括高程点名，高程平差值及其中误差、高差端点点名、高差平差值、高差改正数、高差平差值的中误差。M

14、 文件 printlevledata 的源代码如下：functionprintlevledata(pathName,Filename,m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,pname,height,startP,endP,L,V,P,pvv,m_mu,caheight,dx,Height,qii,Q)%输出结果Filename pathName=uiputfile(*.txt,saving file as);if isequal(Filename,0)|isequal(pathName,0),return;endstr=pathName Filename;

15、%依次输入结果到 txt 文件fid1=fopen(str,wt);fprintf(fid1,Total Observed:%dTotal points:%d Total known points:%dn,m_Lnumber,m_Pnumber,m_Knumber);fprintf(fid1,n Priori unit weight weight mean error:%8.3f,m_Sigma);fprintf(fid1,nn Known Height:n);for i=1:m_Pnumber fprintf(fid1,n%5d%8s,i,pnamei);fprintf(fid1,%8.4f

16、,height(i);endfprintf(fid1,nnThe difference of Observations n);for i=1:m_Lnumber fprintf(fid1,n%5d%8s%8s,i,pnamestartP(i),pnameendP(i);fprintf(fid1,%8.4f%8.4f ,L(i),1.0/P(i);endfprintf(fid1,nn Resluts of Least Squares:n pvv=%8.9f,pvv);fprintf(fid1,n u=%1f,m_mu);fprintf(fid1,nnn =The Height Adjustmen

17、gt and Precision=n);fprintf(fid1,n Pname ApproHeight V AdjustHeight RMSEn);for i=1:m_Pnumber fprintf(fid1,n%6s,pnamei);fprintf(fid1,%10.4f%8.4f,caheight(1,i),dx(1,i),Height(1,i),qii(1,i);end%8.4f%8.4ffprintf(fid1,nnn=The Observations Adjustmengt andPrecision=n);fprintf(fid1,n NO.SPoint EPoint DiffOb

18、servations V);fprintf(fid1,AdjustObservations WObservations RMSEn);for i=1:m_Lnumber k1=startP(i);k2=endP(i);ml=0;if k1m_Knumber ml=qii(1,k2);end if k1m_Knumber&k2m_Knumber&k2m_Knumberml=sqrt(Q(k1-m_Knumber,k1-m_Knumber)+Q(k2-m_Knumber,k2-m_Knumber)-2*Q(k1-m_Knumber,k2-m_Knumber)*m_mu;end fprintf(fi

19、d1,n%5d%5s%5s,i,pnamek1,pnamek2);fprintf(fid1,%8.4f%8.4f%8.4f%8.4f%8.4f,L(1,i),V(1,i),L(1,i)+V(1,i),P(1,i),ml);endfclose(fid1);if exist(str)=2 msgbox(save data successfully);else errordlg(Have not saved data,Failed);return;end5.35.3 水准网平差算例及分析水准网平差算例及分析要进行最小二乘平差计算，还需要编写一个主函数 M 文件，进行数据提取，连接上面提到的函数文件，

20、才能完成特定的平差任务。1 计算步骤（1）已知点提取处理。（2）计算近似高程。（3）计算法方程式。（4）阶段法方程，最小二乘平差。（5）计算残差。（6）最后成果输出，包括高程平差值及其精度、高差平差值及其精度。2 函数源代码%读取起始文件Filename pathName=uigetfile(*.txt,Txt Files|*.txt,choose a File);%打开文件%检查读取文件格式m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,unpnumber,P1,P2,P3,P4,pname=levelcheck(Filename,pathName);%提取已知点

21、数据for i=1:m_Knumber pname1i=P1i;point2(i)=str2num(cell2mat(P2(i);end%提取观测值及每个观测值的起终点号，路线长度for i=1:m_Lnumber cpname1i=P1i+m_Knumber;cpname2i=P2i+m_Knumber;point3(i)=P3(i+m_Knumber);point4(i)=P4(i+m_Knumber);end%按已知点点高程按着该点在点名数组中未知，将高程值赋给对应位置height=zeros(1,m_Pnumber);for i=1:m_Knumber c=strmatch(pname

22、1(i),pname);height(c)=point2(i);end%将观测值高差的起点点名转换成数组位置for i=1:m_Lnumber startP(i)=strmatch(cpname1(i),pname);end%将观测值高差的终点点名转换成数组位置for i=1:m_Lnumber endP(i)=strmatch(cpname2(i),pname);end%获取每条线路的高差值及其对应的权for i=1:m_Lnumber L(i)=point3(i);P(i)=point4(i);P(i)=1.0/P(i);end%计算近似高程caheight=calculateHO(m_K

23、number,m_Pnumber,m_Lnumber,height,startP,endP,pname,L);%计算法方程系数 A，权阵 GP，常数项 GLA GP GL=ca_ATPA(m_Pnumber,m_Lnumber,startP,endP,L,P,caheight,m_Knumber);%最小二乘平差dX=inv(A*GP*A)*A*GP*GL;%未知点的改正数dx=zeros(1,m_Pnumber);%所以点的改正数为定义为0%将改正数为零的未知点的改正数替换成计算值for i=(m_Knumber+1):m_Pnumber dx(1,i)=dX(i-m_Knumber,1);

24、end%计算未知点的最佳估值for i=(m_Knumber+1):m_Pnumber caheight(i)=caheight(i)+dX(i-m_Knumber,1);endHeight=caheight;%计算高差改正数V pvv=ca_V(m_Lnumber,startP,endP,Height,L,P);%计算验后单位权中误差m_mu=sqrt(pvv/(m_Lnumber-(m_Pnumber-m_Knumber);%计算未知点的点位中误差 Q=inv(A*GP*A);qii=zeros(1,m_Pnumber);%所有点的点位中误差定义为零 for i=(1+m_Knumber)

25、:m_Pnumber qii(1,i)=sqrt(Q(i-m_Knumber,i-m_Knumber)*m_mu;end%输出结果printlevledata(pathName,Filename,m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,pname,height,startP,endP,L,V,P,pvv,m_mu,caheight,dx,Height,qii,Q);3 最小二乘算例水准网如图 5-3-1 所示，共有7 个高程点，其中A、B、C、D为已知点，E、F、G为未知点，共观测了6 段高差，已知高程和观测高差分别见表5-3-1 和表 5-3-2，观测值的

26、每千米高差中误差为0.005m。试以此数据进行最小二乘平差。表 5-3-1 已知点数据点名 A BCD高程/m35.41845.71225.27024.678图 5-3-1表 5-3-2 观测高差No.123过程（1）准备数据文件7520.001A1 B1 P1 P2 P3A1 5.016B1 6.016A1 P11.3591.1A1 P22.0091.7B1 P10.3632.3B1 P21.0122.7P1 P20.6572.4起点AEE终点EBFh/m8.2282.0601.515S/km4.04.02.0No.456起点FGG终点P2P2P3h/m7.477S/km4.012.4172

27、.013.0002.0P1 P30.2381.4P3 B1-0.5952.6（2）算例函数Filename pathName=uigetfile(*.txt,Txt Files|*.txt,choose a File);m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,unpnumber,P1,P2,P3,P4,pname=levelcheck(Filename,pathName);for i=1:m_Knumber pname1i=P1i;point2(i)=str2num(cell2mat(P2(i);endfor i=1:m_Lnumber cpname1i=P

28、1i+m_Knumber;cpname2i=P2i+m_Knumber;point3(i)=P3(i+m_Knumber);point4(i)=P4(i+m_Knumber);endheight=zeros(1,m_Pnumber);for i=1:m_Knumber c=strmatch(pname1(i),pname);height(c)=point2(i);endfor i=1:m_Lnumber startP(i)=strmatch(cpname1(i),pname);endfor i=1:m_Lnumber endP(i)=strmatch(cpname2(i),pname);end

29、for i=1:m_Lnumber L(i)=point3(i);P(i)=point4(i);P(i)=1.0/P(i);endcaheight=calculateHO(m_Knumber,m_Pnumber,m_Lnumber,height,startP,endP,pname,L);A GP GL=ca_ATPA(m_Pnumber,m_Lnumber,startP,endP,L,P,caheight,m_Knumber);dX=inv(A*GP*A)*A*GP*GL;dx=zeros(1,m_Pnumber);for i=(m_Knumber+1):m_Pnumber dx(1,i)=d

30、X(i-m_Knumber,1);endfor i=(m_Knumber+1):m_Pnumber caheight(i)=caheight(i)+dX(i-m_Knumber,1);endHeight=caheight;V pvv=ca_V(m_Lnumber,startP,endP,Height,L,P);m_mu=sqrt(pvv/(m_Lnumber-(m_Pnumber-m_Knumber);Q=inv(A*GP*A);qii=zeros(1,m_Pnumber);for i=(1+m_Knumber):m_Pnumber qii(1,i)=sqrt(Q(i-m_Knumber,i-

31、m_Knumber)*m_mu;endprintlevledata(pathName,Filename,m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,pname,height,startP,endP,L,V,P,pvv,m_mu,caheight,dx,Height,qii,Q);（3）计算结果 Total Observed:6 Total points:7 Total known points:4 Priori unit weight weight mean error:0.005 Known Height:1 A1 35.4180 2 B1 45.7120 3

32、 C1 25.2700 4 D1 24.6780 5 E1 0.0000 6 F1 0.0000 7 G1 0.0000The difference of Observations 1 A1 E1 8.2280 4.0000 2 E1 B1 2.0600 4.0000 3 E1 F1 1.5150 2.0000 4 G1 F1 7.4770 4.0000 5 C1 G1 12.4170 2.0000 6 D1 G1 13.0000 2.0000 Resluts of Least Squares:pvv=0.000027000 u=0.003000 =The Height Adjustmengt

33、 and Precision=Pname ApproHeight V AdjustHeight RMSE A1 35.4180 0.0000 35.4180 0.0000 B1 45.7120 0.0000 45.7120 0.0000 C1 25.2700 0.0000 25.2700 0.0000 D1 24.6780 0.0000 24.6780 0.0000 E1 43.6480 0.0020 43.6480 0.0037 F1 45.1620 0.0010 45.1620 0.0045 G1 37.6830 -0.0010 37.6830 0.0028 =The Observatio

34、ns Adjustmengt and Precision=NO.SPointEPointDiffObservationsVAdjustObservationsWObservations RMSE 1 A1 E1 8.2280 0.0020 8.2300 0.25000.0037 2 E1 B1 2.0600 0.0040 2.0640 0.25000.0037 3 E1 F1 1.5150 -0.0010 1.5140 0.50000.0037 4 G1 F1 7.4770 0.0020 7.4790 0.25000.0045 5 C1 G1 12.4170 -0.0040 12.4130 0

35、.50000.0028 6 D1 G1 13.0000 0.0050 13.0050 0.50000.0028各项字母含义注释：Total Observed 观测值总数 Total points 总点数 Total known points 已知点数Priori unit weight weight mean error 验前单位权中误差 Known Height 已知高程The difference of Observations 观测值高差 Pname 点名 V 改正数 WObservations 观测权Resluts of Least Squares 最小二乘平差结果 ApproHeight 近似高程The Height Adjustmengt and Precision 高程平差值及其精度 AdjustHeigh 高程平差值RMSE中误差 The Observations Adjustmengt and Precision 观测平差值及其精度SPoint 起点 EPoint 终点 DiffObservations 观测高差 AdjustObservations 高差平差值