对于线性回归问题,我们说明了利用梯度下降和正规方程进行处理的原理,其中最重要的就是最小二乘的原理,而在测绘科学中,测量平差的法方程就对应着线性回归的正规方程。平差利用最小二乘,却也有自己的不同,下面我们通过一道水准网间接平差例题来理解间接平差,间接平差就是在一个平差问题中,独立参数个数等于必要观测数
t
t
t 时,将每个观测值表达成这t个参数的函数,组成观测方程的平差方法
如题
如果没有学过平差,我们很容易从图中看出,
h
P
1
=
h
A
+
h
1
h_{P1}=h_A+h_1
h P 1 = h A + h 1 ,
h
P
2
=
h
C
+
h
3
h_{P2}=h_C+h_3
h P 2 = h C + h 3 ,那不就太简单了嘛,
P
1
P_1
P 1 就等于
12.003
12.003
1 2 . 0 0 3 m,
P
2
P_2
P 2 就等于
12.511
12.511
1 2 . 5 1 1 m。真这么简单就好了,由于观测的误差,答案肯定是接近于此但是不为此,我们可以看到根据图形和观测值,得到三个‘正确’的公式:
h
P
1
=
h
A
+
h
1
=
12.003
h_{P1}=h_A+h_1=12.003
h P 1 = h A + h 1 = 1 2 . 0 0 3 m,
h
P
1
=
h
B
+
h
4
=
12.005
h_{P1}=h_B+h_4=12.005
h P 1 = h B + h 4 = 1 2 . 0 0 5 m,
h
P
1
=
h
C
+
h
3
−
h
2
=
12.010
h_{P1}=h_C+h_3-h_2=12.010
h P 1 = h C + h 3 − h 2 = 1 2 . 0 1 0 m。同一个
P
1
P_1
P 1 的高度通过三个方向计算却不同,如何确定
P
1
P_1
P 1 最‘正确’的高度,这就需要平差来解决了。不过至少我们能得到这样一个答案,这里至少需要两次观测才能求出
P
1
P_1
P 1 和
P
2
P_2
P 2 ,即要测得
h
3
h_3
h 3 和
h
1
h_1
h 1 或者
h
3
h_3
h 3 和
h
4
h_4
h 4 。那我们就可以得到这道题的必要观测数
t
=
2
t=2
t = 2 。我们设
P
1
P_1
P 1 和
P
2
P_2
P 2 的真实高度为
X
1
^
\hat{X_1}
X 1 ^ 和
X
2
^
\hat{X_2}
X 2 ^ ,前面我们也算得了P1和P2的近似高度,我们就设近似高度为
X
1
X_1
X 1 和
X
2
X_2
X 2 。下面我们设每一条边的改正数设为
v
v
v ,观测值为
L
L
L ,真实值
L
^
=
L
+
v
\hat{L}=L+v
L ^ = L + v ,而对于线性表达式,
L
^
=
B
h
^
+
d
(
y
=
a
x
+
b
)
\hat{L}=B\hat{h}+d(y=ax+b)
L ^ = B h ^ + d ( y = a x + b ) 。这里我们就可以得到
L
+
v
=
B
X
^
+
d
L+v=B\hat{X}+d
L + v = B X ^ + d ,结合
X
^
=
X
+
x
^
\hat{X}=X+\hat{x}
X ^ = X + x ^ 和
l
=
L
^
−
L
=
L
−
(
B
X
+
d
)
l=\hat{L}-L=L-(BX+d)
l = L ^ − L = L − ( B X + d ) 我们就可以得到误差方程
L
+
v
=
B
X
^
+
d
B
X
+
d
+
l
+
v
=
B
X
+
B
x
^
+
d
v
=
B
x
^
−
(
L
^
−
L
)
v
=
B
x
^
−
l
\begin{aligned} L + v &= B\hat{X} + d\\ BX + d + l + v& = BX + B\hat{x} + d\\ v &= B\hat{x} - (\hat{L} - L) \\ v &= B\hat{x} - l \end{aligned}
L + v B X + d + l + v v v = B X ^ + d = B X + B x ^ + d = B x ^ − ( L ^ − L ) = B x ^ − l 现在我们列出误差方程
v
1
=
x
1
^
−
(
h
1
−
(
X
1
−
h
A
)
)
v
2
=
−
x
1
^
+
x
2
^
−
(
h
2
−
(
X
2
−
X
1
)
)
v
3
=
x
2
^
−
(
h
3
−
(
X
2
−
h
C
)
)
v
4
=
x
1
^
−
(
h
4
−
(
X
1
−
h
B
)
)
\begin{aligned} &v_1 = \hat{x_1} - (h_1 - (X_1 - h_A))\\ &v_2 = -\hat{x_1} + \hat{x_2} - (h_2 - (X_2 - X_1))\\ &v_3 = \hat{x_2} - (h_3 - (X_2 - h_C))\\ &v_4 = \hat{x_1} - (h_4 - (X_1 - h_B)) \end{aligned}
v 1 = x 1 ^ − ( h 1 − ( X 1 − h A ) ) v 2 = − x 1 ^ + x 2 ^ − ( h 2 − ( X 2 − X 1 ) ) v 3 = x 2 ^ − ( h 3 − ( X 2 − h C ) ) v 4 = x 1 ^ − ( h 4 − ( X 1 − h B ) )
把数字带进去,化简,就得到如下的式子(这里的l的单位是毫米)
v
1
=
x
1
^
−
0
v
2
=
−
x
1
^
+
x
2
^
−
(
−
7
)
v
3
=
x
2
^
−
0
v
4
=
x
1
^
−
2
\begin{aligned} &v_1 = \hat{x_1} - 0\\ &v_2 = -\hat{x_1} + \hat{x_2} - (-7)\\ &v_3 = \hat{x_2} - 0\\ &v_4 = \hat{x_1} - 2\\ \end{aligned}
v 1 = x 1 ^ − 0 v 2 = − x 1 ^ + x 2 ^ − ( − 7 ) v 3 = x 2 ^ − 0 v 4 = x 1 ^ − 2
对比
v
=
B
x
^
−
l
v = B\hat{x} - l
v = B x ^ − l 的形式,我们可以写出
B
B
B 和
l
l
l 矩阵,根据题目的边长的长度,我们以最长的
2
2
2 km为单位权,就可以列出权重的矩阵
P
P
P
B
=
[
1
0
−
1
1
0
1
1
0
]
l
=
[
0
−
7
0
2
]
P
=
[
2
0
0
0
0
1
0
0
0
0
1
0
0
0
0
2
]
B= \left[ \begin{matrix} 1 & 0\\ -1 & 1\\ 0 & 1\\ 1& 0 \end{matrix} \right] l= \left[ \begin{matrix} 0\\ -7\\ 0\\ 2 \end{matrix} \right] P= \left[ \begin{matrix} 2 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 2 \end{matrix} \right]
B = ⎣ ⎢ ⎢ ⎡ 1 − 1 0 1 0 1 1 0 ⎦ ⎥ ⎥ ⎤ l = ⎣ ⎢ ⎢ ⎡ 0 − 7 0 2 ⎦ ⎥ ⎥ ⎤ P = ⎣ ⎢ ⎢ ⎡ 2 0 0 0 0 1 0 0 0 0 1 0 0 0 0 2 ⎦ ⎥ ⎥ ⎤
平差中的准则是要是
V
T
P
V
=
m
i
n
V^TPV=min
V T P V = m i n (也就是
[
v
v
]
=
∑
v
2
[vv]=∑v^2
[ v v ] = ∑ v 2 要最小),这是进行最小二乘平差计算的一个基本原则,是求解不定线性方程组的一个附加条件。当VTPV越小,平差的效果越好。其中我们也得到了
v
=
B
x
^
−
l
v = B\hat{x} - l
v = B x ^ − l 。现在就可以按列正规方程 的方法(
V
T
P
V
V^TPV
V T P V 关于
x
^
\hat{x}
x ^ 求导,当导数为零时,
V
T
P
V
V^TPV
V T P V 最小),求得
x
^
=
(
B
T
P
B
)
−
1
B
T
P
l
\hat{x}=(B^TPB)^{-1}B^TPl
x ^ = ( B T P B ) − 1 B T P l ,下面是推导
d
(
V
T
P
V
)
/
d
(
x
^
)
=
0
2
V
T
P
(
d
(
V
)
/
d
(
x
^
)
)
=
0
2
V
T
P
(
d
(
B
x
^
−
l
)
/
d
(
x
^
)
)
=
0
V
T
P
B
=
0
B
T
P
V
=
0
B
T
P
(
B
x
^
−
l
)
=
0
B
T
P
B
x
^
−
B
T
P
l
=
0
x
^
=
(
B
T
P
B
)
−
1
B
T
P
l
\begin{aligned} d(V^TPV) / d(\hat{x})& = 0\\ 2V^TP(d(V) / d(\hat{x}))& = 0 \\ 2V^TP(d(B\hat{x} - l) / d(\hat{x}))& = 0\\ V^TPB& = 0\\ B^TPV& = 0\\ B^TP(B\hat{x} - l)& = 0 \\ B^TPB\hat{x} - B^TPl& = 0\\ \hat{x}& = (B^TPB)^{-1}B^TPl \end{aligned}
d ( V T P V ) / d ( x ^ ) 2 V T P ( d ( V ) / d ( x ^ ) ) 2 V T P ( d ( B x ^ − l ) / d ( x ^ ) ) V T P B B T P V B T P ( B x ^ − l ) B T P B x ^ − B T P l x ^ = 0 = 0 = 0 = 0 = 0 = 0 = 0 = ( B T P B ) − 1 B T P l 现在就可以很方便地计算了,我们用python-numpy计算一下
import numpy as np
B = np. array( [ [ 1 , 0 ] , [ - 1 , 1 ] , [ 0 , 1 ] , [ 1 , 0 ] ] )
l = np. array( [ [ 0 ] , [ - 7 ] , [ 0 ] , [ 2 ] ] )
P = np. array( [ [ 2 , 0 , 0 , 0 ] , [ 0 , 1 , 0 , 0 ] , [ 0 , 0 , 1 , 0 ] , [ 0 , 0 , 0 , 2 ] ] )
print ( B. shape)
print ( l. shape)
print ( P. shape)
输出
B
B
B 、
l
l
l 和
P
P
P 矩阵的维度
(4, 2)
(4, 1)
(4, 4)
现在计算,在python-numpy中,np.dot()
是矩阵相乘,np.linalg.inv()
是求逆矩阵,现在输出看看
x_h = np. zeros( ( 2 , 1 ) )
x_h = np. dot( np. dot( np. dot( np. linalg. inv( np. dot( np. dot( B. T, P) , B) ) , B. T) , P) , l)
print ( x_h)
结果是这样的,这里的单位是毫米(mm)
[[ 1.66666667]
[-2.66666667]]
所以P1和P2点的高程平差值应该为
X
^
=
X
+
x
^
\hat{X}=X+\hat{x}
X ^ = X + x ^ ,即
P
1
P_1
P 1 的高为
12.0047
12.0047
1 2 . 0 0 4 7 m和
P
2
P_2
P 2 的高为
12.5083
12.5083
1 2 . 5 0 8 3 m。这个时候吧
x
^
\hat{x}
x ^ 代回去也可以得到所有的
v
v
v 的值,就可以求出所有观测量的平差值。我们可以算一下此时的所有观测量的平差值
h
1
^
=
h
1
+
x
1
^
−
0
=
1.0047
h
2
^
=
h
2
+
−
x
1
^
+
x
2
^
−
(
−
7
)
=
0.5037
h
3
^
=
h
3
+
x
2
^
−
0
=
0.5003
h
4
^
=
h
4
+
x
1
^
−
2
=
0.5047
\begin{aligned} &\hat{h_1} = h_1 + \hat{x_1} - 0 = 1.0047\\ &\hat{h_2} = h_2 + -\hat{x_1} + \hat{x_2} - (-7) = 0.5037\\ &\hat{h_3} = h_3 + \hat{x_2} - 0 = 0.5003\\ &\hat{h_4} = h_4 + \hat{x_1} - 2 = 0.5047 \end{aligned}
h 1 ^ = h 1 + x 1 ^ − 0 = 1 . 0 0 4 7 h 2 ^ = h 2 + − x 1 ^ + x 2 ^ − ( − 7 ) = 0 . 5 0 3 7 h 3 ^ = h 3 + x 2 ^ − 0 = 0 . 5 0 0 3 h 4 ^ = h 4 + x 1 ^ − 2 = 0 . 5 0 4 7
这里我们再算一下
h
P
1
h_{P1}
h P 1 ,就发现三个方向算的答案都基本可以对的上了
h
P
1
=
h
A
+
h
1
^
=
12.0047
h_{P1}=h_A+\hat{h_1}=12.0047
h P 1 = h A + h 1 ^ = 1 2 . 0 0 4 7 m,
h
P
1
=
h
B
+
h
4
^
=
12.0047
h_{P1}=h_B+\hat{h_4}=12.0047
h P 1 = h B + h 4 ^ = 1 2 . 0 0 4 7 m,
h
P
1
=
h
C
+
h
3
^
−
h
2
^
=
12.0046
h_{P1}=h_C+\hat{h_3 }-\hat{h_2}=12.0046
h P 1 = h C + h 3 ^ − h 2 ^ = 1 2 . 0 0 4 6 m 结果里面有省略一些地小数,所以没有完全的相等