最小二乘平差之附有参数的条件平差

在正式开始正题之前，先引述一下模式可分性的Cover定理：
假设空间不是稠密分布的，将复杂的模式分类问题非线性地投射到高维空间将比投射到低维空间更可能是线性可分的。
Cover定理在机器学习中较常被提及，比如：支持向量机。但由于平差的对象空间一般是稠密的，此定理似乎不太适合于这里，但原理是相通的，也就是低维空间的问题一一投射到高维空间往往能得到更好的解决方案。从数学角度上讲，上面的增加空间的维度和本文要讲的增加拟合的参数是等价的。
附有参数的条件平差也可看作是一种间接平差和条件平差的混合模型。增加条件方程是“降维”，它可以把复杂的问题变得简单；添加参数是“升维”，它可以把问题变得理论上可分的概率更大。两者并不矛盾，而且相辅相成。我们中学初等几何中常用的添加辅助线、辅助面的方法其实就是附加参数以增加条件方程。
为了说明条件平差附有参数的情况，在这里先举一个例子，这是书上的一个例子，很多同学对它不陌生。

A、B为已知点，AC为已知边。在网中观测了9个角度，即n=9，为了确定C、D、E三点的坐标，其必要观测数t=5，故r=n−t=4，因此，必须列出4个条件方程。其中3个图形条件是很容易列出的，但要列出第4个条件就比较困难了。为了解决这一问题，选取非观测量∠ABD作为参数X^，这样就要列出c=r+u=4+1=5个条件方程。此时除了3个图形条件外，还可列出1个极条件和1个固定边条件。若以A点为极，则极条件为：

sin(L^5+L^7)sinX^sinL^6sin(L^9−X^)sin(L^6+L^8)sin(L^5)=1

固定边条件（由AC推算到AB）为：

SAB=SACsinL^2sin(L^6+L^8)sinL^3sinX^

从这5个条件方程进行平差，该问题就得到了解决。
一般而言，在某一平差问题中，观测值个数为n，必要观测数为t，多余观测数r=n−t，再增选u个独立参数，0<u<t，则总共应列出c=r+u个条件方程，一般形式为

Fc1(L~,X~)=0

如果条件方程是线性的，其形式为

AL~+BX~+A0=0

或者

AV+BX~+W=0

式中V为观测值L的改正数，x^为参数近似值X0的改正数，即

L^=L+V,X^=X0+x^

这里c=r+u,c<n,u<t,系数矩阵的秩分别为

R(A)=cR(B)=u

即A为行满秩矩阵，B为列满秩矩阵。式中

W=AL+BX0+A0

此平差问题，由于增选了u个参数，条件方程总数由r个增加到c=r+u个，平差自由度，即多余观测数不变，为r。
为了求出原问题的解，依附最小二乘及拉格朗日乘子，组成函数

Φ=VTPV−2KT(AV+Bx^+W)

为求Φ的极小值，将其分别对V和x^求一阶导数并令其为零，则有

∂Φ/∂V=2VTP−2KTA=0

∂Φ/∂x^=−2KTB=0

由两式转置得

PV−ATK=0

BTK=0

于是三个方程组成的方程组

AV+Bx^+W=0V=P−1ATK=QATKBTK=0⎫⎭⎬⎪⎪

称为附有参数的条件平差的基础方程。解算此基础方程，通常是将其中的改正数方程代入条件方程，得到一组包含K和x^的对称线性方程组，即

AQATK+Bx^+W=0BTK=0}

这里令NAA=AQAT，上式可写为

NAAK+Bx^+W=0BTK=0}

上式称为附有参数的条件平差的法方程。由于NAA是对称满秩矩阵，故用其左乘上式法方程的第一式，得

K=−N−1AA(Bx^+W)

又以BTN−1AA左乘法方程第一式，并与第二式相减，得

BTN−1AABx^+BTN−1AAW=0

设

NBB=BTN−1AAB

则有

NBBx^+BTN−1AAW=0

由于NBB可逆，所以

x^=−N−1BBBTN−1AAW

这样，V可由下式给出

V=−QATN−1AA(Bx^+W)

在算出x^后，直接计算改正数V。最后计算平差值。

L^=L+V

下面是具体的代码实现，其中基本的矩阵运算没有在下面给出，在矩阵算法相关代码，如有需要可以下载。

// <summary>   /// 附有参数的条件平差/// </summary>     /// <param name="correction">返回的改正数</param> /// <param name="MatrixA"></param>/// <param name="MatrixB"></param>/// <param name="w"></param>/// <param name="C"></param>/// <param name="N"></param>/// <param name="U"></param>//////////////////////////////////////////////////////////////////////////template<class T>void GetConditionCorrectionWithParameters(T correction[],const T matrixA[],const T matrixB[],const T w[],int C,int N,int U){    T *transposedA=new T[N*C];    T *transposedB=new T[U*C];    T *Naa=new T[C*C];    T *inverseForNaa=new T[C*C];    T *Nbb=new T[U*U];    T *inverseForNbb=new T[U*U];    T *x=new T[U];    MatrixTranspose(matrixA,transposedA,C,N);    MultMatrix(matrixA,transposedA,Naa,C,N,C);    MatrixAnti(Naa,inverseForNaa,C);    MatrixTranspose(matrixB,transposedB,C,U);    T *transit1=new T[U*C];    MultMatrix(transposedB,inverseForNaa,transit1,U,C,C);    MultMatrix(transit1,matrixB,Nbb,U,C,U);    MatrixAnti(Nbb,inverseForNbb,U);    T *transit2=new T[U*C];    MultMatrix(inverseForNbb,transposedB,transit2,U,U,C);    T *transit3=new T[U*C];    MultMatrix(transit2,inverseForNaa,transit3,U,C,C);    MultMatrix(transit3,w,x,U,C,1);    MatrixMinus(x,x,U,1);    T *transit4=new T[N*C];    MultMatrix(transposedA,inverseForNaa,transit4,N,C,C);    T *transit5=new T[C];    MultMatrix(matrixB,x,transit5,C,U,1);    T *transit6=new T[C];    MatrixAdd(transit5,w,transit6,C,1);    MultMatrix(transit4,transit6,correction,N,C,1);    delete [] transposedA;    delete [] transposedB;    delete [] Naa;    delete [] inverseForNaa;    delete [] Nbb;    delete [] inverseForNbb;    delete [] x;    delete [] transit1;    delete [] transit2;    delete [] transit3;    delete [] transit4;    delete [] transit5;    delete [] transit6;}

最后说一下精度评定，单位权方差可由下式估值

σ^20=VTPVr=VTPVc−u

其中，R(A)=c， R(B)=u。
平差值函数的协因数的计算公式为

Qφ^φ^=FTQL^L^F+FTQL^X^Fx+FTxQX^L^F+FTxQX^X^Fx

QL^L^,QL^X^,QX^L^,QX^X^可在参考文献上查表获得，这里不详述。
最后，φ^的中误差为

σ^φ^=σ0^Qφ^φ^−−−−√

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

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