测量平差之附有限制条件的条件平差（概括平差模型）

接下来要写的可看成是对前面讲过的四种经典平差方法的概括模型。
从四种基本平差算法的函数模型来看，主要包括如下两种类型的条件方程：$$\{\begin{array}{c}F(\hat{L},\hat{X})=0\\ \varphi (\hat{X})=0\end{array}$$
第一种方程中同时含有观测值和未知参数，称为一般条件方程；第二种方程则只含有未知参数而无观测值，称为限制条件方程。
一般而言，对于任何一个平差问题，设想观测值个数为$n$，必要观测数为$t$,则多余观测数$r=n+t$。若增选了$u$个参数，不论$u<t,u=t$或$u>t$,也不论参数是否独立，每增加一个参数则相应的多产生一个方程，故总共应列出$r+u$个方程。如果在$u$个参数中有$s$个是不独立的，或者说，在这$u$个参数之间存在着$s$个函数关系式，则应列出$s$个限制条件方程，此外，还应列出$c-r+u-s$个一般条件方程，如此就形成了如下的函数模型：$$\begin{array}{c}\underset{c1}{F}(\hat{L},\hat{X})=0\\ \underset{s1}{\varphi }(\hat{X})=0\end{array}$$这就是附有限制条件的条件平差函数模型。其线性形式可概括如下：$$\begin{array}{c}\underset{cn}{A}\underset{n1}{V}+\underset{cu}{B}\underset{u1}{\hat{x}}+\underset{c1}{W}=0\\ \underset{su}{C}\underset{u1}{\hat{x}}+\underset{s1}{W_x}=0\end{array}$$式中$$W=F(L,X^0),W_x=\varphi (X^0)$$其中，$X^0$是未知参数的初始值，一般选择在收敛值的某个开邻域$U$里，使得$U$必须和3维欧式空间$R^3$中的一个开子集是同胚的。
且$c=r+u-s,c>r,s<u$，系数矩阵的秩分别为$$R(A)=c,R(B)=u,R(C)=s$$随机模型为$$\underset{nn}{D}={\sigma }_0^2\underset{nn}{Q}={\sigma }_0^2{\underset{nn}{P}}^{-1}$$在大多数情况下，参数之间是等价的，不需要加权，这时候$P$是单位矩阵。
上面提到过，该模型是四种平差模型的综合：

1. 系数阵中的$B=0,C=0$时，它就变成了条件平差法的函数模型。
2. 系数阵中的$C=0$，它就变成了附有参数的条件平差的函数模型。
3. 系数阵中的$A=-I$和$C=0$时，它就变成了间接平差法的函数模型。
4. 系数阵中的$A=-I$时，就变成了附有限制条件的间接平差法的函数模型。

由此，也称之为“概括平差模型”。
回到正题来简化数学模型，上面的线性模型用拉格朗日乘子进行复合，组成函数如下：$$\varphi =V^{T}PV-2K^T(AV+B\hat{x}+W)-2K_s^T(C\hat{x}+W_x)$$
注意，这里有两个拉格朗日乘子，对于拉格朗日乘子，可参看我的另外一篇文章：拉格朗日乘子法的由来。
由凸优化理论，要计算该极值，需对$V$，$\hat{x}$，$K$，$K_s$分别取极值，即是求偏导并令其为零，得$$\begin{array}{c}\frac{\partial \Phi }{\partial V}=2V^{T}P-2K^{T}A=0\\ \frac{\partial \Phi }{\partial \hat{x}}=-2K^{T}B-2K_s^{T}C=0\\ AV+B\hat{x}+W=0\\ C\hat{x}+W_x=0\end{array}$$
上面四式包含$n$个改正数，$u$个参数，$c$个对应于一般条件方程的联系数以及$s$个对应于限制条件方程的联系数，而方程个数为$c+s+n+u$，即方程个数等于未知数个数，故有唯一解。以上四式成为附有限制条件的条件平差法的基础方程。
解此基础方程，令$$N_{AA}=AQ{A}^T$$得
$$\begin{array}{c}{N}_{AA}K+B\hat{x}+W=0\\ B^{T}K+C^T{K}_s=0\\ C\hat{x}+W_x=0\end{array}$$
以上三式称为附有限制条件的条件平差的法方程。
由法方程得$$\begin{array}{c}K=-N_{AA}^{-1}(W+B\hat{x})\\ \\ B^T{N}_{AA}^{-1}B\hat{x}-C^{T}K+B^T{N}_{AA}^{-1}W=0\end{array}$$若令$$N_{BB}=B^T{N}_{AA}^{-1}B,W_e=B^T{N}_{AA}^{-1}W$$可进一步化简得
$$N_{BB}\hat{x}-C^T{K}_s+W_e=0$$解出$\hat{x}$，得$$\hat{x}=N_{BB}^{-1}(C^T{K}_s-W_e)$$将上式代入$C\hat{x}+W_x=0$，得$$N_{CC}{K}_s-C{N}_{BB}^{-1}{W}_e+W_x=0$$
其中$N_{CC}=C{N}_{BB}^{-1}{C}^T$。解出$K_s$得$$K_s=-N_{CC}^{-1}(W_x-C{N}_{BB}^{-1}{W}_e)$$经整理即得：$$\hat{x}=-(N_{BB}^{-1}-N_{BB}^{-1}{C}^T{N}_{CC}^{-1}C{N}_{BB}^{-1})W_e-N_{BB}^{-1}{C}^T{N}_{CC}^{-1}{W}_x$$
下面是具体的代码实现，其中基本的矩阵运算没有在下面给出，在矩阵算法相关代码，如有需要可以下载。

// <summary>   
/// 附有限制条件的条件平差
/// </summary>     
/// <param name="correction">返回的改正数</param> 
/// <param name="MatrixA"></param>
/// <param name="MatrixB"></param>
/// <param name="MatrixC"></param>
/// <param name="W"></param>
/// <param name="Wx"></param>
/// <param name="C"></param>
/// <param name="N"></param>
/// <param name="U"></param>
/// <param name="S"></param>
//
template<class T>
void GetConditionCorrectionWithCondition(T correction[],const T matrixA[],const T matrixB[],const T matrixC[],const T W[],const T Wx[],int C,int N,int U,int S)
{
    T *transposedA=new T[N*C];
    T *transposedB=new T[U*C];
	T *transposedC=new T[U*S];
    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 *Ncc=new T[S*S];
    T *inverseForNcc=new T[S*S];
	T *We=new T[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(transposedB,inverseForNaa,transit2,U,C,C);
	MultMatrix(transit2,W,We,U,C,1);
	
	T *transit3=new T[S*U];
	MultMatrix(matrixC,inverseForNbb,transit3,S,U,U);
	MatrixTranspose(matrixC,transposedC,S,U);
	MultMatrix(transit3,transposedC,Ncc,S,U,S);
	MatrixAnti(Ncc,inverseForNcc,S);
	
	
	T *transit4=new T[U*S];
	MultMatrix(inverseForNbb,transposedC,transit4,U,U,S);
	T *transit5=new T[U*S];
	MultMatrix(transit4,inverseForNcc,transit5,U,S,S);
	
	T *transit6=new T[U*U];
	MultMatrix(transit5,matrixC,transit6,U,S,U);
	T *transit7=new T[U*U];
	MultMatrix(transit6,inverseForNbb,transit7,U,U);
	MatrixMinus(inverseForNbb,transit7,transit7,U,U);
	T *transit8=new T[U];
	MultMatrix(transit7,We,transit8,U,U,1);
	
	T *transit9=new T[U];
	MultMatrix(transit5,Wx,transit9,U,S,1);
	
	MatrixMinus(transit9,transit9,correction,U,1);
	MatrixMinus(correction,transit8,correction,U,1);
	MatrixMinus(correction,transit9,correction,U,1);
	
	
    delete [] transposedA;
    delete [] transposedB;
	delete [] transposedC;
    delete [] Naa;
    delete [] inverseForNaa;
    delete [] Nbb;
    delete [] inverseForNbb;
	delete [] Ncc;
    delete [] inverseForNcc;
	delete [] We;
    delete [] x;
    delete [] transit1;
    delete [] transit2;
    delete [] transit3;
    delete [] transit4;
    delete [] transit5;
    delete [] transit6;
	delete [] transit7;
    delete [] transit8;
    delete [] transit9;
}

下面进行精度评定，附有限制条件的条件平差的单位权方差估值也是$V^{T}PV$除以它的自由度，即
$${\hat{\sigma }}_0^2=\frac{V^{T}PV}{r}=\frac{V^{T}PV}{c-u+s}$$
另外，应用协因数传播率，可得$$Q_{\hat{x}\hat{x}}=N_{BB}^{-1}-N_{BB}^{-1}{C}^T{N}_{CC}^{-1}C{N}_{BB}^{-1}$$则平差值函数的中误差很容易得到，在这里不想写的太详细了，前面几篇文章详细讨论过。

参考资料：
[1] 误差理论与测量平差基础 武汉大学测绘学院测量平差学科组编著
[2] 变分学讲义 张恭庆著

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