第三章水准网平差程序设计

1、第三章第三章 水准网平差程序设计水准网平差程序设计本章难点：本章难点：1、近似高程计算2、最短路线的计算3、误差方程及法方程的构建4、直接计算出法方程的系数矩阵BTPB和常数矩阵BTPL 该课程中所采用的平差模型为间接平差，即所选的独立参数的个数等于必要观测数，这样可以将每个观测值表示成这t个参数的函数，组成观测方程。 间接平差的函数模型：1 ,1 ,1 ,nttnndXBL平差时一般对所选参数取近似值 00)(LLdBXLlxXX0代入上式并令由此可得误差方程lxBVttnn1 ,1 ,式中l为误差方程常数项一、水准网间接平差算法概述一、水准网间接平差算法概述采用的平差准则为最小二乘法:mi

2、nPVVT 不难看出总共有n+t个参数，而所列的误差方程只有n个，故只采用误差方程是不能解求n+t个参数的，还需寻找其他的条件。按着最小二乘法法则，对t个参数求偏导，得出式子0PVBT0PlBxPBBTT由误差方程和求偏导后得到的式子，可唯一确定n+t个待求量。并有这两个 式子可以推出法方程，如为：PlBPBBxTT1)(从而可求得,将该式代入到误差即可求得改正数 。 v从而平差值结果为：xXXvLL0, 本课程中以未知点高程值作为未知参数，以高差作为观测值，以观测路线长度的倒数作为观测值的权，采用间接平差模型进行平差计算。其平差的一般步骤如下：其平差的一般步骤如下：1.根据平差问题选定未知参

3、数；2.根据观测值与未知参数之间的函数关系建立误差方程式，若误差方程是非线性方程，还要引入参数近似值，将误差方程线性化；3.由误差方程组成法方程；4.解算法方程，求取未知参数；5.精度评定。 平差程序设计不仅仅是“写程序”，测量程序设计还包括程序功能设计、平差模型选择、算法选择、数据接口设计等内容。一般应按数据处理和计算功能的划分，将网平差程序分为若干独立函数（或过程），每个函数（或过程）完成特定的计算或操作，当需进行某种平差时，再调用这些函数即可。 二、水准二、水准网平差程序网平差程序设计数据组设计数据组织及处理流织及处理流程程2.1 2.1 水准网平差中数据的组织水准网平差中数据的组织(

4、(数据文数据文件格式的设计件格式的设计) ) 水准网平差程序设计中数据的组织是一个关键问题。文件格式的设计就是规定数据文件中包含的内容、各类数据的先后顺序以及各种数据的书写格式。从程序设计者的角度来说，数据文件格式设计相当于告诉用户应当怎样准备数据文件，当然，程序中也必须按照约定的格式从文件中读取数据。数据文件由四个部分组成，其格式如下。1）基本信息部分：该部分仅占一行，其内容为已知点数已知点数N1N1，未知点数，未知点数N2N2，高差观测值个数，高差观测值个数NSNS。2）已知点高程信息：该部分占N1行，每行格式为已知点编号，该点的高程已知点编号，该点的高程（单位：米）。注意：在给控制点编号

5、时，先给待定点编号，然后注意：在给控制点编号时，先给待定点编号，然后给已知点编号。编号从给已知点编号。编号从1开始顺序编号。开始顺序编号。3）高差观测值信息：该部分占NS行，每行格式为测段编号，测段起点编号，测段终点编号，测段观测段编号，测段起点编号，测段终点编号，测段观测高差（单位：米），测段长度（单位：千米）。测高差（单位：米），测段长度（单位：千米）。4）点名信息：该部分占N（N=N1+N2）行，每行格式为水准点编号，该点名称水准点编号，该点名称其中点名长度不超过8个字符，即4个汉字。但最好最好用字符串命名用字符串命名。 以下为一具有两个已知点、三个待定点和七个测段的水准网的示例数据，打

6、开TXT文件文件后，内容如下。2 , 3 , 7 4 , 5.016 5 , 6.016 1 , 4 , 1 , 1.359 , 1.1 2 , 4 , 2 , 2.009 , 1.7 3 , 5 , 1 , 0.363 , 2.3 4 , 5 , 2 , 1.012 , 2.7 5 , 1 , 2 , 0.657 , 2.4 6 , 1 , 3 , 0.238 , 1.4 7 , 3 , 5 ,-0.595 , 2.6 1 ,P12 ,P23 ,P34 ,A5 ,B2.2 2.2 数据存储方案设计数据存储方案设计 这里的数据存储是指数据在内存中的存储。测量平差程序中用到的数据较多，一般先将磁

7、盘中的数据读到内存中，按照不同的类别有序地放在变量或数组中，然后再进行平差计算，最后将计算结果写到磁盘文件中。 要整洁、简练、实用、方便有关数据的输入或输出，必要的话要考虑和用户的交互，尽量不要画蛇添足。三、水准网平差界面设计三、水准网平差界面设计水准网平差输入界面一水准网平差输入界面一水准网平差输入界面二水准网平差输入界面二水准网平差输入界面三水准网平差输入界面三四、高差闭合差的计算与检查四、高差闭合差的计算与检查 为了检查水准高差的质量，水准网平差前一般进行附合水准路线闭合差计算与多边形闭合差计算，并进行闭合差的检核。1.附合水准路线高差闭合差的计算附合水准路线高差闭合差的计算 a）最短路

8、线搜索算法 “邻接点”：若一个点借助另一个点与目标点发生了联系，称另一个点是这个点的“邻接点邻接点”。如果一个点与目标点有直接联系，我们也把目标点称为这个点的邻接点。 显然邻接点不是唯一的，但是我们规定每个点只有一个邻接点；当网中每个点的邻点都唯一确定时，顺着邻接点又可找到邻点的邻接点，这样一层一层地找下去，最终找到目标点，每个点到目标点的路线是唯一确定的。 根据邻接点的概念可知，任何一点沿邻接点到目标点的路线长度是其邻接点到目标点的路线长度再加上该点与其邻接点之间的路线长度。当一个点有多个点可以作为邻接点时，只选择路线最短的一点作为邻接点，这样就可保证沿邻接点到目标点的路线是最短路线。 设p

9、1选择p2作为邻接点，P2到目标点的路线长度为S2, p1到p2的路线长度为S12，则p1到目标点的路线长度S1为S S1 1=S=S2 2+S+S1212 若还有一点P3也是P1的邻接点，则经由P3到目标点的路线长度为 S S1 1=S=S3 3+S+S1313 为了寻求最短路线，显然当S1SS1 1+S+S12 12 ( (条件条件1)1)p2可作为p1的邻接点的条件为：S S1 1SS2 2+S+S12 12 ( (条件条件2)2)若条件1成立时,Neighbor(p2)=p1,S(p2)=S1+S12；若条件2成立时,Neighbor(p1)=p2,S(p2)=S1+S12。（4）当水

10、准网的规模较大时，在循环访问每一个观测值后某些点到目标点的路线还不是最短路线，甚至有些点还没有找到邻接点，需要转到（3）继续为寻找路线最短的邻接点，直到所有的观测边不再满足条件条件1 1和条件和条件2 2，这时表明每个点都已经找到了到目标点最短的邻接点，搜素工作结束。最短路线搜索函数FindShortPath（p as integer, exclude as integer, neighbor() as integer,diff() as double, S() as double） 从一个已知点出发，用观测高差依次推算其从一个已知点出发，用观测高差依次推算其它各点的高程，最后闭合到另外一个已

11、知点上，它各点的高程，最后闭合到另外一个已知点上，闭合点上推算高程和已知高程值之差称之为附合闭合点上推算高程和已知高程值之差称之为附合水准路线闭合差，简称路线闭合差水准路线闭合差，简称路线闭合差。设k1、k2为两个已知点，已知高程值分别为H1、H2，两点之间的推算路线由高差h1,h2hn，各测段路线长度分别为S1,S2 Sn，路线闭合差为式中，正负号取决于高差起始点到终点的方向与推算路线是否一致。b)路线闭合差计算算法路线闭合差计算算法2211HhhhHWn 闭合差的限差为:假如网中共有m个已知点，将m个已知点两两组合，计算m(m+1)/2条附合路线的闭合差，每条路线按最短路线进行计算。 差。

12、为每千米观测高差中误限其中，SSSWn212 在水准网中，由观测高差相连接可以形成闭合多边形，理论上构成闭合多边形的各观测高差之和应该为0，由于观测误差的存在，高差之和一般不等于0，闭合多边形的观测之和称为环闭合多边形的观测之和称为环闭合差。闭合差。 一般情况下，设构成闭合环的观测高差h1,h2hn，环闭合差为 式中，正负号取决于高差起点到终点的方向与推算路线的方向一致。c)环闭合差计算算法环闭合差计算算法nhhhW21环闭合差的限差为实际计算时，一般只计算最小独立环的闭合差。最小闭合环应该满足如下条件：（1）多边形环应该是相互独立（线性无关）的，即任何一个多变边不能由其它多边形线性组合而得到

13、。满足独立条件可以避免重复计算，也满足独立条件可以避免重复计算，也可以避免遗漏。可以避免遗漏。设水准网中有n个观测高差、t个高程点，那么独立闭合环的个数为r=n-t+1保证闭合环独立的方法是，每个新环找到的闭合环都有前面找到的闭合环中不曾使用的观测值。差。为每千米观测高差中误限其中，SSSWn212 （2）在多边形环独立的情况下，闭合环的边长最短。 关于最小环问题可以用最短路径的方法解决。若仅通过一个观测值即可将两个点p1到p2连接起来，构成一条路线，可把这种路线称为p1到p2的直接路线。如果直接路线能同别的观测值构成闭合环，那么从p1到p2应有两条或两条以上的路线。 如果选定一条直接路线，再

14、从非直接路线中找到p1到p2的最短路线，把这条最短路线与直接路线连接起来一定是包含直接路线的最小环。为保证找到的闭合环是独立的，只能将未参加过前面环闭合差计算的观测作为直接路线。五、待定点近似高程计算五、待定点近似高程计算 待定点高程计算是为了得到未知点的高程近似值，在求取未知点近似高程的时候，要注意测段的方向及相关正负号。既可以进行人工计算后输入到程序中，也可以直接将原始数据直接读入到程序中，根据相关算法让程序自动计算。计算原理如下式：测段已知近似hHH水网观测数据（测段）的组织：水网观测数据（测段）的组织：测段号测段号 起点起点 终点终点 测段高差测段高差(m) 测段长测段长(km) 测站

15、数测站数1 A B 5.835 1.52 E C 1.006 0.83 D C 7.384 2.14 B C 3.782 3.25 A D 2.270 1.76 A C 9.64 1.37 A F 0.003 4.1水准网近似高程计算算法水准网近似高程计算算法1、定义存储水准网高程的数组H()，将已已知点（原有高程和计算出的近似高程都可知点（原有高程和计算出的近似高程都可视为已知点）视为已知点）存到H()数组中；2、遍历所有测段（即从第一个测段到最后一个测段）； a）判断测段中起点为已知点，判断终点为未知点，若同时成立，则利用公式计算该测段中未知点的高程，并存储到H()中； b）再次遍历测段，

16、找到满足这些条件的测段 （其起点未知点，终点是已知点），若同时满足，则进行起点近似高程的计算，把计算出的近似高程存到H()数组中；3、判断计算出近似高程的个数是否小于未计算出近似高程的个数是否小于未知点高程点的个数知点高程点的个数，若是再次进行循环，若否，计算完毕（即水准网中未知高程点的近似高程已全部计算出来），水准网中所有未知点的高程存放到H()数组中了。六、六、误差方程的组成误差方程的组成 对于水准网平差，关键就是误差方程的构建，即相关矩阵的生成，这里所指的矩阵是误差方程系数矩阵B和常数矩阵L。如何来组成这两个矩阵呢？还是一个方程一个方程的建立，最后将相应值存储到矩阵B和L的对应元素中。系

17、数矩阵系数矩阵B B构建算法描述：构建算法描述：1、定义一个存储系数B的数组B()，其行数为观测测段数，列数为所选必要参数个数；定义一个存储常数阵L的数组L()，其行数为测段个数，列数为1；2、取出第j个测段（测段数据存储在测段信息数组中的），首先判断该测段起点是否为已知点（也即其在已知点数组H()中的下标位于起算高程点个数内，如：起算高程点个数为5，起点在已知数组中的下标为i，当i小于等于5时，说明该起点为已知点，若i大于5，则说明该起点为未知点，同样终点是什么类型的点也通过该方法判断），记下其在a()中的下标i1（原始已知点数为i0）；判断该测段终点是什么类型的点记下其在a()中的下标i2

18、； a)、若该测段的起点为未知点，且终点也为未知点，则矩阵B的第j行、第（i1-i0）列的元素为-1，第j行、第（i2 i0）列的元素为1；那么第j行其它元素都为 0； b)、若该测段的起点为已知点，且终点为未知点，则矩阵B的第j行、第（i2-i0）列的元素为1，该行其它元素都为0； c)、若该测段的起点为未知点，且终点为已知点，则矩阵B的第j行、第（i1-i0）列的元素为-1，该行其它元素都为0；3、根据每测段起点和终点点的类型，就可以计算出这个误差方程的常数项，把计算出的数值放到L()中去就可以了；4、遍历完所有的测段，即可构建出矩阵B和常数项矩阵L，即其所有的元素都存储到数组B()和L(

19、)中了。七、法方程的组建七、法方程的组建有了B和L，基于最小二乘原理，对所选参数求偏导后的式子，推出法方程，即：0PlBxPBBTT直接计算法方程系数阵和常数阵（累加法）直接计算法方程系数阵和常数阵（累加法） 组成法方程是间接平差的关键步骤，由间接平差模型可知，在误差方程系数矩阵B，误差方程常数项以及观测值的权P确定时，法方程系数矩阵和自由项向量分别为BTPB和BTPL，组成法方程就是矩阵的乘积。但在水准网平差程序设计时，为了节省存储单元，可不必为误差方程矩阵B开辟存储数组，权阵也只保存对角线上的元素，组成法方程的计算就不是矩阵乘积运算了，而是根据观测值的起始点号、终点号、高差值和观测值的权直

20、接计算出法方程的系数阵BTPB和自由项向量BTPL。000000kkkkkKTkkkkkkKTklplplPBPpppBPB 应将应将p pk k、- -p pk k累加到法方程系数矩阵的什么位累加到法方程系数矩阵的什么位置？若第置？若第K K个观测高差的起点点号和终点点号分个观测高差的起点点号和终点点号分别为别为i i、j j，则需将，则需将p pk k累加到法方程系数矩阵的第累加到法方程系数矩阵的第i i行第行第i i列和第列和第j j行第行第j j列，将列，将- -p pk k累加到法方程系累加到法方程系数矩阵的第数矩阵的第j j行和第行和第i i列。同理，应将列。同理，应将p pk k

21、l lk k和和- -p pk kl lk k分别累加到法方程常数项的第分别累加到法方程常数项的第j j行和第行和第i i行中。行中。 组成法方程的过程就是以观测值总数为循环组成法方程的过程就是以观测值总数为循环界，逐个观测值计算累加项的过程，在第界，逐个观测值计算累加项的过程，在第k 个循个循环中，所做的工作如下：环中，所做的工作如下：（1）获得高差的起点点号）获得高差的起点点号i和终点号和终点号j；（2）获得起点和终点的高程值）获得起点和终点的高程值Hi和和Hj；（3）计算误差方程的自由项）计算误差方程的自由项lk=Hj-Hi-hk；（4）将）将pk、pklk累加到法方程系数阵和法方程自累

22、加到法方程系数阵和法方程自由项的相应单元中。由项的相应单元中。 可以将法方程系数矩阵和法方程常数项矩阵的可以将法方程系数矩阵和法方程常数项矩阵的计算定义成一个函数计算定义成一个函数Ca_ATPA()。void CLevelingAdjust:ca_ATPA() int t=m_Pnumber; for(int i=0; it*(t+1)/2; i+) ATPAi=0.0; for(i=0; it; i+) ATPLi=0.0;for(int k=0; km_Lnumber; k+) int i=StartPk;int j=EndPk; double Pk=Pk; double Lk=Lk-(H

23、eightj-Heighti); ATPLi-=Pk*Lk;ATPLj+=Pk*Lk;ATPAij(i,i)+=Pk;ATPAij(j,j)+=Pk;ATPAij(i,j)-=Pk; 八、水准网平差定权八、水准网平差定权 水准网平差中还有涉及到权阵P的确定，有了权阵P ，和我们上面所构建出来的B和l ，就可以很容易的进行参数改正值和观测值改正值的解算了。 我们一般怎样来定权?设有观测值Li(i=1,2,n)它们的方差 为若选定一常数 , 定义 ，并称P为观测值Li的权。对于一组观测值，其权之比等于相应方差的倒数之比。这表明，方差越小，其权越大；或者说精度越高，其权越大。因此权可以作为比较观测值

24、之间的精度高低的一种指标。), 2 , 1(2nii0220iip结论:（1） 选择了一个 ，即对应一组权；（2） 一组观测值的权，其大小是随 的不同而异，但不论选何值，权之间的比例关系始终不变；（3）为了使全能够起到比较精度的作用，在同一问题中只能选取一个 值，不能同时选取几个 ，否则就破坏了权之间的比例关系；（4）只要给出了一定的观测条件，不一定要知道每公里观测高差精度的具体数值，就可以确定出权的数值。0000 水准网平差中定权，一般二种方式，即以距离距离定权定权和测站数测站数定权。 一般观测量之间是相互独立的，故在构建权矩阵P时，只有对角线上的元素 非零，非对角元素都为零。 注意采用以上

25、方法定权时，必须满足“每测每测站观测高差精度相等站观测高差精度相等”或者“每公里观测高差每公里观测高差精度相等精度相等”这一前提条件。九、精度的评定九、精度的评定 间接平差采用了与条件平差不同的数学模型，但它们是在相同的平差原理下进行的，所以这两种方法的计算出的结果是相同的。中误差的估值为： 而对于 的计算，将误差方程代入即可。平差参数的协方差阵 ，而对于平差函数，其权函数式为： ，协因数 方差为：tnPVVT0PVVT12020BBXXXXNQDxFdTFNFFQFQBBTXXT1QD20水准网如图水准网如图1 1所示，共所示，共7 7个高程点、个高程点、1010个测段观个测段观测高差，观测

26、值的每千米观测中误差为测高差，观测值的每千米观测中误差为0.001m0.001m，已知高程和观测高差分别见表，已知高程和观测高差分别见表1 1及表及表2 2。水准网平差程序设计实验中用例水准网平差程序设计实验中用例水准网路线水准网路线 已知数据和观测数据已知数据和观测数据点名点名高程高程/m/mA0.000F11.414表表1 1 已知高程点已知高程点表表2 2 观测高差与路线长度观测高差与路线长度编号编号起点点号起点点号终点点号终点点号观测高差观测高差/m/m距离距离/m/m1AB73.79520.42AD14.00518.83AG14.16715.44CB71.9498.95DB59.78

27、014.26CD12.15912.87CE15.3649.88FE5.79719.69GE3.04415.110DG0.16910.0该水准网平差数据文件（该水准网平差数据文件（data.txtdata.txt）如下：）如下：10,7,2,0.00110,7,2,0.001A,0.000A,0.000F,11.414F,11.414A,B,73.795,20.4A,B,73.795,20.4A,D,14.005,18.8A,D,14.005,18.8A,G,14.167,15.4A,G,14.167,15.4C,B,71.949,8.9C,B,71.949,8.9D,B,59.780,14.2

28、D,B,59.780,14.2C,D,12.159,12.8C,D,12.159,12.8C,E,15.364,9.8C,E,15.364,9.8F,E,5.797,19.6F,E,5.797,19.6G,E,3.044,15.1G,E,3.044,15.1D,G,0.169,10.0D,G,0.169,10.01 1、水准网程序设计中变量及数组说明、水准网程序设计中变量及数组说明 n网中已知点数：m_knPnumberm_knPnumbern网中总点数：m_Pnumberm_Pnumbern先验权中误差：M Mn单位权中误：m_mum_mun存储网中观测测段高差个数：m_GaochaNm_G

29、aochaNn存储测段高差数组：gaocha()gaocha()n测段起点编号数组：StartP()StartP()n测段终点编号数组：EndP()EndP()n测段长度数组：Li()Li()n存储已知高程与近似高程值数组：H()n误差方程系数矩阵数组：B()n误差方程常数项数组：L()n权阵数组：P()n存储值的数组：p_vvn高程改正值数组：dx()n法方程系数矩阵数组：Q()n观测高差改正值数组：v()2、数据文件读入数据文件读入 为了从数据文件中读取数据，首先利用VB6.0中的Open函数打开该数据文件，接着，按着数据文件中数据内容的顺序和格式顺序和格式，先读取网的概况信息数据，再读取

30、已知高程数据，最后读取观测数据。1）读取网的概况数据，将总点数、已知点数、观测值总数分别存储到变量m_knPnumberm_knPnumber、 m_Pnumberm_Pnumber及及m_GaochaNm_GaochaN中。中。2）动态定义数组。观测高差起点点号数组StartP()、高差终点数组EndP()、观测高程数组gaocha()、残差数组v()、观测值权P()与观测高差数有关；高程数组H()注意：动态数组是先声明，再动态定义其大小。注意：动态数组是先声明，再动态定义其大小。3 3、近似高程的推算（、近似高程的推算（ca_Hca_H（）（）近似高程的计算算法可参照该章节前面的内容。在编

31、写程序过程中，最好完成特有任务和功能的语句，定义成函数（或过程），以方便其调用。Sub ca_H() 水准网中待定点近似高程计算水准网中待定点近似高程计算Dim k1 As Integer, H_N As Integer H_N为计算出近似高程点个数为计算出近似高程点个数Dim c1 As Integer, c2 As IntegerH_N = 0Do While H_N (m_Pnumber - m_knPnumber) For k1 = 1 To m_GaochaN c1 = StartP(k1) 观测高差起点观测高差起点 c2 = EndP(k1) 观测高差终点观测高差终点 If H(c

32、1) = 0 And H(c2) 0 Then H(c1) = H(c2) - gaocha(k1) H_N = H_N + 1 End If If H(c1) 0 And H(c2) = 0 Then H(c2) = H(c1) + gaocha(k1) H_N = H_N + 1 End If NextLoopEnd Sub4 4、误差方程系数矩阵、误差方程系数矩阵B B和常数项矩阵和常数项矩阵L L的生成的生成Sub GouJianBL(H() As Double, StartP() As Integer, EndP() As Integer) 构建误差方程系数矩阵构建误差方程系数矩阵B B与与L L * * * * * * * * * * * * * * * * * * * * * * *H()H()为近似高程数组，为近似高程数组， StatPStatP()()为起点点号数组，为起点点号数组，EndPEndP()()为终点点号数组为终点点号数组 * * * * * * * * * * * * * * * * * * * * * * *Dim k1 As Integer, k2 As Integer, k3 As Integer, c1 As