序贯平差算法程序设计与实现

摘要近些年来，空间科学、信息科学和计算机技术迅速发展，测绘工程技术学科逐渐融合了空间技术以及计算机技术等现代技术的优点，其测量技术和手段也随之得到创新。在工程项目建设过程中，测绘工程作为基础性工作直接影响着工程施工的整体质量，但误差总是存在着的。在测量平差中，当观测数据量较少，则采用整体平差，若观测数据量较大，且于不同时期得到，整体平差需反复解算陈旧的观测数据，必然会效率低下，而序贯平差的提出刚好解决了这在一问题。序贯平差公式简单且规律性强，该方法利用前期的平差结果进行数据的处理，不需要存储大量的历史观测数据，也不需求大矩阵的逆就可获得与整体间接平差相同的平差效果。本文就针对这种平差方法，以C语言作为编程语言，由一个单独的函数即可实现，函数的功能单一，程序结构简单。关键词：近代平差；序贯平差；矩阵逆；间接平差ABSTRACTInrecentyears,withtherapiddevelopmentofspacescience,informationscienceandcomputertechnology,surveyingandmappingengineeringtechnologyhasgraduallyintegratedtheadvantagesofspacetechnologyandcomputertechnology,anditsmeasurementtechnologyandmeanshavealsobeeninnovated.Intheprocessofprojectconstruction,surveyingandmappingengineering,asabasicwork,directlyaffectstheoverallqualityofprojectconstruction,buttherearealwayserrors.Inthemeasurementadjustment,whentheamountofobservationdataissmall,thewholeadjustmentisadopted.Iftheamountofobservationdataislargeandobtainedindifferentperiods,thewholeadjustmentneedstosolvetheoldobservationdatarepeatedly,whichisboundtobeinefficient.Theproposalofsequentialadjustmentjustsolvestheproblem.Thesequentialadjustmentformulaissimpleandregular,thismethodusesthepreviousadjustmentresultsfordataprocessing,doesnotneedtostorealargenumberofhistoricalobservationdata,anddoesnotneedtheinverseofthelargematrixtoobtainthesameadjustmenteffectasthewholeindirectadjustment.Inthispaper,aimingatthisadjustmentmethod,Clanguageisusedastheprogramminglanguage,whichcanberealizedbyasinglefunction.Thefunctionofthefunctionissingleandtheprogramstructureissimple.Keywords：Modernadjustment，Sequentialadjustment，Matrixinverse，Indirectadjustment目录1绪论 11.1研究背景 11.2研究意义 11.3国内外发展现状 21.4研究内容与章节安排 22近代平差方法 42.1秩亏自由网平差 42.2附加参数系统的平差 122.3最小二乘滤波 152.4序贯平差 213序贯平差算法程序设计与实现 273.1水准网概况 273.2整体法平差与序贯法平差 273.3程序实现 314结论与展望 38参考文献 39致谢 40错误!超链接引用无效。1绪论 3错误!超链接引用无效。1.1研究背景 3错误!超链接引用无效。1.2研究意义 3错误!超链接引用无效。1.3国内外发展现状 4错误!超链接引用无效。1.4研究内容与章节安排 4错误!超链接引用无效。2近代平差方法基本理论知识 6错误!超链接引用无效。2.1观测误差 6错误!超链接引用无效。2.2函数模型和数学模型 8错误!超链接引用无效。2.3参数估计与最小二乘原理 10错误!超链接引用无效。3近代平差方法 15错误!超链接引用无效。3.1秩亏自由网平差 15错误!超链接引用无效。3.2附加参数系统的平差 23错误!超链接引用无效。3.3最小二乘滤波 26错误!超链接引用无效。3.4序贯平差 32错误!超链接引用无效。4算法程序设计与实现 38错误!超链接引用无效。4.1水准网概况 38错误!超链接引用无效。4.2整体法平差与序贯法平差 38错误!超链接引用无效。4.3程序实现 42错误!超链接引用无效。5结论与展望 46错误!超链接引用无效。参考文献 47错误!超链接引用无效。致谢 48山西大同大学建筑与测绘工程学院2020届本科生毕业设计1绪论1.1研究背景在当前社会，人们对于计算机技术和程序设计的认识和发展已经在不断提高，与此同时，效率低下的纯人工或手工作业的测绘工作正逐渐向着高效率的计算机方式稳步迈进。电子计算机的应用，使得测量平差作业模式发生了翻天覆地的变化，把平差计算者从繁重的、枯燥乏味的数字计算机中彻底解放了出来。本课题即利用实例数据通过C语言来实现测绘中的平差工作，并且采用序贯平差的方法来进行平差计算。C语言是一种过程式语言，有十分丰富的运算符集合和控制语句，支持多种数据类型，运算范围广，因此广泛用于底层开发，具有高效、灵活、可移植等特点，还能以简易的方式编译、运行。用C语言开发的项目，文件结构和程序结构都非常清晰，生成目标程序指令质量高，程序的执行效率也高，不仅适用于设计和编写编译器这类大型系统程序，又适用于编写小的控制程序。目前，C语言编译器普遍存在于各种不同的操作系统中，例如Windows下广泛使用的VisualStudio集成的编译器和Linux下广泛使用的GCC。近25年来，C语言已经成为最广泛使用的编程语言。1.2研究意义测量平差的理论和方法已经发展到了一个新的阶段，形成了内容丰富的近代平差。其主要特点是，观测值的概念广义化了，扩展了经典平差的数学模型，从处理随机独立的观测数据，发展到处理随机相关的观测数据；从只对满秩问题的平差，发展到具有任何秩的平差问题；从仅处理随机变量，发展到一并处理随机变量和具有各态历经性的平稳随机函数；从侧重于函数模型的研究，发展到也重视随机模型的研究；从不顾及模型误差，发展到顾及模型误差；从无偏估计，发展到也重视有偏估计的研究；从几何数据和物理数据分开处理，发展到几何数据和物理数据综合处理；从二维平差，发展到三维平差；从静态平差，将逐步发展到动态平差；从按经验设计大地网，发展到最优的设计大地网。在一个工程项目中，随着工程项目的不断开展，需要不断地布设高程点来满足测量工作的需要。对于传统的平差方法就是每次布设完毕后把网中的所有点进行整体平差，水准网不断扩大，于是，观测的数据也越来越多。计算时方程矩阵的阶数越来越大，由此，计算起来是十分麻烦的。序贯平差可以分组以进行平差，使得即使对于超大型的网，也不会出现观测数据庞大，计算繁琐困难的情况，这对于工程项目的建设有着必要的意义，本文通过一个简单的C语言程序，有针对性地实现了通过程序算法就可以直接得到待估参数的平差结果和单位权中误差，对于减少工作人员的计算量有着很大的研究意义。1.3国内外发展现状经典测量平差是C.F.高斯1794年提出的。随着测量手段的逐渐精密和现代化,从20世纪40年代以来,测量平差的理论和方法有了很大的发展,出现了一些新的测量平差方法,形成了内容丰富的近代测量平差。序贯平差是柯尔莫哥（A.H.Kolmogorov）和维纳（N.Wiener）提出的，也叫做逐次相关间接平差。其是利用了新旧观测值的全部信息，通过规律性很强的递推公式把新旧观测数据对未知参数估计的影响视为等影响原则进行参数估计。1960年，卡尔曼进一步发展了这一方法，称为卡尔曼滤波。而在工程测量领域，则主要应用的是参数不随时间变化的静态卡尔曼滤波，即序贯平差。2018年，山西建筑中李军亮在“某高程控制网的静态卡尔曼滤波计算”一文中以某高程控制网为例，进行了网中待定点高程的计算。分别采用两种平差方法，即经典的间接平差方式以及序贯平差的方式，进行了数据计算。1.4研究内容与章节安排本文主要介绍了经典平差的基本内容以及几种近代平差方法的基本理论，应用序贯平差方法设计了一个算法程序，实现了通过序贯平差方法来计算参数平差值以及单位权中误差。第一章，简要地介绍了测量平差的发展情况，以及近代平差方法的发展以及意义，序贯平差方法的用于高程网计算的研究意义。第二章，简单地列出了几种近代平差方法。如秩亏自由网平差，附加系统参数的平差，最小二乘滤波以及序贯平差方法的原理。第三章，利用VisualStudio2019平台基于C语言实现由序贯平差方法，求得单位权中误差，并与整体法平差进行比较，对比分析结果。第四章，结论与展望，总结得出结论，并对不足和未来进行展望。

。2近代平差方法2.1秩亏自由网平差在经典的平差内容里中，均以已知的起算数据作为基本起算条件，把控制网限制在已知的数据范围内。譬如在水准网里最基本的已知条件就是已知的水准网中某一个点的高程，平面控制网基础的已知条件是一个点的坐标、一条边的边长和一条边的方位角。如若控制网内有充足的起算数据，那么这时候的误差方程就可以写成：（3（2-1-1）其式子系数矩阵是列满秩矩阵，具体地说秩为。根据最小二乘准则，可获得法方称为（3（2-1-2）式中，。因其系数矩阵的秩为。故此是满秩矩阵，就是可逆矩阵，存在逆。因此具有唯一解，即（3（2-1-3）当控制网内没有任何起算数据或者是缺少部分的起算数据时，这时网中全部点都是待定点，假定未知参数的个数是，未知参数不是观测量时，误差方程就是（3（2-1-4）式中，，就是必要的基本起算数据的个数。的秩。这种情况下的列矩阵不满秩，即秩亏阵，秩亏数就是。这时，依据最小二乘准则就可以组成其法方程式：（3（2-1-5）式中，，并且。故此也是秩亏阵，其秩亏数即。由此可见，在进行经典的平差的时候必要的起算数据的个数就是不一类型控制网的秩亏数。因而，经常视平差问题的基准即为这种起算数据。把没有起算数据或起算数据不足，并且以待定点坐标为未知参数进行平差的控制网，由于缺少确定控制网坐标的基准而导致误差方程系数阵不列满秩，这种自由网被称作是秩亏自由网。设控制网中基准数据的个数是，则可得到：测边网、导线网、边角网：，一个点的坐标，一个方位；测角网：，一个方位，网中两个点的坐标或一个点的坐标，以及一条基线；高程网：，需要一个点的高程数据；GPS网：如果所采用的是GPS基线向量隐含的方位和尺度基准，那么，，要求已知一个点的三维坐标；三维控制网：，需三个已知定向角，一点的坐标和一个已知边。若一个控制网没有任何起算数据，并且平差的时候将控制网中的所有待定点都作为未知参数，列出误差方程式，这时未知参数的个数比间接平差相应的参数多了个。这样的情况下，误差方程式中的系数阵不列满秩，它的秩亏数为，像这种没有足够起算数据的平差问题就叫做秩亏自由网平差。在控制网秩亏的情况下，法方程有解但不唯一，换句话说，秩亏自由网平差和经典平差的本质区别就是，若只满足了最小二乘准则，仍然无法求出参数的唯一解，如果想求出唯一解，就需要增加新的约束条件。1.秩亏自由网平差的类型秩亏自由网平差就是在符合最小二乘和最小范数的条件下，求出参数一组最佳估值的平差方法。秩亏自由网平差可以根据最小范数的不同选择被分为以下三种类型。普通秩亏自由网平差它是在最小范数的限制条件下选择基准条件，来获得未知参数唯一估值的方法。拟稳平差把网中的未知参数分成两类，也就是设，它是在部分参数最小范数的要求下选择基准条件,来获得未知参数唯一估值的办法。加权秩亏自由网平差此平差是在加权最小范数的要求下选择基准条件，以此获得未知参数唯一估值的方法。可以看到，当，加权秩亏自由网平差就可以被认为是普通秩亏自由网平差，若取，加权秩亏自由网平差就退化成了拟稳平差。因此，加权秩亏自由网平差是秩亏自由网平差的一种普通形式，包括了普通秩亏自由网平差和拟稳平差。下面将分别介绍三种类型的秩亏自由网平差。普通秩亏自由网平差普通秩亏自由网平差的函数模型和平差原则为（3（2-1-6）其中，其秩亏数。且必要观测值个数为，未知参数个数为。依据最小二乘原理，可得法方程为：（3（2-1-7）其中。因秩亏，故秩亏。，秩亏数为，法方程的解不唯一。为了求出唯一解，需要再利用参数的约束条件。事实上，依据广义逆矩阵理论，相容方程组虽然具有无穷多组解，但有唯一的最小范数解，即：（3（2-1-8）其中，，叫做矩阵的最小范数广义逆。称为矩阵的逆。代入式（3（2-1-8）得：（3（2-1-9）上式即在广义逆矩阵理论的支撑下直接求解参数的最小唯一范数解的公式。鉴于广义逆进行计算时十分繁琐，下面会将公式继续简单化令：（3（2-1-10）得出（3（2-1-11）鉴于为行满秩矩阵，也就是，而，因此就是满秩方阵。以降阶法为依据获得矩阵广义逆的方法，也就是说如果有矩阵：此中，，有凯利逆，则有的阵：（3（2-1-12）根据上式得：（3（2-1-13）代入式（3（2-1-9）得：（3（2-1-14）或写成（3（2-1-15）未知参数的协因数阵为：（3（2-1-16）加权秩亏自由网平差加权秩亏自由网平差的函数模型和平差原则为：（3（2-1-17）其中，，其秩亏数。为必要观测值个数，为未知参数个数，利用最小二乘原理可以得到法方程：（3（2-1-18）其中，。由于秩亏，所以也秩亏。，秩亏数，法方程的解不唯一。为了求得唯一解，需要再利用参数约束条件。由：（3（2-1-19）可以得到（3（2-1-20）其中，，叫做矩阵的加权最小范数广义逆。称为矩阵的逆。代入上式得加权最小范数解为： （3（2-1-21）（3（2-1-22） 同样的，利用降阶法求矩阵广义逆的方法，就得到下面简化公式：（3（2-1-23）未知参数的协因数为：（3（2-1-24）拟稳平差在监测自由网中，假设其中有一部分点相较另外一部分点是相对稳定的。把网中全部点的坐标或高程当做未知数，那么这个坐标未知数就可以有两种划分，一种是稳定坐标未知数，一种是不稳定坐标未知数。假定其近似值为和，可以得到误差方程式：（3（2-1-25）式中，。依据最小二乘原理，由上式组成法方程为：（3（2-1-26）式中，且，非奇异，，奇异，由上式不能得到唯一解。消去不稳定未知参数，获得稳定未知参数的约化法方程，即： （3（2-1-27）式中，（3（2-1-28）由于经过约化后秩亏不会消除，仍为秩亏，式子（3（2-1-27）的解不唯一。为获得最优解，附加的最小范数条件。于是由（3（2-1-27）式可得：（3（2-1-29）式中，。将（3（2-1-29）代入（3（2-1-26）的第一式，得到：（3（2-1-30）式中，（3（2-1-31）由（3（2-1-29）和（3（2-1-30）可以得到：（3（2-1-32）残差以及单位权方差估计为：（3（2-1-33）（3（2-1-34）（4）秩亏自由网平差的附加条件法秩亏自由网平差就是在满足最小二乘和最小范数的条件下，求参数一组最佳估值的平差方法，实际上就是求相容方程组的最小范数解。附加条件法的基本思想是，由于网中没有初始数据，因此在平差中选择更多的未知参数，使得参数之间必须有附加条件，即在原平差模型中需要加入个未知参数间的限制条件方程，从而可以按附有条件的间接平差法进行求解。目前最大的问题是以什么方式如何推导出等价于最小范数的限制条件方程的具体形式。其次，给出限制条件方程推导平差计算公式，然后再证明，在给定的限制条件方程下所求得的解，就是相容方程组的最小范数解。2.加权秩亏自由网平差的解若用基准条件代替加权最小范数条件，则加权秩亏自由网平差的函数模型和平差原则为：（3（2-1-35）其中，，其秩亏数。这里的为必要观测值个数，为未知参数个数。为阶矩阵，叫做附加阵，并且，，表示附加的个基准条件与误差方程相互独立。在最小二乘原理的依据下，按照附有条件的间接平差法进行平差，可以获得法方程：（3（2-1-36）其中，，并且，。出于为秩亏阵，为解决秩亏这一问题，将上式的第一式左乘，得：（3（2-1-37）顾及到，故有。因此，有：（3（2-1-38）又因为，所以一定为零。由此可见，最小二乘准则和未知参数的附加基准条件并不相关，为不变量。附加基准条件仅代替加权最小范数条件以求得参数的唯一解。将式子（3（2-1-36）的第二个式子左乘和（3（2-1-36）的第一个式子相加，并且顾及到得：（3（2-1-39）由于满秩，正常逆存在。令：（3（2-1-40）则（3（2-1-41）由此可以看出，此方法被称为附加条件法实际上是通过增加未知参数间满足的个附加条件，而后依照附有条件的间接平差法而求得的。由此可见，其与经典的附有条件的间接平差法最本质的区别就是：附加的阵满足，并且。解为加权最小范数解下面将通过证明为加权的最小范数广义逆以此证明按附加条件法求得的解也是加权最小范数解。根据最小范数广义逆满足的矩阵方程可以得知，要证明为的加权最小范数广义逆，必须证明满足如下矩阵方程：；事实上，由，可得：；由，或者是（3（2-1-42）并且（3（2-1-43）由，得：（3（2-1-44）（3（2-1-45）显而易见，按附加条件法所得的解就是加权最小范数解。尽管加权最小范数广义逆不唯一，但是加权最小范数解唯一。因此下面求加权秩亏自由网平差参数解的模型等价，即： 等价于（3（2-1-46）换句话说，加权最小范数条件等价于下列基准条件：等价于（3（2-1-47）精度评定单位权方差估值的计算：（3（2-1-48）（3（2-1-49）按照协因数传播律，未知参数的协因数阵为：（3（2-1-50）由（3（2-1-43）知：（3（2-1-51）故：（3（2-1-52）又由式（3（2-1-42）得，，故：（3（2-1-53）计算时，标准化矩阵，使之满足：（3（2-1-54）又根据，得（3（2-1-55）从而： （3（2-1-56）2.2附加参数系统的平差在现实工程项目测量过程中，会有随机误差，系统误差和粗差的存在，测量实践表明，即使经过一系列的观测措施的运用以及各种预处理改正，残余的系统误差仍然是不可避免的。在经典的平差处理中往往会假设观测值中不包括系统误差，其实应消除或减弱这种残余的系统误差。因此就出现了一种平差的方法，就是通过在经典平差模型中附加参数系统参数对系统误差进行补偿，这种平差方法被称作附加系统参数的平差方法。附加系统参数的选择可以分为两种情况。一种是多项式型的附加参数，例如一般多项式，球谐函数中的系数，正形多项式作为附加参数。还有一种是考虑到系统误差特点的附加参数，例如卫星多普勒定位中的频偏、时延等未知数，三角高程网平差的折光未知数，测边网平差中的尺度比未知数。系统误差可以分为不变的系统误差和可变的系统误差。不变的系统误差就是误差符号和大小都固定不变的系统误差，可变的系统误差包括线性变化的系统误差和非线性变化的系统误差。非线性变化的系统误差又有多项式变化、周期性变化和按复杂规律变化的系统误差等。因不变的系统误差是线性变化的系统误差的特殊情况，非线性误差可用线性误差近似，故下只列出附加线性系统参数的平差方法。附加系统参数的平差模型经典高斯—马尔可夫模型为（3（2-2-1）第一种情况是，若假设观测值仅含偶然误差，则。第二种情况是，当观测值中除了偶然误差外，还包含有系统误差时，则。对于第二种情况，需要对高斯—马尔可夫模型进行补充。假设观测误差包含偶然误差和系统误差，即： （3（2-2-2）由于常系统误差是线性系统误差的特殊情况，而非线性又可用线性近似，因此，考虑附加线性系统参数，即设，于是有：（3（2-2-3）可以看出，。将（3（2-2-3）代入到（3（2-2-2），就得到附加系统参数的平差函数模型为：（3（2-2-4）（3（2-2-5）写出误差方程式：（3（2-2-6）式中，—系统参数；—系统参数的平差值；—系统误差的影响项，即对平差模型的补偿性；—系数矩阵，且，即为列满秩阵。下面详细介绍平差原理，并假设之间都是相互独立的。将（3（2-2-6）改写成矩阵形式：（3（2-2-7）根据最小二乘准则，得到其法方程的系数和常数项为：因此，法方程为：（3（2-2-8）记：。故上式可以简写成：（3（2-2-9）根据分块矩阵求逆公式得：（3（2-2-10）式中，。如果平差模型中不含有系统误差，即，则有：（3（2-2-11）故，由（3（2-2-10）可以得到：（3（2-2-12）（3（2-2-13）由间接平差可以得知，平差参数估值的协因数阵就是法方程系数阵的逆阵，因此，（3（2-1-14）即和的协因数阵为：（3（2-1-15）（3（2-1-16）其单位权中误差为：（3（2-1-17）利用附加系统参数的平差方法可以有效地补偿观测数据之中的系统误差，但是在实际选择和处理附加参数时，还必须考虑实际会遇到的问题，为了尽可能地补偿系统误差，人们总是希望用全面且包含了大量附加参数的系统误差模型来扩展平差的函数模型。但是结果往往会造成参数过度化，使得附加参数间或附加参数与基本参数之间存在近似的线性关系，使得法方程系数系数阵病态化，换句话说，系数阵接近奇异，最终的结果是导致最小二乘平差求得的参数估值变坏。因此，必须采取有效措施来避免这种情况的发生，附加系统参数的引入，改变了原平差模型。需要对附加系统参数的显著性进行检验以此保证平差模型的正确性。若系统参数不存在，或者是系统参数存在但是和列入模型的项不相符，但是仍然采用模型（3（2-2-6）式进行平差，就必然会影响平差结果的正确性。所以附加系统参数的平差必须对模型的正确性、加入项的显著性等进行检验。总之，平差前对附加系统参数做统计检验是必要的，以便采用正确的平差模型，获得可信的平差结果。2.3最小二乘滤波在经典最小二乘平差中，通常假定平差参数是非随机变量；或者是没有先验统计性质的随机变量，在参数估计时无法顾及其随机性；或者是具有先验统计性质，但是在参数估值的时候不考虑随机性的随机变量，按照最小二乘原理求定最佳估值。但是在有些实际问题中，首先需要求定最佳估值的参数是具有先验统计性质的随机参数，并且在估计这些参数的时候必须考虑到这些统计性质，以及模型中随机参数与非随机参数同时存在，且需要顾及随机参数先验统计性的参数估计问题。这类附有随机参数的平差问题即滤波、推估和配置。在测量平差中，滤波就是利用含有误差的观测值，求定参数最佳估值的方法。滤波的函数模型由以下观测方程表示：（3（2-3-1）式中，—已知系数矩阵；—观测值向量；—观测误差向量；—随机参数；其中，分两种情况，一种是与观测值之间有函数关系的已测点参数，称为滤波信号，用向量表示，另一种是和观测值之间没有函数关系的未测点参数，称为推估信号，用向量表示，并且和统计相关，即：通过滤波可以求得和的最佳估值和。通常将求定滤波信号的最佳估值的过程称为滤波，将求定推估信号的最佳估值的过程称为推估。推估还可以分为两种，一种是内插或平滑，也就是说，对于静态系统，如果在的范围内，就是平滑或内插；第二种就是外推或预报。如果在的范围以外，就是外推或预报。配置也被称作拟合推估，最初是组合各种资料研究地球形状和重力场的一种数学方法。配置的普遍形式是在其函数模型中，既有非随机参数部分（倾向参数），又有随机参数部分（滤波信号和推估信号），也就是说，在拟合推估中既包含求信号估值的任务，又包含求倾向参数估值的任务。将高斯—马尔可夫模型加以扩展，得到最小二乘配置的函数模型为：（3（2-3-2）式中，为观测向量，为非随机参数，为随机参数。同上，又分为两种情况，第一种是与观测值之间有函数关系的已测点参数，第二种是与观测值之间没有函数关系的未测点参数。很明显，最小二乘配置的函数模型是一个既包含非随机参数，又含有随机参数的模型。当或时，模型就变为，即高斯—马尔可夫模型。当或时，模型就变为，即无系统参数的配置模型，也称作滤波和推估模型。当时，模型变为，称为滤波模型。1.最小二乘滤波与推估函数模型和随机模型滤波的函数模型为：（3（2-3-3）式中，—观测向量；—随机参数其中，（3（2-3-4）由和的先验信息构成的随机模型为：（3（2-3-5）其中，若，即表示观测误差与随机参数相互独立。最小二乘滤波的误差方程按照广义最小二乘原理，将作为虚拟观测值，其权为，方差为。令，其权则为（3（2-3-6）则虚拟观测方程为：（3（2-3-7）式中，为对应的观测误差。于是与对应的误差方程为： （3（2-3-8）这时，随机参数已经转化为非随机参数，可以按非随机参数的处理方式进行处理。将滤波的函数模型也改写为误差方程的形式：（3（2-3-9）式中，。由上述两式即构成最小二乘滤波的总误差方差：（3（2-3-10）改写成矩阵形式为：（3（2-3-11）即：（3（2-3-12）式中，。则将最小二乘滤波问题转化成了一般间接平差问题。最小二乘滤波和推估的计算公式分两种情况来分别说明最小二乘滤波和推估。1.根据最小二乘原理，即：（3（2-3-13）利用间接平差计算公式，可以得到法方程：（3（2-3-14）将代入到上式，得到：（3（2-3-15）解得：（3（2-3-16）取（其中，取），则上式就写成了：（3（2-3-17）根据矩阵反演公式：（3（2-3-18）由式（3（2-3-17）可得：（3（2-3-19）由于：（3（2-3-20）所以：（3（2-3-21）将上式展开，可以得到：（3（2-3-22）由间接平差可以知道，参数的权逆阵就是法方程系数阵的逆阵，因为，当取，有：（3（2-3-23）由矩阵反演公式（3（2-3-18）可得:（3（2-2-24）将（3（2-2-20）代入（3（2-2-24），可以得到：（3（2-2-25）即：（3（2-2-26）式子（3（2-3-22）和式子（3（2-2-26）就是滤波和推估求和的估值和以及其方差、协方差的公式。2.取，根据最小二乘原理，即：（3（2-3-27）根据间接平差计算公式，得到法方程为：（3（2-3-28）式中，。类似第一种情况，即可得出下面结论：（3（2-3-28）2.4序贯平差在工程实际施工场地上，随着工程项目的不断开展，需要不断地布设高程点来满足测量工作的需要。传统的方法是每次布设完毕后将全网中的所有点进行整体平差，这样随着水准网的不断扩大，观测数据越来越多，就使得计算时的方程矩阵的阶数越来越大，从而给计算造成了较大的繁琐。序贯平差作为一种新的平差方法，可以解决这类问题，尤其对于超大型的网，从而使得网不断扩展的同时而计算可以分步进行，但是不会出现整体平差时带来的观测数据超多、计算超麻烦的情况，降低法方程阶数后，就减轻了计算强度。序贯平差也称为逐次相关间接平差，或者也叫做静态卡尔曼滤波，是近代出现的一种最小二乘分解法。它是将观测值分成两组或多组，按组的顺序分别做相关间接平差，不用考虑前一个阶段的观测值，但是利用前期的平差结果，达到与两期或多期网一起整体平差同样的结果，由于序贯平差有一套规律性很强的递推公式，无论观测数据如何增加，无需解算高阶逆矩阵，就能够不断对原平差参数进行改进。同时，由于法方程阶数不高，因此解算时数值稳定性较好，应用范围就相当广泛。平差原理将观测值分为两组，记作和，它们的权阵分别记为，假设这两组观测值互不相关，即有（3（2-4-1）其中，是必要观测数。当参数之间不存在约束条件时，其误差方程为：（3（2-4-2）（3（2-4-3）式子中，。将（3（2-4-3）式子单独平差，得到：（3（2-4-4）式子中为（也就是）的协因数阵，故有：（3（2-4-5）或者是：（3（2-4-6）表示的是由第一组观测值平差所得的值。将（3（2-4-2）和（3（2-4-3）式联合解算，即由两组观测值作为整体平差，可以组成法方程为：（3（2-4-7）其解为：它是两组观测值整体平差的结果，按间接平差可知，其法方程系数的逆阵就是的协因数阵，故有：（3（2-4-8）故上式可写成：（3（2-4-9）再由（3（2-4-4）式和（3（2-4-8）式可得：（3（2-4-10）（3（2-4-11）则（3（2-4-9）式可以写为：（3（2-4-12）令（3（2-4-13）（3（2-4-14）则上式可以简写为：（3（2-4-15）将（3（2-4-11）式两边左乘，又根据式（3（2-4-13）可得：（3（2-4-16）再对两边右乘，得（3（2-4-17）由矩阵反演公式得知（3（2-4-18）相较（3（2-4-17）和（3（2-4-18）式得知：（3（2-4-19）阵称为卡尔曼滤波增益矩阵或称为序贯平差的增益矩阵。将（3（2-4-19）式代入（3（2-4-15）式就得到序贯平差的递推计算式。当时，（3（2-4-19）式中的和由（3（2-4-14）式计算的都是纯量，前者的逆也是一个数，在这种情况下，序贯平差的计算就非常简单。这是应用序贯平差的主要场合，也就是说第二组总是假定，逐次递推，求出整体平差的最后结果。精度评定单位权中误差为（3（2-4-20）当两组观测值整体平差时，有（3（2-4-21）由于：式中（3（2-4-22）即为仅用第一组观测值单独平差时所算得的改正数。顾及（3（2-4-6）式以及单独平差第一组观测值时应满足，可得：（3（2-4-23）故：（3（2-4-24）又因为：（3（2-4-25）由于（3（2-4-14）式，可将（3（2-4-25）式写为：（3（2-4-26）将（3（2-4-19）式代入到上式，得：（3（2-4-27）故：（3（2-4-28）由于：故：将（3（2-4-29）式代入（3（2-4-24）式可得的递推公式：（3（2-4-29）参数的协因数阵递推公式为：（3（2-4-30）由此，可以得到序贯平差的一组递推公式：（3（2-4-31）其中：（3（2-4-32）43序贯平差算法程序设计与实现43.1水准网概况在下图所示的水准网中，A、B点为已知点，第一次同精度独立观测第二次同精度独立观测观测值为分别为。（图34-1）43.2整体法平差与序贯法平差1.整体法平差由于必要观测数选取C、D、E三点高程为参数，取未知参数近似值，依据图形列出平差值条件方程式，计算误差方程如下：改写成矩阵形式，即的形式：，组成法方程为：求出参数，计算参数的平差值，由误差方程式计算V,求出观测量平差值：计算单位权中误差为：序贯法平差因必要观测数选取C、D、E三点高程为参数，取未参数近似值，列出第一期误差方程：同样地，改写成矩阵形式为：组成法方程式：即可解出参数的第一次改正数、第一次平差值以及权阵和第一期观测值的第一次改正数：列立第二期误差方程，可以用第一期平差后的参数平差值直接列立，这时的误差方程常数项就是，即：改写成矩阵形式：考虑第一次的平差结果，组成法方程：求解参数的第二次改正数以及平差值计算第二期观测值的改正数：计算单位权中误差43.3程序实现1.数据格式以本文的算例为例，在数据文件中，第一行分别为观测值总数和参数总数；第二行开始至第七行为误差方程式，前三列为系数阵，最后一列为自由项；第八行至第十四行是观测值的权矩阵，采用下三角存储方式；第十五行至第十七行是参数的先验值；第十八行至第二十行是参数的先验权逆阵，同样地采用下三角存储方式。具体数据文件见下图：具体数据见下图：（图34-1数据文件图3）2.平差公式以及具体计算步骤假设误差方程式为（3-3-1）式中：为n维残差向量；为阶系数矩阵；为t维参数向量；为n维观测值向量，其权矩阵为对称正定矩阵；参数的先验值为，先验权逆阵为。对于序贯平差而言，为已知矩阵（或向量），是待求量。根据最小二乘原理，有下列计算公式： （3-3-2）（3-3-3）（3-3-4）（3-3-5）（3-3-6）（3-3-7）（3-3-8）其中，，可以证明,则（3-3-9）下面列出具体计算步骤（1）由（3-2-2）计算预估残余向量；（2）由（3-2-3）计算预估残差的权逆阵；（3（3）由（3-2-4）计算矩阵M；（4）由（3-2-5）计算增益矩阵J；（5）由（3-2-6）计算参数平差值的权逆阵；（6）由（3-2-7）计算参数平差值；（7）由（3-2-9）计算单位权中误差。3.功能模块实现（1）（1）对称矩阵的下标计算由于此程序中涉及到对称矩阵的计算，例如观测值的权矩阵、权逆阵、法方程系数矩阵等都是对称矩阵，为了节省存储空间和避免重复计算，故仅存储下三角矩阵（包含对角线元素）的存储方案存储对称矩阵，也就是说仅将对称矩阵的下三角元素按顺序存到数组中。对称矩阵的下标计算，就是当对称矩阵仅存下三角矩阵时根据矩阵元素的行号和列号确定相应元素在数组中的存储位置。其函数的源代码如下：intij(inti,intj){return(i>=j)?i*(i+1)/2+j:j*(j+1)/2+i;}（2）对称正定矩阵求逆boolinverse(doublea[],intn)a——函数调用前为待求逆矩阵的元素，调用结束后为逆矩阵元素；n——矩阵的阶数；函数返回值——若计算成功返回true；若计算失败返回false。{double*a0=newdouble[n];for(intk=0;k<n;k++){doublea00=a[0];if(a00+1.0==1.0){delete[]a0;returnfalse;}for(inti=1;i<n;i++){doubleai0=a[i*(i+1)/2];if(i<=n-k-1)a0[i]=-ai0/a00;elsea0[i]=ai0/a00;for(intj=1;j<=i;j++){a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+ai0*a0[j];}}for(inti=1;i<n;i++){a[(n-1)*n/2+i-1]=a0[i];}a[n*(n+1)/2-1]=1.0/a00;}delete[]a0;returntrue;}（3）权逆阵传播计算设B为r×n矩阵，Q为n×n对称矩阵，以N表示它们的乘积，即,则N的第i行，第j列元素的计算公式就为voidCalculate_BQBT(doubleB[],doubleQ[],intr,intn,doubleN[])B——B矩阵，数组长度为r×n；Q——Q矩阵，数组长度为n×（n+1）/2；N——N矩阵，数组长度为r×（r+1）/2；r——B矩阵的行数；n——B矩阵的列数。{for(inti=0;i<r;i++)for(intj=0;j<=i;j++){doublenij=0.0;for(intk=0;k<n;k++)for(ints=0;s<n;s++)nij+=B[i*n+k]*Q[ij(k,s)]*B[j*n+s];N[ij(i,j)]+=nij;}}（4）序贯平差函数原型doubleSuccessive_Adjustment(intn,intt,doubleA[],doubleL[],doubleP[],doubleX[],doubleQX[],doubleV[])n——观测值个数；t——参数个数；A——误差方程系数矩阵，数组长度为n×t，数组内容已知；L——观测值向量，数组长度为n，数组内容已知；P——观测值权矩阵的下三角元素，数组长度为n×（n+1）/2，数组内容已知；X——参数向量，函数调用之前为参数先验值，函数调用后为参数平差值，数组长度为t；QX——参数权逆阵，函数调用之前为参数先验权逆阵，函数调用之后为参数平差值，仅保存下三角矩阵元素，数组长度为t×（t+1）/2；V——残差向量，数组长度为n，数组内容待计算；函数返回值——若计算成功，返回单位权中误差μ；若计算失败，返回-1.0。函数源代码就不在此处列出。4.计算结果分析经过程序运算，得到本文算例的计算结果如下图所示：图三计算结果图（1）整体法平差与序贯法平差得到的计算结果对比分析如下表：表3-1整体法与序贯法结果对比整体法平差87.363m88.046m86.167m8.6mm序贯法平差87.363m88.046m86.167m8