秩亏网平差若干计算方法.docx

秩亏网平差若干计算方法1.概述在测量平差中，控制网中除了必要起算数据外还有多余起算数据的是附合网，仅有必要起算数据的是自由网，这两种控制网在间接平差时误差方程系数矩阵B都是满秩的，由此得到的法方程系数阵N=BTPB也是满秩的，即法方程有唯一解。这是经典平差的范畴。自由网中有一种具有特殊用途的控制网，就是秩亏自由网，这种自由网没有起始数据参与平差并且以待定点的坐标为待定参数。此时的误差方程的系数阵B是列亏阵，由此所得的法方程系数阵N=BTPB也是秩亏阵。一般设网中全部的待定坐标个数为u，必要观测数为t，全部观测数为n，B为nu阶矩阵，相应的法方程系数阵N是uu阶矩阵， RB=RN=tu，秩亏数都为d=u-t，所以法方程有无穷组解。这里产生秩亏的原因是控制网中没有起算数据，所以d就是网中必要的起算数据个数。对于水准网，必要起算数据是一个点的高程，故d=1；对于测角网，必要起算数据是两个点的坐标，故d=4；对于测边网或是边角网，必要起算数据是一个点的坐标和一条边的方位，故d=3。2.秩亏网平差模型以间接平差为例，令u个坐标参数的平差值为X，观测向量为L，则秩亏网的误差方程为：V=Bx-l （1）式中，RB=tu，d=u-t，X=X0+x，l=L-L0随机模型是：D=2Q=2P-1 （2）根据最小二乘原理，在VTPV=min下，可组成发方程如下：BTPBx-BTPl=0 （3）若是按照直接解法用如下的方程组来解求x的解：V=Bx-lBTPBx-BTPl=0VTPV=min （a）容易得到|BTPB|=0，即该方程组有解但不唯一，虽然满足最小二乘准则，但有无穷多组x的解，无法求得唯一的x，因为参数x必须在一定的坐标基准下才能唯一确定。为了得到x的唯一解，增加d个坐标基准约束条件，即：ST x=0 (4)在限制条件VTPV+2KT(ST x)=min下，得到法方程如下：BTPBx+SK=BTPlST x=0 （5）由此可以根据下面的方程组解得x的唯一解：V=Bx-lBTPBx+SK=BTPlST x=0VTPV+2KT(ST x)=min （b）由上述方程组（b），可以得到：x=（BTPB+SST ）-1BTPl Q=（BTPB+SST ）-1BTPB（BTPB+SST ）-1 (7)3.矩阵分解应用于秩亏网平差3.1 奇异值分解用于秩亏网平差可以看出，上面提到的这种计算秩亏网平差的方式很复杂，现在我们不妨把秩亏自由网平差看成在满足最小二乘VTPV=min和最小范数xTx=min的条件下，求参数一组最佳估值的平差方法，也就是通过对如下的方程组来解求x的唯一解：V=Bx-lBTPBx-BTPl=0VTPV=minxTx=min (c)这是个复杂的方程组，如果按照正常求解的方法是很困难的，下面我们把矩阵的奇异值分解融合进来。矩阵的奇异值分解是线性代数中一种重要的矩阵分解，在最优化问题、特征值问题和广义逆问题及统计学问题中都有重要的应用。对秩为r的mn阶矩阵A进行奇异值分解的步骤是：1）求得AHA的特征值1，2,n,及对应的特征向量并正交单位化，得矩阵V，使得VHAHAV=M2000,M=diag(1，2,n)；2）将V的前r列作为V1，令U1=AV1M-1，再扩张成m阶的矩阵U；3）那么A=UM000VH。根据上述矩阵奇异分解的原理求得矩阵B的奇异值分解：B=UM000VH,在此基础上令矩阵G=VM-1000UH。通过矩阵理论的学习我们知道，可以通过如下的方式来验证G就是B的广义逆：（1）BGB= UM000VH VM-1000UHUM000VH=UM000VH=B(2) GBG=VM-1000UH UM000VHVM-1000UH=VM-1000UH=G(3) (BG)H=(UM000VHVM-1000UH)H=BG(4) (GB)H=(VM-1000UHUM000VH)H=GB对于不相容方程组Bx=b,使得x=Gb为极小范数最小二乘的充要条件是G为B的广义逆，这就说明G是满足该方程式的极小范数最小二乘解，也就是说，我们得到未知参数的估值x=Gl=VM-1000UHl。通过这种计算方法，我们求解方程组（c）就简单多了。3.2 满秩分解用于秩亏网平差满秩分解，是指把秩为r的mn阶矩阵A分解成A=FG，其中F是秩为r的mr阶矩阵，G是秩为r的rn阶矩阵。满秩矩阵的解求利用Hermite行标准形来完成。把矩阵A经过处等行变换成为Hermite行标准形B，B的j1,j2,jr列为单位矩阵Im的前r列，令A的第j1,j2,jr列为矩阵F，B的前r行为矩阵G，则有A=FG。在这里，秩亏网平差的数学模型还是前面提到的方程组（c），根据广义逆理论，相容方程组Nx-W=0虽有无穷多组解，但它有唯一的最小范数解，即：x=Nm-W （8）式中，Nm-=NT（NNT）-，代入上述（8）中，有：x=NT（NNT）-W （9）这就是根据广义逆理论直接解求参数的唯一最小范数解的公式。由于广义逆计算较为复杂，且在计算机语言编程中不易实现，因此可以应用满秩分解来进行改进。同前面提到的一样，首先根据前面的最大秩分解原理，得到NNT的满秩分解为NNT=FG。令C=GR-1FL-1，式中，FL-1=（FTF）-1FT，GR-1=GT（GTG）-1。另外，只要NNT和C满足以下四个条件时，即：（NNT）C（NNT）=NNT、C（NNT）C=C、(NNTC）T=NNTC和(CNNT）T=CNNT，就说明C是NNT的广义逆。则（NNT）-=GR-1FL-1，从而得到未知参数估值：x=NTGR-1FL-1W。经过这样的改进，秩亏网平差的算法更利于计算机语言编程。 4. 遗传算法在秩亏网平差中的应用为了确定秩亏网中未知数估值的唯一性，需要加入基准条件，并与法方程联立来进行参数估计，这会涉及大量又复杂的矩阵运算，遗传算法的引入可以使计算简单而有效。遗传算法模仿的机制是一切生命与智能的产生与进化过程，具有智能性、良好的并行性、可操作性和全局优化性等特点。秩亏自由网平差的法方程组用下式表示：NX=W （10）令FXNX-W，其中FX=（f1X，f2X，ftX）T，为了方便用遗传算法来解求奇异方程组（10），将之转换成一个多维优化问题来解决。令目标函数ObjFX=FT（X）FX，则无约束多维化问题为：Min(ObjFX) （11）方程组（10）的最小二乘解等价于上面多维无约束优化问题（11）的全局近似最优解。在此过程中涉及的问题有：1）参数编码方案。对于问题（11），把参数向量的每个分量先进行二进制编码得到各个子串，再把各个子串参数的顺序排列在一起形成一个解个体。2）适应度函数的设计。对于问题（11），采用下式计算个体染色体的适应值：即适应度函数：FitnessX=exp(-ObjFX) (12)该适应度函数变化趋势平缓，有利于趋近全局最优解。3）遗传算法基本参数选择。根据众多的研究成果并结合所研究的问题，取群体规模L为40-100，交叉概率Pc 为0.4-0.8，变异概率Pm 为0.001-0.01。4）遗传操作算子设计。按适应度把个体从大到小排列，然后把序号在前的e个个体复制两份，覆盖序号在后面e个个体，而序号在中间的个体复制一份，其中e的大小视群体规模而定，一般为1-5。另外，交叉算子采用单点交叉，变异算子采用基本位变异。5）迭代终止条件。采用目标函数值ObjFX3(m+n)。令m个测站的三维坐标为（Xi，Yi，Zi）， n个点的坐标为（xj，yj，zj），对于任意测站i和点j，有误差方程: Vij=(Xi-xj)2+(Yi-yj)2+(Zi-zj)2 (13)对所有的测站i和点j，有误差方程组如下： V=AX-l （14）其中，A为mn3(m+n)阶矩阵，未知参数t=3(m+n)，观测量个数为mn。由于mn3(m+n)，所以必要观测数据足够，但由于没有起始数据，只能应用秩亏网平差来解决该问题。普通秩亏网的附加约束法的解算模型为： V=Bx-l STx=0 （15）利用上述涉及到的求解秩亏网的各种方法都可以解出方程组（15）中的未知参数 ，即：x=BTPB+SST-1BTPl 。可以看出，秩亏网的应用，很简单的解