测量平差之条件平差

测量平差是一套完整的理论，在前面我分别简要的叙述并实现了间接平差和附有限制条件的间接平差，这是两类最常用的测量平差算法，希望能给这一方面的初学者一点点有意义的参考。下面我准备再写几篇文章来叙述一下余下经典平差算法，这一篇就说一下条件平差。
在拟合和测量工作中，为了能达到平差目的，一般要作多的观测。如果观测数据确实有限，添加条件方程一般也可以达到目的。添加条件方程本质上是降维，增加了“偏执”使得必要观测数减少。同时收敛路径也变得单调，就好像过“独木桥”，一旦独木桥无法通过，算法将无法收敛到规定误差范围内。所以添加条件方程是把“双刃剑”。但其在工业测量及最小二乘拟合中依然是一个重要的手段，很多情况下维它不行。以条件方程为函数模型的平差方法，就是条件平差。
一般而言，如果有$n$个观测值$\underset{n1}{L}$，$t$个必要观测，则应列出$r=n-t$个条件方程，即$$F(\tilde{L})=0$$如果条件方程为线性形式，可直接写为$$\underset{rn}{A}\widetilde{\underset{n1}{L}}+\underset{r1}{A_0}=0$$将$$\tilde{L}=L+V$$代入上式，并令$$W=AL+A_0$$则可把公式简化为$$AV+W=0$$此即为条件平差的函数模型。其自由度即为条件方程的个数。另外，条件平差也可看成是附加限制条件的间接平差的特例($B$矩阵为单位矩阵，$l$为0)，虽然这样理解有些别扭，但在工业测量中附加限制条件的间接平差会被更多的用到。另外，条件平差与间接平差的区别：条件平差不含有独立参数，并且一般表示为观测量的隐函数形式。
依附最小二乘原理，平差准则为$$V^{{}^T}PV=min$$其中$P$是观测向量的权阵，大多数情况下可设为单位矩阵。
按拉格朗日乘子理论，设乘子$\underset{r1}{K}={\left[\begin{array}{cccc}{k}_a& k_b& \cdots & k_r\end{array}\right]}^T$，组成函数$$\Phi =V^{T}PV-2K^T(AV+W)$$
将$\Phi$对$V$求一阶导数，并令其为0，得$$\frac{d\Phi }{dV}=2V^{T}P-2K^{T}A=0$$两边转置，得$$PV=A^{T}K$$再用$P^{-1}$左乘上式两端，得改正数$V$的计算公式为$$V=P^{-1}{A}^{T}K=Q{A}^{T}K$$上式为改正数方程。将上式代入$AV+W=0$，得$$AQ{A}^{T}K+W=0$$令$$N_{aa}=AQ{A}^T=A{P}^{-1}{A}^T$$则有$$N_{aa}K+W=0$$为法方程，其中$N_{aa}$必须是非奇异矩阵，由此可得$K$的唯一解$$K=-N_{aa}^{-1}W$$再把它代回到上面的改正数计算公式，得$$V=-Q{A}^T{N}_{aa}^{-1}W$$求出改正数$V$值，再求平差值$\hat{L}=L+V$，这样就完成了按条件平差求平差值的工作。
下面是具体的代码实现，其中基本的矩阵运算没有在下面给出，在矩阵算法相关代码，如有需要可以下载。

template<class T1,class T2,class T3>
void  __declspec(dllexport) MatrixMinus(T1* M,T2* M1,int nrows,int ncols)
{
	for(int i=0;i<nrows;++i)
		for(int j=0;j<ncols;++j)
			M1[i*ncols+j]=-1*M[i*ncols+j];   
}

// <summary>   
/// 条件平差
/// </summary>     
/// <param name="correction">返回的改正数</param> 
/// <param name="MatrixA"></param>
/// <param name="w"></param>
/// <param name="R"></param>
/// <param name="N"></param>
//
template<class T>
void GetConditionCorrection(T correction[],const T matrixA[],const T w[],int R,int N)
{
	T *transposedA=new T[R*N];
	T *Naa=new T[R*R];
	T *inverseForNaa=new T[R*R];
	T *K=new T[R];
	MatrixTranspose(matrixA,transposedA,R,N);
	MultMatrix(matrixA,transposedA,Naa,R,N,R);
	MatrixAnti(Naa,inverseForNaa,R);
	MultMatrix(inverseForNaa,w,K,R,R,1);
	MatrixMinus(K,K,R,1);
	MultMatrix(transposedA,K,correction,N,R,1);
	delete [] transposedA;
	delete [] Naa;
	delete [] inverseForNaa;
	delete [] K;
}

最后进行误差分析，设平差值函数为$$\hat{\phi }=f({\hat{L}}_1,{\hat{L}}_1,\cdots ,{\hat{L}}_n)$$对两边进行全微分$$d\hat{\phi }=(\frac{\partial f}{\partial {\hat{L}}_1}{)}_{0}d{\hat{L}}_1+(\frac{\partial f}{\partial {\hat{L}}_2}{)}_{0}d{\hat{L}}_2+\cdots +(\frac{\partial f}{\partial {\hat{L}}_n}{)}_{0}d{\hat{L}}_n$$令$(\frac{\partial f}{\partial {\hat{L}}_i}{)}_0=f_i$，则上式为$$d\hat{\phi }=f_{1}d{\hat{L}}_1+f_{2}d{\hat{L}}_2+\cdots +f_{n}d{\hat{L}}_n$$
设$f={\left[\begin{array}{cccc}{f}_1& f_2& \cdots & f_n\end{array}\right]}^T$，$d\hat{L}={\left[\begin{array}{cccc}d{\hat{L}}_1& d{\hat{L}}_2& \cdots & d{\hat{L}}_n\end{array}\right]}^T$，则上式可简写为$d\hat{\phi }=f^{T}d\hat{L}$，由此即得$$Q_{\hat{\phi }\hat{\phi }}=f^T{Q}_{\hat{L}\hat{L}}f$$式中，$Q_{\hat{L}\hat{L}}$为平差值$\hat{L}$的协因数阵。并且可得$$Q_{\hat{L}\hat{L}}=Q-Q_{VV}=Q-Q{A}^T{N}_{aa}^{-1}AQ$$代入$Q_{\hat{\phi }\hat{\phi }}$，得$$Q_{\hat{\phi }\hat{\phi }}=f^{T}Qf-(AQf{)}^T{N}_{aa}^{-1}AQf$$当平差值函数已知，即可按照上式计算函数$\hat{\phi }$协因数。
平差值函数的中误差为$${\hat{\sigma }}_{\hat{\phi }}={\hat{\sigma }}_0\sqrt{Q_{\hat{\phi }\hat{\phi }}}$$

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

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