在之前我们讨论了条件平差,但在有时我们无法列出足够的无关的条件方程,这时我们就需要添加参数,于是就有了附有参数的条件方程
附有参数的条件平差也可看作是一种间接平差和条件平差的混合模型。增加条件方程是“降维”,它可以把复杂的问题变得简单;添加参数是“升维”,它可以把问题变得理论上可分的概率更大。两者并不矛盾,而且相辅相成。我们中学初等几何中常用的添加辅助线、辅助面的方法其实就是附加参数以增加条件方程(这段话摘自这里,感觉挺有意思的)。该文提到了一个Cover定理来理解,也就是某些问题在低维不好解决就可以放到高维去解决,增加维度和这里的增加拟合参数意思差不多。比如下图,我们要分开这两个颜色的球很难
但是我们增加一个维度,将这个问题放到三维来解决可能就会更容易一些
所以我们附有参数的条件平差也是这样,我们通过一个简单的题感受一下
如题
我们很直观就能看出,总的观测数
n
=
6
n=6
n=6,必要观测值为
t
=
3
t=3
t=3,就可以知道多余观测数
r
=
3
r=3
r=3,此时我们根据图形,只能列出如下两个图像条件
∠
1
+
∠
2
+
∠
3
=
180
∠
4
+
∠
5
+
∠
6
=
180
\angle1+\angle2+\angle3=180\\ \angle4+\angle5+\angle6=180
∠1+∠2+∠3=180∠4+∠5+∠6=180
这时候还差一个条件,于是我们连接
A
B
AB
AB,设
∠
B
A
D
=
X
\angle BAD=X
∠BAD=X,现在我们多了一个参数,使得参数的个数
u
=
1
u=1
u=1,现在我们就需要列出
c
=
r
+
u
=
4
c=r+u=4
c=r+u=4个条件方程。我们可以列出一个极条件(看这里)
s
i
n
(
∠
5
)
s
i
n
(
∠
2
+
∠
4
)
s
i
n
X
s
i
n
(
∠
1
−
X
)
s
i
n
(
∠
4
)
s
i
n
(
∠
3
+
∠
5
)
=
1
\frac{sin(\angle5)sin(\angle2 + \angle4)sinX}{sin(\angle1 - X)sin(\angle4)sin(\angle3 + \angle5)}=1
sin(∠1−X)sin(∠4)sin(∠3+∠5)sin(∠5)sin(∠2+∠4)sinX=1
还有一个固定边条件,这个就是个正弦定理
A
B
s
i
n
X
B
D
s
i
n
(
∠
3
+
∠
5
)
=
1
\frac{ABsinX}{BDsin(\angle3 + \angle5)}=1
BDsin(∠3+∠5)ABsinX=1
现在就有四个无关的条件方程了,当然上面只是一些关系,角度都是表示的平差后的值,我们正规的列出条件方程应该是这样的
v
1
+
v
2
+
v
3
+
w
a
=
0
v
4
+
v
5
+
v
6
+
w
b
=
0
s
i
n
(
L
^
5
)
s
i
n
(
L
^
2
+
L
^
4
)
s
i
n
X
^
s
i
n
(
L
^
1
−
X
^
)
s
i
n
(
L
^
4
)
s
i
n
(
L
^
3
+
L
^
5
)
=
1
A
B
s
i
n
X
^
B
D
s
i
n
(
L
^
3
+
L
^
5
)
=
1
v_{1}+ v_{2}+ v_{3}+w_{a}=0\\ v_{4}+ v_{5}+ v_{6}+w_{b}=0\\ \frac{sin(\hat L_{5})sin(\hat L_{2} + \hat L_{4})sin\hat X}{sin(\hat L_{1} - \hat X)sin(\hat L_{4})sin(\hat L_{3} + \hat L_{5})}=1\\ \frac{ABsin\hat X}{BDsin(\hat L_{3} + \hat L_{5})}=1
v1+v2+v3+wa=0v4+v5+v6+wb=0sin(L^1−X^)sin(L^4)sin(L^3+L^5)sin(L^5)sin(L^2+L^4)sinX^=1BDsin(L^3+L^5)ABsinX^=1
这里附有参数的条件平差模型为
A
V
+
B
X
^
+
W
=
0
AV+B\hat X+W=0
AV+BX^+W=0,目的仍然是使得改正数的平方和最小。根据拉格朗日乘子法,新的代价函数为
Φ
=
V
T
P
V
−
2
K
T
(
A
V
+
B
x
^
+
W
)
Φ=V^{T}PV−2K^{T}(AV+B\hat x+W)
Φ=VTPV−2KT(AV+Bx^+W)。同样求导使得
∂
(
Φ
)
∂
(
V
)
=
0
\frac{\partial(\Phi)}{\partial(V)}=0
∂(V)∂(Φ)=0和
∂
(
Φ
)
∂
(
x
^
)
=
0
\frac{\partial(\Phi)}{\partial(\hat x)}=0
∂(x^)∂(Φ)=0,最后得到解为
x
^
=
−
N
B
B
−
1
B
T
N
A
A
−
1
W
V
=
−
P
−
1
A
T
N
A
A
−
1
(
B
x
^
+
W
)
\hat x=−N^{−1}_{BB}B^{T}N^{−1}_{AA}W\\ V=−P^{-1}A^{T}N^{−1}_{AA}(B\hat x+W)
x^=−NBB−1BTNAA−1WV=−P−1ATNAA−1(Bx^+W)
这里
N
A
A
=
A
Q
A
T
N_{AA}=AQA^{T}
NAA=AQAT,
N
B
B
=
B
T
N
A
A
−
1
B
N_{BB}=B^{T}N^{−1}_{AA}B
NBB=BTNAA−1B,详细的推导过程可以参考这里。现在可以开始列矩阵了。其中需要将非线性条件(极条件和固定边条件)线性化,这里用到的就是使用泰勒公式展开,这里大概取个
X
0
X^{0}
X0是30度,具体怎么展开不写了,百度一堆
s
i
n
x
sinx
sinx的泰勒展开式。现在我们根据展开后的式子就可以开始列出条件方程的
A
A
A阵、
B
B
B阵和
W
W
W阵如下,由于是等精度观测,
P
=
I
P=I
P=I
A
=
[
1
1
1
0
0
0
0
0
0
1
1
1
1.732
0.577
−
0.577
1.115
−
1.115
0
0
0
0.577
0
0.577
0
]
B
=
[
0
0
−
3.464
1.732
]
W
=
[
9
−
8
5.196
−
6.051
]
A= \left[ \begin{matrix} 1 & 1 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1\\ 1.732 & 0.577 & -0.577 & 1.115 & -1.115 & 0\\ 0 & 0 & 0.577 & 0 & 0.577 & 0 \end{matrix} \right]\\ B= \left[ \begin{matrix} 0\\ 0\\ -3.464\\ 1.732 \end{matrix} \right] W= \left[ \begin{matrix} 9\\ -8\\ 5.196\\ -6.051 \end{matrix} \right]
A=⎣⎢⎢⎡101.7320100.577010−0.5770.577011.115001−1.1150.5770100⎦⎥⎥⎤B=⎣⎢⎢⎡00−3.4641.732⎦⎥⎥⎤W=⎣⎢⎢⎡9−85.196−6.051⎦⎥⎥⎤
现在我们就可以代入公式计算答案了
A = np.array([[1,1,1,0,0,0],
[0,0,0,1,1,1],
[1.732,0.577,-0.577,1.155,-1.155,0],
[0,0,0.577,0,0.577,0]])
B = np.array([0,0,-3.464,1.732]).T.reshape(4,1)
W = np.array([9,-8,5.196,-6.051]).T.reshape(4,1)
P = np.eye(6)
NAA = np.dot(np.dot(A, np.linalg.inv(P)), A.T)
NBB = np.dot(np.dot(B.T, np.linalg.inv(NAA)), B)
x_h = -np.dot(np.dot(np.dot(np.linalg.inv(NBB), B.T), np.linalg.inv(NAA)), W)
V = -np.dot(np.dot(np.dot(np.linalg.inv(P), A.T), np.linalg.inv(NAA)), (np.dot(B, x_h)+W))
print(x_h)
print(V)
得到的结果为
[[5.40626195]]
[[ 2.40314236]
[-5.70157118]
[-5.70157118]
[ 8.07214805]
[-0.03958256]
[-0.03256549]]
现在我们就得到
x
^
=
30
°
+
5.4
″
=
30
°
00
′
05.4
″
\hat x=30°+5.4″=30°00′05.4″
x^=30°+5.4″=30°00′05.4″,对应的观测值的平差值就按照
L
^
=
L
0
+
v
\hat L = L^{0}+v
L^=L0+v依次算出就行了
这里我们四种平差模型都复习了一遍,没有刚开始感觉的那么难,其实平差就是一个优化的问题,核心准则是最小二乘法。可以简单的看成找到平差值和观测值之间的改正数,再使得改正数的平方和最小,这时候再改正每一个观测值。有点像我们开始写的线性回归,只是做起来麻烦些而已