第五章水准网程序设计

第五章水准网程序设计5.1概述水准网是为了确定地面点的高程而布设的控制网，网中的观测值是高程控制点之间的高差。为了统一全国的高程系统，我国采用黄海平均海水面作为全国高程系统的基准面，在该面上的任一点，其高程为零。水准网中的任一点高程以及点与点之间的高差应属于正常高系统。但是，水准测量中水准仪逐站测量并累加得到的原始高差并非正常高差，对于高精度的水准网，应在原始数据基础上添加尺长改正、正常水准网面不平行、球气差改正等系统改正，才能得到正常高差，受篇幅所限，本章不讨论观测值的归算问题，本章提到的观测值均假定是已经归算后的正常高差观测值。由于观测值存在误差，实际工作中需要若干已知点作为平差的基准控制点。水准网平差的目的就是求解各观测值的最佳估值及评定未知量的精度。本章开始，我们将开始接触到测量数据处理程序的设计与开发。测量数据因其数据量大、繁杂等特殊性，往往需要频繁对矩阵进行处理计算。然而在MATLAB软件中，矩阵的计算变得格外简单，这为广大测绘工作者提供了实用性强、简单方便的数据处理平台。平差程序就是将平差计算过程程序化，综合考虑，选择参数平差模型作为水准网平差的主要模型。测量程序设计一般包含程序功能设计、平差模型选择、算法选择等内容。在本章的结构的组织上，先对水准原理进行简要分析，在此基础上按不同功能模块设计函数，最后对程序进行分析验证。5.2水准路线处理程序设计5.2.1水准网平差函数类设计1.类设计Main%函数主体程序calculateHO（）%计算近似高程ca_ATPA()%计算法方程系数A，权阵GP，常数项GLca_V()%计算高差改正数Printlevledata()%输出结果2.函数说明具体各个成员变量的含义在后面程序中会有注释，在此对主函数main做简要强调。主函数main中包括数据的读取、检查数据格式、提取相关观测值、各类成员函数和最小二平差。是本程序运行的主要M文件。5.2.2原始数据文件格式设计水准网平差所需的数据需要从txt文本中提取，本节主要介绍原始数据文件的内容与格式，需要说明的是，原始数据文件的设计并没有统一的、严格的标准，在方便程序设计的基础下，程序设计者可以自由设计数据文件的格式。我们将网的数据分为三类：网的概况信息、已知数据、观测数据。网的概况信息包括总点数、已知点点数、观测者总数、验前单位权中误差。已知数据包括已知点名、已知高程值。观测数据包括高差起点点名、高差终点点名、观测值高差值、线路长度。下面结合实例说明原始文件的具体数据格式图5-2-1为典型的水准网表5-2-1已知点数据图5-1水准网示意图表5-2-2观测高差No.起点终点h/mS/kmNo.起点终点h/mS/km1AP11.3591.15P1P20.6572.42AP22.0091.76P1P30.2381.43BP10.3632.37P3B-0.5952.64BP21.0122.7利用以上数据进行水准网平差，数据文件的内容如下：7520.001A1B1P1P2P3A15.016B16.016A1P11.3591.1A1P22.0091.7B1P10.3632.3B1P21.0122.7P1P20.6572.4P1P30.2381.4P3B1-0.5952.6格式说明：（1）第一行为网的概况信息：观测值总数、总点数、已知点总数、验前单位权中误差。点名高程/mA5.160B6.016（2）第二行是所有点点名依次书写。（3）第三行是对应已知点点名、高程（单位为m），依次往下写、列出所有一直点名和高程。（4）接下来是观测高差：高差起始点名、高差终点点名、高差观测子和高差线路长度。网中点名需要字母和数字结合，否则点名读不出来，点名之间空格；实际平差时，水准网规模和数据可能会不同，但是只要按照上面的格式和顺序输入相关数据即可用本章的程序。5.2.2数据的存储1点数据的存储平差程序用到的数组中，点名、高程值、高程改正数等数组与高程点名一一对应，这点数据称为点数据。点数据数组的总长度为总点数。高程值用数组Height保存，高差改正数用数据dX存放，近似值存放在数组caheight中。点名用数组pname数组存放。pname数组不仅存放数组名，同时也可看做点名地址数组。点数据数组存放顺序必须保持严格的一致，例如某点点名地址存在pname数组中的第3单元，即pname（3）中，该点的高程值也就存放在Height（3）中，高程改正数也存放在dX（3）中。当点数据存放一致时，可以更快捷的进行点数据存储，更利于程序设计。2观测数据的存储观测数据包括高差观测值、高差起始点的点名、高差的权、高差的改正数，它们与观测值一一对应，数组长度都等于总点数。高差起点点名和终点点名在程序中的作用是缺点观测值的起始点号和终点点号，有了点号才能访问点数据数组，进行相应的计算。利用函数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判断某点是否需要计算近似高程，当该点近似高程计算出来之后就会被取代。第二步，计算高程值。近似高程计算的基本思路是，遍历观测值，找到每一个观测值起点号、终点号，根据点号从Height数组中提取起点和终点的高程值,当高差的一端是高程已知点，另一端是未知点时，就由已知点高程和观测高差算出未知点高程，并将计算结果赋值给高程数组，如果两端都是未知点，就跳到下一循环。在遍历观测值时，排在前面的观测值总是最先用于高程近似值计算，当观测值的排序不同时，近似高程计算用到的观测值也不同，计算结果也不同。有时候一次遍历观测值可能还有部分点的高程值仍未计算出来，还需要再次遍历观测值，直至算出所有高程值为止。2函数文件calculateHO（）代码functiony=calculateHO(a,b,c,d,e,f,g,h)%计算近似高程%caheight=calculateHO(m_Knumber,m_Pnumber,m_Lnumber,height,startP,endP,pname,L);jj=0;fori=(a+1):bd(1,i)=-9999.9;%未知点的初始高程为-9999.9endfori=1:(b-a)%循环的次数不超过（m_Pnumber-m_Knumber）forj=1:c%按照观测组进行循环K1=e(j);%起点点号K2=f(j);%终点点号ifd(1,K1)-9999.0&&d(1,K2)-9999.0%起点未知，终点已知d(1,K2)=d(1,K1)+h(1,j);jj=jj+1;elseifd(1,K1)-9999.0&&d(1,K2)-9999.0%起点已知，终点未知d(1,K1)=d(1,K2)-h(1,j);jj=jj+1;endendy=d;ifjj==(b-a)breakendifi(b-a)%生成产生错误的txtfid=fopen('resultfid.txt','a+');fprintf(fid,'\n\n下列点无法计算出概略高程:\\n');fori=1:bifd(i)-9999.0fprint(fid,'\n%s',f{i});endendfclose(fid)errordlg('computationFailed!');breakendEnd5.2.3组成法方程组成法方程式是参数平差的关键步骤，由参数平差模型可知，在误差方程系数矩阵A，误差方程自由项l及观测值权阵P确定时，法方程系数矩阵和自由项分别为PAAT和PlAT。计算法方程的ca_ATPA（）代码如下：function[xyz]=ca_ATPA(a,b,c,d,e,f,g,h)%[AGPGL]=ca_ATPA(m_Pnumber,m_Lnumber,startP,endP,L,P,caheight,m_Knumber)%A为系数矩阵，GP为权阵，GL为常系数矩阵A=zeros(b,a);R=zeros(b,b);L=zeros(b,1);%获取观测值的常系数阵fork=1:bi=c(1,k);j=d(1,k);L(k,1)=e(1,k)-(g(1,j)-g(1,i));end%获取观测值的权阵fori=1:bP(i,i)=f(1,i);end%获取法方程的系数矩阵fork=1:bi=c(1,k);j=d(1,k);ifih&&jhA(k,i)=-1;A(k,j)=1;endifi=h&&jhA(k,j)=1;endifih&&j=hA(k,i)=-1;endendx=A(:,(h+1):a);%除去已知点的系数，使其成为满秩矩阵y=P;z=L;5.2.3残差计算残差也叫观测值的平差改正数，在参数平差中观测值的平差值可以直接用参数平差值计算出来，所以计算残差的目的不是用来改正观测值，而是主要用来进行精度估计。参加计算由M文件ca_V完成，计算结果存在数据V中，并返回[pvv]值（用于计算单位权中误差）。1算法以观测值序号k为循环变量，观测值循环，由式（5-2-1）计算各观测值的残差kv，第k次循环中所要进行的工作包括：（1）获得高差的起点号i和终点号j；（2）获得起点和终点的高程值iH和jH（3）计算残差kijkhHHv，存在kV中；（4）计算pvv。2ca_V函数源代码function[xy]=ca_V(a,b,c,d,e,f)%计算高差改正数%[Vpvv]=ca_V(m_Lnumber,startP,endP,Height,L,P)pvv=0.0;fori=1:ak1=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文件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)%输出结果[FilenamepathName]=uiputfile('*.txt','savingfileas');ifisequal(Filename,0)||isequal(pathName,0),return;endstr=[pathNameFilename];%依次输入结果到txt文件fid1=fopen(str,'wt');fprintf(fid1,'TotalObserved:%dTotalpoints:%dTotalknownpoints:%d\n',m_Lnumber,m_Pnumber,m_Knumber);fprintf(fid1,'\nPrioriunitweightweightmeanerror:%8.3f',m_Sigma);fprintf(fid1,'\n\nKnownHeight:\n');fori=1:m_Pnumberfprintf(fid1,'\n%5d%8s',i,pname{i});fprintf(fid1,'%8.4f',height(i));endfprintf(fid1,'\n\nThedifferenceofObservations\n');fori=1:m_Ln