今年3月份的时候我在公司写一个算法接触到测量平差理论,感慨于它的算法的优美,利用最小二乘而又独成体系,在测量拟合中发挥着巨大的作用。但网上关于它的描述甚少,更不用说好的实现代码,我将写两篇文章来说一下间接平差和附加限制条件的间接平差,这一篇先说间接平差。
在测量拟合中,由于观测结果不可避免地存在着误差,因此,如何处理带有误差的观测值,找出待求量的最佳估值,是最小二乘平差所研究的内容。在一个平差问题中,当所选的独立参数
X
^
\hat{X}
X^的个数等于必要观测数
t
t
t时,可将每个观测值表达成这
t
t
t个参数的函数,组成观测方程,这种以观测方程为函数模型的平差方法,就是间接平差。用数学语言描述就是:通过在映射
B
B
B的像上的平差拉回到原像来间接地反映原像上的平差过程的算法。
一般而言,如果某平差问题有
n
n
n个观测值,
t
t
t个必要观测值,选择
t
t
t个独立量作为平差参数
X
^
t
1
\underset{t1}{\hat{X}}
t1X^,则每个观测值必定可以表达成这
t
t
t个参数的函数,即有:
L
^
n
1
=
F
(
X
^
)
\underset{n1}{\hat{L}}=F\left ( \hat{X} \right )
n1L^=F(X^)如果表达式是线性的,一般为:
L
^
n
1
=
B
n
t
X
^
t
1
+
d
n
1
\underset{n1}{\hat{L}}=\underset{nt}{B}\underset{t1}{\hat{X}}+\underset{n1}{d}
n1L^=ntBt1X^+n1d其中自由度是:
r
=
n
−
t
r=n-t
r=n−t。
平差时,一般对参数
X
^
\hat{X}
X^都要取近似值
X
0
X^{0}
X0,令
X
^
=
X
0
+
x
^
\hat{X}=X^{0}+\hat{x}
X^=X0+x^代入上式,并令
l
=
L
−
(
B
X
0
+
d
)
=
L
−
L
0
l=L- \left ( BX^{0}+d \right ) =L-L^{0}
l=L−(BX0+d)=L−L0
L
0
=
B
X
0
+
d
L^{0}=BX^{^{0}}+d
L0=BX0+d为观测值的近似值,所以
l
l
l是观测值与其近似值之差,由此可得误差方程
V
=
B
x
^
−
l
V=B\hat{x}-l
V=Bx^−l平差准则为
V
T
P
V
=
m
i
n
V^{T}PV=min
VTPV=min以上是基本原理,下面就逐步推到得到平差公式。
设有
n
n
n个观测值方程为:
L
1
+
v
1
=
a
1
X
1
^
+
b
1
X
2
^
+
⋯
+
t
1
X
t
^
+
d
1
L
2
+
v
2
=
a
2
X
1
^
+
b
2
X
2
^
+
⋯
+
t
2
X
t
^
+
d
2
⋯
L
n
+
v
n
=
a
n
X
1
^
+
b
n
X
2
^
+
⋯
+
t
n
X
t
^
+
d
n
\begin{matrix} L_{1}+v_{1}=a_{1}\hat{X_{1}}+b_{1}\hat{X_{2}}+\cdots +t_{1}\hat{X_{t}}+d_{1}\\ L_{2}+v_{2}=a_{2}\hat{X_{1}}+b_{2}\hat{X_{2}}+\cdots +t_{2}\hat{X_{t}}+d_{2}\\ \cdots \\ L_{n}+v_{n}=a_{n}\hat{X_{1}}+b_{n}\hat{X_{2}}+\cdots +t_{n}\hat{X_{t}}+d_{n} \end{matrix}
L1+v1=a1X1^+b1X2^+⋯+t1Xt^+d1L2+v2=a2X1^+b2X2^+⋯+t2Xt^+d2⋯Ln+vn=anX1^+bnX2^+⋯+tnXt^+dn令
X
j
^
=
X
j
0
+
x
j
^
(
j
=
1
,
2
,
⋯
,
t
)
\hat{X_{j}}=X_{j}^{0}+\hat{x_{j}}\left ( j=1,2,\cdots ,t \right )
Xj^=Xj0+xj^(j=1,2,⋯,t)
l
i
=
L
i
−
(
a
i
X
1
^
+
b
i
X
2
^
+
⋯
+
t
i
X
t
^
+
d
i
)
(
i
=
1
,
2
,
⋯
,
n
)
l_{i}=L_{i}-\left ( a_{i}\hat{X_{1}}+b_{i}\hat{X_{2}}+\cdots +t_{i}\hat{X_{t}}+d_{i} \right )\left ( i=1,2,\cdots ,n \right )
li=Li−(aiX1^+biX2^+⋯+tiXt^+di)(i=1,2,⋯,n)可得最终的误差方程为
v
i
=
a
i
x
1
^
+
b
i
x
2
^
+
⋯
+
t
i
x
t
^
−
l
i
(
i
=
1
,
2
,
⋯
,
n
)
v_{i}= a_{i}\hat{x_{1}}+b_{i}\hat{x_{2}}+\cdots +t_{i}\hat{x_{t}}-l_{i}\left ( i=1,2,\cdots ,n \right )
vi=aix1^+bix2^+⋯+tixt^−li(i=1,2,⋯,n)其中:
B
n
t
=
[
a
1
b
1
⋯
t
1
a
2
b
2
⋯
t
2
⋮
⋮
⋮
a
n
b
n
⋯
t
n
]
\underset{nt}{B}=\begin{bmatrix} a_{1}& b_{1} & \cdots & t_{1}\\ a_{2}& b_{2}& \cdots& t_{2}\\ \vdots& \vdots & & \vdots\\ a_{n}& b_{n}& \cdots& t_{n} \end{bmatrix}
ntB=⎣⎢⎢⎢⎡a1a2⋮anb1b2⋮bn⋯⋯⋯t1t2⋮tn⎦⎥⎥⎥⎤
V
n
1
=
[
v
1
v
2
⋯
v
n
]
T
\underset{n1}{V}=\begin{bmatrix} v_{1} & v_{2} & \cdots & v_{n} \end{bmatrix}^{T}
n1V=[v1v2⋯vn]T
x
t
1
^
=
[
x
^
1
x
^
2
⋯
x
^
t
]
T
\hat{\underset{t1}{x}}=\begin{bmatrix} \hat{x}_{1}& \hat{x}_{2}& \cdots & \hat{x}_{t} \end{bmatrix}^{T}
t1x^=[x^1x^2⋯x^t]T
l
n
1
=
[
l
1
l
2
⋯
l
n
]
T
\underset{n1}{l}=\begin{bmatrix} l_{1}& l_{2}& \cdots & l_{n} \end{bmatrix}^{T}
n1l=[l1l2⋯ln]T按最小二乘原理,可得:
∂
V
T
P
V
∂
x
^
=
2
V
T
P
∂
V
∂
x
^
=
V
T
P
B
=
0
\frac{\partial V^{T}PV}{\partial \hat{x}}=2V^{T}P\frac{\partial V}{\partial \hat{x}}=V^{T}PB=0
∂x^∂VTPV=2VTP∂x^∂V=VTPB=0
上述方程与
V
=
B
x
^
−
l
V=B\hat{x}-l
V=Bx^−l联立,解得:
B
T
P
B
x
^
−
B
T
P
l
=
0
B^{T}PB\hat{x}-B^{T}Pl=0
BTPBx^−BTPl=0令
N
B
B
=
B
T
P
B
,
W
=
B
T
P
l
N_{BB}=B^{T}PB,W=B^{T}Pl
NBB=BTPB,W=BTPl上式可简写为:
N
B
B
x
^
−
W
=
0
N_{BB}\hat{x}-W=0
NBBx^−W=0式中系数阵
N
B
B
N_{BB}
NBB为满秩,
x
^
\hat{x}
x^有唯一解,上式称为间接平差的法方程。解之得:
x
^
=
N
B
B
−
1
W
\hat{x}=N_{BB}^{-1}W
x^=NBB−1W从而平差结果为:
X
^
=
X
0
+
x
^
\hat{X}=X^{0}+\hat{x}
X^=X0+x^其中,
P
P
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;
}
最后再说一下精度评定,首先通过下面的公式得到单位权中误差:
σ
^
0
=
V
T
P
V
n
−
t
\hat{\sigma }_{0}=\sqrt{\frac{V^{T}PV}{n-t}}
σ^0=n−tVTPV其中
n
n
n为方程数,
t
t
t为拟合向量的维数。平差值的协因数阵根据下面公式求得:
Q
x
^
x
^
=
N
B
B
−
1
Q_{\hat{x}\hat{x}}=N_{BB}^{-1}
Qx^x^=NBB−1假定间接平差问题中有
t
t
t个参数,设参数的函数为:
φ
^
=
ϕ
(
X
1
^
,
X
2
^
,
⋯
,
X
t
^
)
\hat{\varphi }=\phi \left ( \hat{X_{1}},\hat{X_{2}},\cdots ,\hat{X_{t}} \right )
φ^=ϕ(X1^,X2^,⋯,Xt^)为求函数$\phi
的
中
误
差
,
先
对
上
式
全
微
分
得
权
函
数
式
为
:
的中误差,先对上式全微分得权函数式为:
的中误差,先对上式全微分得权函数式为:
d
φ
^
=
f
1
d
X
1
^
+
f
2
d
X
2
^
+
⋯
+
f
t
d
X
t
^
d\hat{\varphi }=f_{1}d\hat{X_{1}}+f_{2}d\hat{X_{2}}+\cdots+ f_{t}d\hat{X_{t}}
dφ^=f1dX1^+f2dX2^+⋯+ftdXt^
令
令
令F^{T}=\begin{bmatrix}
f_{1} & f_{2} & \cdots & f_{t}
\end{bmatrix}
,
则
函
数
,则函数
,则函数\hat{\varphi }
的
协
因
数
阵
为
:
的协因数阵为:
的协因数阵为:
Q
φ
^
φ
^
=
F
T
Q
x
^
x
^
F
=
F
T
N
B
B
−
1
F
Q_{\hat{\varphi }\hat{\varphi }}=F^{T}Q_{\hat{x}\hat{x}}F=F^{T}N_{BB}^{-1}F
Qφ^φ^=FTQx^x^F=FTNBB−1F$
最后由下面公式可算得平差值函数的中误差:
σ
^
φ
^
=
σ
^
0
×
Q
φ
^
φ
^
\hat{\sigma }_{\hat{\varphi}} = \hat{\sigma }_{0}\times \sqrt{Q_{\hat{\varphi }\hat{\varphi }}}
σ^φ^=σ^0×Qφ^φ^
参考文献:
[1] 误差理论与测量平差基础 武汉大学测绘学院测量平差学科组编著
转载请注明出处:http://blog.csdn.net/fourierFeng/article/details/47272167