最小二乘平差之间接平差

今年3月份的时候我在公司写一个算法接触到测量平差理论，感慨于它的算法的优美，利用最小二乘而又独成体系，在测量拟合中发挥着巨大的作用。但网上关于它的描述甚少，更不用说好的实现代码，我将写两篇文章来说一下间接平差和附加限制条件的间接平差，这一篇先说间接平差。
在测量拟合中，由于观测结果不可避免地存在着误差，因此，如何处理带有误差的观测值，找出待求量的最佳估值，是最小二乘平差所研究的内容。在一个平差问题中，当所选的独立参数X^的个数等于必要观测数t时，可将每个观测值表达成这t个参数的函数，组成观测方程，这种以观测方程为函数模型的平差方法，就是间接平差。用数学语言描述就是：通过在映射B的像上的平差间接地反映原像上的平差过程的算法。
一般而言，如果某平差问题有n个观测值，t个必要观测值，选择t个独立量作为平差参数X^t1，则每个观测值必定可以表达成这t个参数的函数，即有：

L^n1=F(X^)

如果表达式是线性的，一般为：

L^n1=BntX^t1+dn1

其中自由度是：r=n−t。
平差时，一般对参数X^都要取近似值X0,令

X^=X0+x^

代入上式，并令

l=L−(BX0+d)=L−L0

L0=BX0+d为观测值的近似值，所以l是观测值与其近似值之差，由此可得误差方程

V=Bx^−l

平差准则为

VTPV=min

以上是基本原理，下面就逐步推到得到平差公式。
设有n个观测值方程为：

L1+v1=a1X1^+b1X2^+⋯+t1Xt^+d1L2+v2=a2X1^+b2X2^+⋯+t2Xt^+d2⋯Ln+vn=anX1^+bnX2^+⋯+tnXt^+dn

令

Xj^=X0j+xj^(j=1,2,⋯,t)

li=Li−(aiX1^+biX2^+⋯+tiXt^+di)(i=1,2,⋯,n)

可得最终的误差方程为

vi=aix1^+bix2^+⋯+tixt^−li(i=1,2,⋯,n)

其中：

Bnt=⎡⎣⎢⎢⎢⎢a1a2⋮anb1b2⋮bn⋯⋯⋯t1t2⋮tn⎤⎦⎥⎥⎥⎥

Vn1=[v1v2⋯vn]T

xt1^=[x^1x^2⋯x^t]T

ln1=[l1l2⋯ln]T

按最小二乘原理，可得：

∂VTPV∂x^=2VTP∂V∂x^=VTPB=0

上述方程与V=Bx^−l联立，解得：

BTPBx^−BTPl=0

令

NBB=BTPB,W=BTPl

上式可简写为：

NBBx^−W=0

式中系数阵NBB为满秩，x^有唯一解，上式称为间接平差的法方程。解之得：

x^=N−1BBW

从而平差结果为：

X^=X0+x^

其中，P是观测数据的协因数阵的逆矩阵，一般可认为是单位矩阵。
下面是具体的代码实现，其中基本的矩阵运算没有在下面给出，在矩阵算法相关代码，如有需要可以下载。

template<class T>void GetWu1(T Wu1[],const T *matrixB,const T *l,int N,int U){    T *transposedB=new T[N*U];    MatrixTranspose(matrixB,transposedB,N,U);    MultMatrix(transposedB,l,Wu1,U,N,1);    delete [] transposedB;}template<class T>void GetNBB(T nbb[],const T *matrixB,int N,int U){    T *transposedB=new T[N*U];    MatrixTranspose(matrixB,transposedB,N,U);    MultMatrix(transposedB,matrixB,nbb,U,N,U);    delete [] transposedB;}/// <summary>   /// 间接平差/// </summary>     /// <param name="correction">返回的改正数</param> /// <param name="matrixB"></param>/// <param name="l"></param>/// <param name="N"></param>/// <param name="U"></param>//////////////////////////////////////////////////////////////////////////template<class T>void GetCorrection(T correction[],const T *matrixB,const T *l,int N,int U){    T *nbb=new T[U*U];    T *inverseForNbb=new T[U*U];    T *Wu1=new T[U];    GetNBB(nbb,matrixB,N,U);    MatrixAnti(nbb,inverseForNbb,U);    GetWu1(Wu1,matrixB,l,N,U);    MultMatrix(inverseForNbb,Wu1,correction,U,U,1);    delete [] nbb;    delete [] inverseForNbb;    delete [] Wu1;}

最后再说一下精度评定，首先通过下面的公式得到单位权中误差：

σ^0=VTPVn−t−−−−−−√

其中n为方程数，t为拟合向量的维数。平差值的协因数阵根据下面公式求得：

Qx^x^=N−1BB

假定间接平差问题中有t个参数，设参数的函数为：

φ^=ϕ(X1^,X2^,⋯,Xt^)

为求函数ϕ的中误差，先对上式全微分得权函数式为：

dφ^=f1dX1^+f2dX2^+⋯+ftdXt^

令FT=[f1f2⋯ft]，则函数φ^的协因数阵为：

Qφ^=FTQx^x^F=FTN−1BBF

如果要得到x^的中误差，可令FT=[11⋯1]，则

Qx^=FTQx^x^F=FTN−1BBF

最后由下面公式可算得平差值函数的中误差：

σ^φ^=σ^0×

Qφ^φ^−−−−√

参考文献：
[1] 误差理论与测量平差基础 武汉大学测绘学院测量平差学科组编著

转载请注明出处：http://blog.csdn.net/fourierFeng/article/details/47272167