测量平差之附有限制条件的间接平差

上一篇文章中我简单的讲了一下间接平差，这篇文章我接着说一下附有限制条件的情况。
所谓附加的限制条件，就是在平差拟合的基础上，规范收敛路径，使其的收敛过程完全按照某一规范进行，这样就需要一个函数去“描述”这样的规范。附加限制条件是把双刃剑。一方面，它可以非常严格的规范拟合过程；另一方面，它又会拖累收敛速度，甚至影响收敛结果，使得收敛变得不稳定。
从原理上讲，附有限制条件的间接平差就是在间接平差的基础上利用拉格朗日乘子。说到拉格朗日乘子理论，我就多说几句，国内首先对其做研究的应该是钱伟长，在数学上，它涉及到泛函分析和偏微分方程的相关理论。虽然大部分理工科本科生都能用拉格朗日乘子计算条件极值（考研数学中常考的知识点），但仅仅是会用远远不够，甚至可以说是不必需的。能了解由来才是最重要的。对此我还是存有疑问，网上有几篇文章对它在工程上讲的很详细，但也没有涉及到它的由来问题。目前来说，可参考卡罗需-库恩-塔克条件（英文原名：Karush-Kuhn-Tucker Conditions常见别名：Kuhn-Tucker，KKT条件，Karush-Kuhn-Tucker最优化条件，Karush-Kuhn-Tucker条件，Kuhn-Tucker最优化条件，Kuhn-Tucker条件）和钱伟长的《广义变分原理》。对拉格朗日乘子我不在这里赘述，希望不久我有机会单独写一篇文章来说说拉格朗日乘子理论。
下面回到正题，一般而言，附有限制条件的间接平差可组成下列方程：
$$\begin{array}{c}\underset{n1}{\tilde{L}}=F(\underset{n1}{\tilde{X}})\\ \underset{s1}{\Phi }(\tilde{X})=0\end{array}\}$$
这里，我们只讨论线性形式，上面的方程可具体化为：$$\begin{array}{c}\underset{n1}{\tilde{V}}=\underset{nu}{B}\underset{u1}{\hat{x}}-\underset{n1}{l}\\ \underset{su}{C}\underset{u1}{\hat{x}}-\underset{s1}{W_x}=0\end{array}\}$$
其中：$$R(B)=u,R(C)=s,u<n,s<u$$即$B$为列满秩矩阵，$C$为行满秩矩阵。
按拉格朗日乘数法组成函数：$$\Phi =V^{T}PV+2K_s^T(C\hat{x}-W_x)$$
式中$K_s^T$是拉格朗日乘子，前面系数2是为了让结果更简洁。
$V$是$\hat{x}$的显函数，为求$\Phi$的极小，将其对$\hat{x}$取偏导并令其为零，则有
$$\frac{\partial \Phi }{\partial \hat{x}}=2V^{T}P\frac{\partial V}{\partial \hat{x}}+2K_s^{T}C=2V^{T}PB+2K_s^{T}C=0$$
转置后得$$B^{T}PV+C^T{K}_s=0$$
将$\underset{n1}{\tilde{V}}=\underset{nu}{B}\underset{u1}{\hat{x}}-\underset{n1}{l}$代入上式，得：$$B^{T}PB\hat{x}+C^T{K}_s-B^{T}Pl=0$$
令：$$N_{BB}=B^{T}PB,W_{u1}=B^{T}Pl$$
故上式可写成：$$N_{BB}\hat{x}+C^T{K}_s-W_{u1}=0$$再加上$$C\hat{x}-W_x=0$$
称为附有限制条件的间接平差法的法方程。
再令：$N_{cc}=C{N}_{BB}^{-1}{C}^T$，则上面两式可写为：
$$N_{cc}{K}_s-（C{N}_{BB}^{-1}{W}_{u1}-W_x）=0$$
由于$N_{cc}$是一$s$阶的满秩对称方阵，于是$$K_s=N_{cc}^{-1}（C{N}_{BB}^{-1}{W}_{u1}-W_x）$$
经整理得：
$$\hat{x}=（N_{BB}^{-1}-N_{BB}^{-1}{C}^T{N}_{cc}^{-1}C{N}_{BB}^{-1}）W_{u1}+N_{BB}^{-1}{C}^T{N}_{cc}^{-1}{W}_x$$
还有的文献上是这样推导的：
$$\left(\begin{array}{cc}{B}^{T}PB& C^T\\ C& 0\end{array}\right)\left(\begin{array}{c}\hat{x}\\ K_s\end{array}\right)=\left(\begin{array}{c}{B}^{T}Pl\\ W_x\end{array}\right)$$
这样其实等价于把条件方程组加到平差方程后面共同组成间接平差，条件方程组就不再是死的条件，而变成了平差的一部分；所以$N_{BB}$不再要求必须可逆，这样在观测点不足的情况下也能做平差，不好的地方是构造了一个大矩阵，增大了时间复杂度和算法的不确定性。
最后即可求出：$$\hat{X}=X^0+\hat{x}$$
其中，$P$是观测数据的协因数阵的逆矩阵，一般可认为是单位矩阵。
下面是具体的代码实现，其中基本的矩阵运算没有在下面给出，在矩阵算法相关代码，如有需要可以下载。

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;
}

template<class T>
void GetNcc(T Ncc[],const T matrixA[],const T nbb[],int S,int U)
{
	T *inverseForNbb=new T[U*U];
	T *transposedA=new T[S*U];
	T *transit=new T[S*U];
	MatrixTranspose(matrixA,transposedA,S,U);
	MatrixAnti(nbb,inverseForNbb,U);
	MultMatrix(matrixA,inverseForNbb,transit,S,U,U);
	MultMatrix(transit,transposedA,Ncc,S,U,S);

	delete [] inverseForNbb;
	delete [] transposedA;
	delete [] transit;
}

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;
}



/// <summary>   
/// 附有限制条件的间接平差
/// </summary>     
/// <param name="correction">返回的改正数</param> 
/// <param name="MatrixA"></param>
/// <param name="w"></param>
/// <param name="matrixB"></param>
/// <param name="l"></param>
/// <param name="N"></param>
/// <param name="U"></param>
/// <param name="S"></param>
//
template<class T>
void GetCorrectionWithCondition(T correction[],const T matrixA[],const T w[],const T *matrixB,const T *l,int N,int U,int S)
{
	T *nbb=new T[U*U];
	T *inverseForNbb=new T[U*U];
	T *transposedA=new T[S*U];
	T *Ncc=new T[S*S];
	T *inverseForNcc=new T[S*S];
	T *Wu1=new T[U];

	GetNBB(nbb,matrixB,N,U);
	MatrixAnti(nbb,inverseForNbb,U);
	MatrixTranspose(matrixA,transposedA,S,U);
	GetNcc(Ncc,matrixA,nbb,S,U);
	MatrixAnti(Ncc,inverseForNcc,S);
	GetWu1(Wu1,matrixB,l,N,U);

	T *transit1=new T[S*U];

	MultMatrix(inverseForNbb,transposedA,transit1,U,U,S);

	T *transit2=new T[S*U];

	MultMatrix(transit1,inverseForNcc,transit2,U,S,S);

	T *transit3=new T[U*U];

	MultMatrix(transit2,matrixA,transit3,U,S,U);

	T *transit4=new T[U*U];

	MultMatrix(transit3,inverseForNbb,transit4,U,U);

	T *transit5=new T[U*U];

	MatrixMinus(inverseForNbb,transit4,transit5,U,U);

	T *transit6=new T[U];

	MultMatrix(transit5,Wu1,transit6,U,U,1);

	T *transit7=new T[S*U];

	MultMatrix(inverseForNbb,transposedA,transit7,U,U,S);

	T *transit8=new T[S*U];

	MultMatrix(transit7,inverseForNcc,transit8,U,S,S);

	T *transit9=new T[U];

	MultMatrix(transit8,w,transit9,U,S,1);

	MatrixMinus(transit6,transit9,correction,U,1);

	delete [] nbb;
	delete [] inverseForNbb;
	delete [] transposedA;
	delete [] Ncc;
	delete [] inverseForNcc;
	delete [] Wu1;
	delete [] transit1;
	delete [] transit2;
	delete [] transit3;
	delete [] transit4;
	delete [] transit5;
	delete [] transit6;
	delete [] transit7;
	delete [] transit8;
	delete [] transit9;
}

最后做一下误差分析，单位权中误差通过下面公式求得：
$$\widehat{{\sigma }_0}=\sqrt{\frac{V^{T}PV}{n-u+s}}$$
其中：$n$为$V$的维数，$u$为拟合向量的维数，$s$为条件方程的个数。
平差值得协因数阵根据下面的公式求得：
$$Q_{\hat{x}\hat{x}}=N_{BB}^{-1}-N_{BB}^{-1}{C}^T{N}_{cc}^{-1}C{N}_{BB}^{-1}$$
把$Q_{\hat{x}\hat{x}}$所有的元素相加得到平差值$\hat{x}$的协因数：
$$Q_{\hat{x}}=F^T{Q}_{\hat{x}\hat{x}}F$$
其中：$F={\left(\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right)}^T$
最后由下面公式可算得平差值函数的中误差：
$${\hat{\sigma }}_{\hat{x}}=\widehat{{\sigma }_0}\sqrt{Q_{\hat{x}}}$$

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

[2] 工业测量拟合 王解先 季凯敏著

[3]广义变分原理 钱伟长著

[4]关于拉格朗日乘子法与KKT条件

[5]维基百科：卡羅需－庫恩－塔克條件

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