测量平差之间接平差

今年3月份的时候我在公司写一个算法接触到测量平差理论，感慨于它的算法的优美，利用最小二乘而又独成体系，在测量拟合中发挥着巨大的作用。但网上关于它的描述甚少，更不用说好的实现代码，我将写两篇文章来说一下间接平差和附加限制条件的间接平差，这一篇先说间接平差。
在测量拟合中，由于观测结果不可避免地存在着误差，因此，如何处理带有误差的观测值，找出待求量的最佳估值，是最小二乘平差所研究的内容。在一个平差问题中，当所选的独立参数$\hat{X}$的个数等于必要观测数$t$时，可将每个观测值表达成这$t$个参数的函数，组成观测方程，这种以观测方程为函数模型的平差方法，就是间接平差。用数学语言描述就是：通过在映射$B$的像上的平差拉回到原像来间接地反映原像上的平差过程的算法。
一般而言，如果某平差问题有$n$个观测值，$t$个必要观测值，选择$t$个独立量作为平差参数$\underset{t1}{\hat{X}}$，则每个观测值必定可以表达成这$t$个参数的函数，即有：$$\underset{n1}{\hat{L}}=F\left(\hat{X}\right)$$如果表达式是线性的，一般为：$$\underset{n1}{\hat{L}}=\underset{nt}{B}\underset{t1}{\hat{X}}+\underset{n1}{d}$$其中自由度是：$r=n-t$。
平差时，一般对参数$\hat{X}$都要取近似值$X^0$,令$$\hat{X}=X^0+\hat{x}$$代入上式，并令$$l=L-\left(B{X}^0+d\right)=L-L^0$$$L^0=B{X}^{{}^0}+d$为观测值的近似值，所以$l$是观测值与其近似值之差，由此可得误差方程$$V=B\hat{x}-l$$平差准则为$$V^{T}PV=min$$以上是基本原理，下面就逐步推到得到平差公式。
设有$n$个观测值方程为：$$\begin{array}{c}{L}_1+v_1=a_1\widehat{X_1}+b_1\widehat{X_2}+\cdots +t_1\widehat{X_t}+d_1\\ L_2+v_2=a_2\widehat{X_1}+b_2\widehat{X_2}+\cdots +t_2\widehat{X_t}+d_2\\ \cdots \\ L_n+v_n=a_n\widehat{X_1}+b_n\widehat{X_2}+\cdots +t_n\widehat{X_t}+d_n\end{array}$$令$$\widehat{X_j}=X_j^0+\widehat{x_j}\left(j=1,2,\cdots ,t\right)$$$$l_i=L_i-\left(a_i\widehat{X_1}+b_i\widehat{X_2}+\cdots +t_i\widehat{X_t}+d_i\right)\left(i=1,2,\cdots ,n\right)$$可得最终的误差方程为$$v_i=a_i\widehat{x_1}+b_i\widehat{x_2}+\cdots +t_i\widehat{x_t}-l_i\left(i=1,2,\cdots ,n\right)$$其中：$$\underset{nt}{B}=\left[\begin{array}{cccc}{a}_1& b_1& \cdots & t_1\\ a_2& b_2& \cdots & t_2\\ ⋮& ⋮& & ⋮\\ a_n& b_n& \cdots & t_n\end{array}\right]$$$$\underset{n1}{V}={\left[\begin{array}{cccc}{v}_1& v_2& \cdots & v_n\end{array}\right]}^T$$$$\widehat{\underset{t1}{x}}={\left[\begin{array}{cccc}{\hat{x}}_1& {\hat{x}}_2& \cdots & {\hat{x}}_t\end{array}\right]}^T$$$$\underset{n1}{l}={\left[\begin{array}{cccc}{l}_1& l_2& \cdots & l_n\end{array}\right]}^T$$按最小二乘原理，可得：
$$\frac{\partial V^{T}PV}{\partial \hat{x}}=2V^{T}P\frac{\partial V}{\partial \hat{x}}=V^{T}PB=0$$
上述方程与$V=B\hat{x}-l$联立，解得：$$B^{T}PB\hat{x}-B^{T}Pl=0$$令$$N_{BB}=B^{T}PB,W=B^{T}Pl$$上式可简写为：$$N_{BB}\hat{x}-W=0$$式中系数阵$N_{BB}$为满秩，$\hat{x}$有唯一解，上式称为间接平差的法方程。解之得：$$\hat{x}=N_{BB}^{-1}W$$从而平差结果为：$$\hat{X}=X^0+\hat{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;
}

最后再说一下精度评定，首先通过下面的公式得到单位权中误差：$${\hat{\sigma }}_0=\sqrt{\frac{V^{T}PV}{n-t}}$$其中$n$为方程数，$t$为拟合向量的维数。平差值的协因数阵根据下面公式求得：$$Q_{\hat{x}\hat{x}}=N_{BB}^{-1}$$假定间接平差问题中有$t$个参数，设参数的函数为：$$\hat{\phi }=\varphi \left(\widehat{X_1},\widehat{X_2},\cdots ,\widehat{X_t}\right)$$为求函数$\phi $的中误差，先对上式全微分得权函数式为：$$d\hat{\phi }=f_{1}d\widehat{X_1}+f_{2}d\widehat{X_2}+\cdots +f_{t}d\widehat{X_t}$$令$F^{T}=\begin{bmatrix}
f_{1} & f_{2} & \cdots & f_{t}
\end{bmatrix}$，则函数$\hat{\varphi }$的协因数阵为：$$Q_{\hat{\phi }\hat{\phi }}=F^T{Q}_{\hat{x}\hat{x}}F=F^T{N}_{BB}^{-1}F$$
最后由下面公式可算得平差值函数的中误差：$${\hat{\sigma }}_{\hat{\phi }}={\hat{\sigma }}_0\times \sqrt{Q_{\hat{\phi }\hat{\phi }}}$$

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

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