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

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

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

最后做一下误差分析，单位权中误差通过下面公式求得：
σ0^=VTPVn−u+s\hat{\sigma _{0}}=\sqrt{\frac{V^{T}PV}{n-u+s}}σ0^=n−u+sVTPV
其中：nnn为VVV的维数，uuu为拟合向量的维数，sss为条件方程的个数。
平差值得协因数阵根据下面的公式求得：
Qx^x^=NBB−1−NBB−1CTNcc−1CNBB−1Q_{\hat{x}\hat{x}}=N_{BB}^{-1}-N_{BB}^{-1}C^{T}N_{cc}^{-1}CN_{BB}^{-1}Qx^x^=NBB−1−NBB−1CTNcc−1CNBB−1
把Qx^x^Q_{\hat{x}\hat{x}}Qx^x^所有的元素相加得到平差值x^\hat{x}x^的协因数：
Qx^=FTQx^x^FQ_{\hat{x}}=F^{T}Q_{\hat{x}\hat{x}}FQx^=FTQx^x^F
其中：F=(11⋯1)TF=\begin{pmatrix} 1 & 1& \cdots & 1 \end{pmatrix}^{T}F=(11⋯1)T
最后由下面公式可算得平差值函数的中误差：
σ^x^=σ0^Qx^\hat{\sigma }_{\hat{x}}=\hat{\sigma _{0}}\sqrt{Q_{\hat{x}}}σ^x^=σ0^Qx^

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

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

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

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

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

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

测量平差之附有限制条件的间接平差相关推荐

测量平差之附有限制条件的条件平差（概括平差模型）
接下来要写的可看成是对前面讲过的四种经典平差方法的概括模型. 从四种基本平差算法的函数模型来看,主要包括如下两种类型的条件方程: { F ( L ^ , X ^ ) = 0 ϕ ( X ^ ) = 0 ...
附有限制条件的间接平差类(+抗差估计)
附有限制条件的间接平差类(+抗差估计) 间接平差模型,有详细的注释,需要的自取头文件是上个博客写的头文件这是个 FX_adj.h文件文件 #if !defined(CHDADJ_FX_ADJ_H_ ...
计算机考研845大纲,中国地质大学2018考研大纲：845测量平差
<中国地质大学2018考研大纲:845测量平差>由会员分享,可在线阅读,更多相关<中国地质大学2018考研大纲:845测量平差(4页珍藏版)>请在人人文库网上搜索. 1.中国地 ...
java二维数组两个框代表什么_在java语言中，二维数组的两个中括号[][]分别表示（）和（）。...
[判断题]一元统计分析是研究一个随机变量统计规律的学科. [判断题]合伙创业的成功率一般低于独资创业的成功率. [单选题]镜检时呈"竹节状"排列的是 [判断题]酵母菌的菌落与放线菌 ...
基于Qt、Eigen的四种平差模型计算器
可实现 1.条件平差 2.间接平差 3.附有限制条件的间接平差 4.附有参数的条件平差效果: 直接上代码 .h #pragma once#include <QtWidgets/QMainWin ...
matlab解算平差实例,MATLAB软件在测量平差解算中的应用
MAT LAB 软件在测量平差解算中的应用胡远新,赵奋军 (浙江省第七地质大队, 浙江丽水市 323000) 摘要:阐述了4种经典的测量平差的解算过程,并结合平差实例,说明了MAT LAB 在测 ...
matlab边角网间接平差计算,12.2测边网与边角网间接平差
§12.2测边网与边角网间接平差测边网或边角网(其实是测方向),都是把所测方向和边作为平差元素,因此这类网只要按边长观测值.方向观测值列出误差方程式,就可组成法方程式,-- 误差方程式的总数等于网中 ...
测量平差MATLAB方法
条件平差 P=input('请输入权阵 '); A=input('请输入矩阵A '); w=input('请输入矩阵w '); n=A*inv(P)*A.'; k=-inv(n)*w; v=inv(P ...
协方差公式性质证明过程_论文推荐 | 刘志平：等价条件平差模型的方差-协方差分量最小二乘估计方法...
<测绘学报> 构建与学术的桥梁拉近与权威的距离等价条件平差模型的方差-协方差分量最小二乘估计方法刘志平1, 朱丹彤1, 余航1, 张克非1,2 1. 中国矿业大学环境与测绘学院, 江 ...

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

测量平差之附有限制条件的间接平差相关推荐

最新文章

热门文章