单像空间后方交会
1. 概要
单像后方交会的基本思想是至少利用三个已知地面控制点的坐标 \(A(X_A,Y_A,Z_A)\)，\(B(X_B,Y_B,Z_B)\)，\(C(X_C,Y_C,Z_C)\) ，与影像上对应的三个像点的影像坐标 \(a(x_a,y_a)\)，\(b(x_b,y_b)\)，\(c(x_c,y_c)\) ，再根据共线方程，反求该像片的外方位元素 \(X_s,Y_s,Z_s,\varphi,\omega,\kappa\)。这种解法是以单张像片为基础，所以也称做单像空间后方交会。
2. 推导
2.1. 基本关系式
共线方程： \[ \begin{array}{c} x=-f\frac{a_1\left( X-X_s \right) +b_1\left( Y-Y_s \right) +c_1\left( Z-Z_s \right)}{a_3\left( X-X_s \right) +b_3\left( Y-Y_s \right) +c_3\left( Z-Z_s \right)} \\ y=-f\frac{a_2\left( X-X_s \right) +b_2\left( Y-Y_s \right) +c_2\left( Z-Z_s \right)}{a_3\left( X-X_s \right) +b_3\left( Y-Y_s \right) +c_3\left( Z-Z_s \right)} \\ \end{array} \]
为了方便计算，需要把非线性的共线方程线性化。\(x,y\)
是关于六个外方位元素的多元函数，进行一次泰勒展开： \[\label{共线方程}
\begin{array}{c}
x = \tilde{x} + \frac{\partial x}{\partial X_s}dX_s + \frac{\partial
x}{\partial Y_s}dY_s+\frac{\partial x}{\partial Z_s}dZ_s+\frac{\partial
x}{\partial \varphi}d\varphi +\frac{\partial x}{\partial \omega}d\omega
+\frac{\partial x}{\partial \kappa}d\kappa \\
y = \tilde{y} +\frac{\partial y}{\partial X_s}dX_s+\frac{\partial
y}{\partial Y_s}dY_s+\frac{\partial y}{\partial Z_s}dZ_s+\frac{\partial
y}{\partial \varphi}d\varphi +\frac{\partial y}{\partial \omega}d\omega
+\frac{\partial y}{\partial \kappa}d\kappa
\end{array}
\] \(\tilde{x},\tilde{y}\)
是函数的近似值，是把外方位元素的初值 \(X_{s0},Y_{s0},Z_{s0},\varphi_0,\omega_0,\kappa_0\)
带入共线方程中得到的结果；
\(dX_s,dY_x,dZ_s,d\varphi_,d\omega,d\kappa\)
是对应外方位元素近似值的改正数（也就是增量）；
\(\frac{\partial x}{\partial
X_s},\dots,\frac{\partial y}{\partial \kappa}\)
为偏导数，是外方位元素改正数的系数。
为了方便表示，式 \(\eqref{共线方程}\) 记为： \[\label{共线方程简化} \begin{array}{c} x = \tilde{x} + a_{11}\ dX_s + a_{12}\ dY_s+a_{13}\ dZ_s+a_{14}\ d\varphi +a_{15}\ d\omega +a_{16}\ d\kappa= \tilde{x} + dx \\ y = \tilde{y} +a_{21}\ dX_s+a_{22}\ dY_s+a_{23}\ dZ_s+a_{24}\ d\varphi +a_{25}\ d\omega +a_{26}\ d\kappa=\tilde{y} + dy \end{array} \] \(dx,dy\) 表示近似值的改正值。
因为上述近似是一阶的，结果较为粗略，所以需要不断迭代趋近，直到改正值小于某一阈值。
对于一个控制点就能由共线方程列写两个式子，所以一个像片中有三个点就能列出六个方程，求出六个外方位元素的改正值。
2.2. 误差方程与法方程
为了提高解算外方位元素的精度，常有多余观测方程。在空间后方交会中，通常是在像片的四个角上选取四个或者更多的地面控制点，因此要用最小二乘法平差计算。
计算中，通常将控制点的地面坐标视为真值，而把相应的像点坐标视为观测值，加入相应的改正数
\(v_x,v_y\)，按照
观测值 + 观测值改正数 = 近似值 + 近似值改正数
的改正原则，得： \[
\begin{array}{c}
x+v_x=\tilde{x}+dx \\
y+v_y=\tilde{y}+dy
\end{array}
\Longrightarrow
\begin{array}{c}
v_x=\tilde{x}+dx-x \\
v_y=\tilde{y}+dy-y
\end{array}
\] 这里 \(v_x,v_y\)
的意义就比较明显了，代表 近似值经过改正后与真实值之间的误差。
结合式 \(\eqref{共线方程简化}\)，我们可以用矩阵形式表示这个方程组 \[ V=AX-l \] 其中 \[ \begin{equation*} \begin{aligned} V&=\left[v_x, v_y\right]^T \in\mathbb{R}^{2\times 1} \\\\ A&=\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & a_{15} & a_{16} \\ a_{21} & a_{22} & a_{23} & a_{24} & a_{25} & a_{26} \end{bmatrix} \in\mathbb{R}^{2\times 6} \\\\ X&=\left[dX_s,dY_s,dZ_s,d\varphi,d\omega,d\kappa\right]^T \in\mathbb{R}^{6\times 1} \\\\ l&=\left[l_x,l_y\right]^T=\left[x-\tilde{x}, y-\tilde{y}\right]^T \in\mathbb{R}^{2\times 1} \end{aligned} \end{equation*} \] 上面只是一个控制点的形式，当存在 \(n\) 个控制点的时候，则有： \[\label{误差方程} V=AX-L \] 其中 \[ \begin{equation*} \begin{aligned} V&=\left[V_1^T , V_2^T,\dots,V_n^T\right]^T \in \mathbb{R}^{2n\times 1} \\\\ A&=\left[A_1^T,A_2^T,\dots,A_n^T \right] \in\mathbb{R}^{2n\times 6} \\\\ X&=\left[dX_s,dY_s,dZ_s,d\varphi,d\omega,d\kappa\right]^T \in\mathbb{R}^{6\times 1} \\\\ l&=\left[l_1^T,l_2^T,\dots,l_n^T \right]^T \in\mathbb{R}^{2n\times 1} \end{aligned} \end{equation*} \] 注意矩阵的大小，写代码的时候要用到。
根据最小二乘间接平差原理，可列出法方程： \[ A^T P A X=A^T L \] 其中，\(P\) 为像点观测值的权值矩阵，表示观测值测量的相对精度，一般认为是等精度测量，故取 \(P\) 为单位阵，因此可解得： \[\label{封闭解} X=(A^T A)^{-1}A^T L \] 从而求得外方位元素近似值的改正数 \(dX_s,dY_s,dZ_s,d\varphi,d\omega,d\kappa\) 。
一些解释
这里的计算过程其实就是很常见的线性模型，最小化误差即 \[\nonumber X=arg\ \underset{X}{min}\|V^TV\|_{2}^{2}\] 这是一个经典的凸优化问题，求偏导，令偏导等于 0 即可。
根据上面的内容，我们可以得到第 \(i\) 次计算的外方位元素的改正值 \(X_i\)，直到第 \(k\) 次收敛的时，得到最后的外方位元素： \[ \begin{equation*} \begin{aligned} X_s &= X_{s0}+dX_{s1}+dX_{s2}+\dots+dX_{sk} \\ Y_s &= Y_{s0}+dY_{s1}+dY_{s2}+\dots+dY_{sk} \\ Z_s &= Z_{s0}+dZ_{s1}+dZ_{s2}+\dots+dZ_{sk} \\ \varphi &= \varphi_0 + d\varphi_1 + d\varphi_2+\dots+d\varphi_k \\ \omega &= \omega_0 + d\omega_1+d\omega_2 + \dots+d\omega_k \\ \kappa &= \kappa_0+d\kappa_1+d\kappa_2+\dots+d\kappa_k \end{aligned} \end{equation*} \]
即 \[ X = X_0 + X_1 + X_2 +\dots+X_k \]
观察式 \(\eqref{封闭解}\)
可知，目前仅 \(A\) 未知，下面对 \(A\) 进行求解。 求解 \(A\) 为书写方便，令共线方程中的分子、分母如下： \[\begin{equation}\label{共线方程分子分母}
\begin{aligned}
\overline{X} = a_1\left( X-X_s \right) +b_1\left( Y-Y_s \right)
+c_1\left( Z-Z_s \right) \\
\overline{Y} = a_2\left( X-X_s \right) +b_2\left( Y-Y_s \right)
+c_2\left( Z-Z_s \right) \\
\overline{Z} = a_3\left( X-X_s \right) +b_3\left( Y-Y_s \right)
+c_3\left( Z-Z_s \right) \\
\end{aligned}
\end{equation}
\] 由 \(\eqref{共线方程}\) 和 \(\eqref{共线方程简化}\) 可知 \[
\begin{equation*}
\begin{aligned}
a_{11} = \frac{\partial x}{\partial X_s} &= \frac{\partial
\left(-f\frac{\overline{X}}{\overline{Z}}\right)}{\partial X_s} &
分式求导\\
&=-f \frac{\frac{\partial \overline{X}}{\partial X_s}\overline{Z} -
\overline{X}\frac{\partial \overline{Z}}{\partial X_s}}{\overline{Z}^2}
&x=-f\frac{\overline{X}}{\overline{Z}}\\
&=-f\left(-\frac{a_1}{\overline{Z}} - \frac{1}{f\overline{Z}}x
a_3\right) \\
&=\frac{1}{\overline{Z}}\left(a_1 f + a_3 x\right)
\end{aligned}
\end{equation*}
\] 同理可得 \[
\begin{align}
&
\begin{cases}\label{对XYZ偏导}
a_{11}=\frac{\partial x}{\partial X_{s}}=\frac{1}{\overline{Z}}(a_{1}f +
a_{3}x)\\\\
a_{12}=\frac{\partial x}{\partial Y_{s}}=\frac{1}{\overline{Z}}(b_{1}f +
b_{3}x)\\\\
a_{13}=\frac{\partial x}{\partial Z_{s}}=\frac{1}{\overline{Z}}(c_{1}f +
c_{3}x)\\\\
a_{21}=\frac{\partial y}{\partial X_{s}}=\frac{1}{\overline{Z}}(a_{2}f +
a_{3}y)\\\\
a_{22}=\frac{\partial y}{\partial Y_{s}}=\frac{1}{\overline{Z}}(b_{2}f +
b_{3}y)\\\\
a_{23}=\frac{\partial y}{\partial Z_{s}}=\frac{1}{\overline{Z}}(c_{2}f +
c_{3}y)
\end{cases}
\\\nonumber \\
&
\begin{cases}\label{对角度偏导}
a_{14}=\frac{\partial x}{\partial
\varphi}=-\frac{f}{(\overline{Z})^{2}}(\frac{\partial
\overline{X}}{\partial \varphi}\overline{Z}-\frac{\partial
\overline{Z}}{\partial \varphi}\overline{X})\\\\
a_{15}=\frac{\partial x}{\partial
\omega}=-\frac{f}{(\overline{Z})^{2}}(\frac{\partial
\overline{X}}{\partial \omega}\overline{Z}-\frac{\partial
\overline{Z}}{\partial \omega}\overline{X})\\\\
a_{16}=\frac{\partial x}{\partial
\kappa}=-\frac{f}{(\overline{Z})^{2}}(\frac{\partial
\overline{X}}{\partial \kappa}\overline{Z}-\frac{\partial
\overline{Z}}{\partial \kappa}\overline{X})\\\\
a_{24}=\frac{\partial y}{\partial
\varphi}=-\frac{f}{(\overline{Z})^{2}}(\frac{\partial
\overline{Y}}{\partial \varphi}\overline{Z}-\frac{\partial
\overline{Z}}{\partial \varphi}\overline{Y})\\\\
a_{25}=\frac{\partial y}{\partial
\omega}=-\frac{f}{(\overline{Z})^{2}}(\frac{\partial
\overline{Y}}{\partial \omega}\overline{Z}-\frac{\partial
\overline{Z}}{\partial \omega}\overline{Y})\\\\
a_{26}=\frac{\partial y}{\partial
\kappa}=-\frac{f}{(\overline{Z})^{2}}(\frac{\partial
\overline{Y}}{\partial \kappa}\overline{Z}-\frac{\partial
\overline{Z}}{\partial \kappa}\overline{Y})
\end{cases}
\end{align}
\] 由式 \(\eqref{共线方程分子分母}\)
可知 \[
\begin{aligned}
\left[ \begin{array}{c}
\overline{X} \\
\overline{Y} \\
\overline{Z} \\
\end{array} \right]
&=\left[\begin{array}{c}
a_1 & a_2 & a_3 \\
b_1 & b_2 & b_3 \\
c_1 & c_2 & c_3 \\
\end{array}\right]
\left[\begin{array}{c}
X-X_s\\
Y-Y_s\\
Z-Z_s
\end{array}\right]
=R^T\left[\begin{array}{c}
X-X_s\\
Y-Y_s\\
Z-Z_s
\end{array}\right] \\
&=R_{\kappa}^T R_{\omega}^T R_{\varphi}^T \left[\begin{array}{c}
X-X_s\\
Y-Y_s\\
Z-Z_s
\end{array}\right]
=R_{\kappa}^{-1} R_{\omega}^{-1} R_{\varphi}^{-1} \left[\begin{array}{c}
X-X_s\\
Y-Y_s\\
Z-Z_s
\end{array}\right]
\end{aligned}
\] 所以 \[
\begin{equation}
\begin{aligned}
\frac{\partial}{\partial \varphi}
\begin{bmatrix}
\overline{X} \\
\overline{Y} \\
\overline{Z}
\end{bmatrix}
&= R_{\kappa}^{-1} R_{\omega}^{-1} \frac{\partial
R_{\varphi}^{-1}}{\partial \varphi}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&= R_{\kappa}^{-1} R_{\omega}^{-1} R_{\varphi}^{-1} R_{\varphi}
\frac{\partial R_{\varphi}^{-1}}{\partial \varphi}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&= R^{-1} R_{\varphi} \frac{\partial R_{\varphi}^{-1}}{\partial
\varphi}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}
\end{aligned}
\end{equation}
\] 而 \[
R_{\varphi}^{-1} = R_{\varphi}^{\mathrm{T}}
=\begin{bmatrix}
\cos\varphi & 0 & \sin\varphi \\
0 & 1 & 0 \\
-\sin\varphi & 0 & \cos\varphi
\end{bmatrix}
\] 则 \[
R_{\varphi} \frac{\partial R_{\varphi}^{-1}}{\partial \varphi}
=\begin{bmatrix}
\cos\varphi & 0 & -\sin\varphi \\
0 & 1 & 0 \\
\sin\varphi & 0 & \cos\varphi
\end{bmatrix}
\begin{bmatrix}
-\sin\varphi & 0 & \cos\varphi \\
0 & 0 & 0 \\
-\cos\varphi & 0 & -\sin\varphi
\end{bmatrix}
=\begin{bmatrix}
0 & 0 & 1 \\
0 & 0 & 0 \\
-1 & 0 & 0
\end{bmatrix}
\] 代入上式，得： \[
\begin{equation}\label{varphi偏导}
\begin{aligned}
\frac{\partial}{\partial \varphi}
\begin{bmatrix}
\overline{X} \\
\overline{Y} \\
\overline{Z}
\end{bmatrix}
&=\begin{bmatrix}
a_{1} & b_{1} & c_{1} \\
a_{2} & b_{2} & c_{2} \\
a_{3} & b_{3} & c_{3}
\end{bmatrix}
\begin{bmatrix}
0 & 0 & 1 \\
0 & 0 & 0 \\
-1 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&=\begin{bmatrix}
a_{1} & b_{1} & c_{1} \\
a_{2} & b_{2} & c_{2} \\
a_{3} & b_{3} & c_{3}
\end{bmatrix}
\begin{bmatrix}
0 & 0 & 1 \\
0 & 0 & 0 \\
-1 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
a_{1} & a_{2} & a_{3} \\
b_{1} & b_{2} & b_{3} \\
c_{1} & c_{2} & c_{3}
\end{bmatrix}
\begin{bmatrix}
\overline{X} \\
\overline{Y} \\
\overline{Z}
\end{bmatrix}\\
&=\begin{bmatrix}
0 & -b_{3} & b_{2} \\
b_{3} & 0 & -b_{1} \\
-b_{2} & b_{1} & 0
\end{bmatrix}
\begin{bmatrix}
\overline{X} \\
\overline{Y} \\
\overline{Z}
\end{bmatrix}
\end{aligned}
\end{equation}
\] 按相仿的方法，得： \[
\begin{equation}\label{omega偏导}
\begin{aligned}
\frac{\partial}{\partial \omega}
\begin{bmatrix}
\overline{X} \\
\overline{Y} \\
\overline{Z}
\end{bmatrix}
&= R_{\kappa}^{-1} \frac{\partial R_{\omega}^{-1}}{\partial
\omega} R_{\varphi}^{-1}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&= R_{\kappa}^{-1} \frac{\partial R_{\omega}^{-1}}{\partial
\omega} R_{\omega} R_{\kappa} R_{\kappa}^{-1} R_{\omega}^{-1}
R_{\varphi}^{-1}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&=R_{\kappa}^{-1}
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 1 \\
0 & -1 & 0
\end{bmatrix}
R_{\kappa}R^{-1}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&=\begin{bmatrix}
\overline{Z}\sin\kappa \\
\overline{Z}\cos\kappa \\
-\overline{X}\sin\kappa - Y\cos\kappa
\end{bmatrix}
\end{aligned}
\end{equation}
\] \[
\begin{equation}\label{kappa偏导}
\begin{aligned}
\frac{\partial}{\partial \kappa}
\begin{bmatrix}
\overline{X} \\
\overline{Y} \\
\overline{Z}
\end{bmatrix}
&=\frac{\partial R_{\kappa}^{-1}}{\partial
\kappa}R_{\kappa}R_{\kappa}^{-1}R_{\omega}^{-1}R_{\varphi}^{-1}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&=\begin{bmatrix}
0 & 1 & 0 \\
-1 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
a_{1} & b_{1} & c_{1} \\
a_{2} & b_{2} & c_{2} \\
a_{3} & b_{3} & c_{3}
\end{bmatrix}
\begin{bmatrix}
X - X_{s} \\
Y - Y_{s} \\
Z - Z_{s}
\end{bmatrix}\\
&=\begin{bmatrix}
\overline{Y} \\
-\overline{X} \\
0
\end{bmatrix}
\end{aligned}
\end{equation}
\] 将 \(\eqref{varphi偏导} \eqref{omega偏导}
\eqref{kappa偏导}\) 带入 \(\eqref{对角度偏导}\)，整理得： \[
\begin{equation}\label{对角度偏导化简}
\begin{aligned}
a_{14} &= y\sin\omega - \left[\frac{x}{f}(x\cos\kappa - y\sin\kappa)
+ f\cos\kappa\right]\cos\omega\\
a_{15} &= -f\sin\kappa - \frac{x}{f}(x\sin\kappa + y\cos\kappa)\\
a_{16} &= y\\
a_{24} &= -x\sin\omega - \left[\frac{y}{f}(x\cos\kappa -
y\sin\kappa) - f\sin\kappa\right]\cos\omega\\
a_{25} &= -f\cos\kappa - \frac{y}{f}(x\sin\kappa + y\cos\kappa)\\
a_{26} &= -x
\end{aligned}
\end{equation}
\]
在竖直摄影情况下，角元素都是小角\((\lt 3°)\)，可用 \(\varphi=\omega=\kappa=0\) 及 \(Z-Z_s=-H\) 代替，得到各系数的近似值： \[ \begin{equation} \begin{aligned} a_{11} &= -\frac{f}{H} & a_{12} &= 0 & a_{13} &= -\frac{x}{H}\\ a_{14} &= -f\left(1 + \frac{x^2}{f^2}\right) & a_{15} &= -\frac{xy}{f} & a_{16} &= y\\ a_{21} &= 0 & a_{22} &= -\frac{f}{H} & a_{23} &= -\frac{y}{H}\\ a_{24} &= -\frac{xy}{f} & a_{25} &= -f\left(1 + \frac{y^2}{f^2}\right) & a_{26} &= -x \end{aligned} \end{equation} \]
3. 解算步骤
综上所述，空间后方交会的求解过程如下：
- 获取已知数据：从摄影资料中查取像片比例尺\(1/m\)，平均航高，内方元素\(x_0\)、\(y_0\)、\(f\)；从外业测量成果中，获取控制点的地面测量坐标\(X_t\)、\(Y_t\)、\(Z_t\)，并转化成地面摄影测量坐标\(X\)、\(Y\)、\(Z\)。
- 量测控制点的像点坐标：将控制点标刺在像片上，利用立体坐标量测仪量测控制点的像框标坐标，并经像点坐标改正，得到像点坐标\(x\)、\(y\)。
- 确定未知数的初始值：在竖直摄影情况下，角元素的初始值为\(0\)，即\(\varphi_0 = \omega_0 = \kappa_0 = 0\)；线元素中，\(Z_{s0} = H = mf\)，\(X_{s0}\)、\(Y_{s0}\)的取值可用四个角上控制点坐标的平均值，即：\(X_{s0} = \frac{1}{4} \sum_{i = 1}^{4} X_i\)，\(Y_{s0} = \frac{1}{4} \sum_{i = 1}^{4} Y_i\)。
- 计算旋转矩阵\(R\)：利用角元素的近似值计算方向余弦值，组成\(R\)阵。
- 逐点计算像点坐标的近似值：利用未知数的近似值按共线方程式计算控制点像点坐标的近似值\(\tilde{x}\)、\(\tilde{y}\)。
- 组成误差方程式：按式 \(\eqref{对XYZ偏导}\eqref{对角度偏导化简}\eqref{误差方程}\) 逐点计算误差方程式的系数和常数项。
- 组成法方程式：计算法方程的系数矩阵\(A^{\mathrm{T}}A\)与常数项\(A^{\mathrm{T}}L\)。
- 解求外方位元素：根据法方程，按式\(\eqref{封闭解}\)解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。
- 检查计算是否收敛：将求得的外方位元素的改正数与规定的限差比较，小于限差则计算终止，否则用新的近似值重复第4至第8步骤的计算，直到满足要求为止。
4. 精度
由平差理论可知，法方程系数的逆矩阵 \((A^{\mathrm{T}}A)^{-1}\) 等于未知数的协因数阵 \(Q_{x}\)，因此可按下式计算未知数的中误差： \[ m_{i}=m_{0} \cdot \sqrt{Q_{ii}} \] 式中，\(i\) 表示相应的未知数，\(Q_{ii}\) 为 \(Q_{x}\) 阵中的主对角线元素，\(m_{0}\) 称为单位权中误差，计算公式为： \[ m_{0}=\pm\sqrt{\frac{\|V^TV\|_{2}^{2}}{2n - 6}} \] 这里，\(n\) 表示控制点的总数。
5. python 实现
以 [1] 的课后题为例
目录结构：
1 | . |
1 | import numpy as np |
1 | NO. x(mm) y(mm) X(m) Y(m) Z(m) |
1 | 外方位元素及其精度： |
如果你觉得这篇博客对你有所帮助，不妨在 GitHub 上给它点个
star🌟，这将是对我莫大的鼓励。
Github
地址
6. 参考资料
[1] 王佩军，徐亚明等编著.摄影测量学[M].武汉：武汉大学出版社，2016.5