单像空间后方交会是摄影测量领域的基础性课题^[1-2]，它是指利用一定数量的地面控制点及像点求解像片外方位元素的过程。摄影测量中常借助快速、正确地恢复图像获取时的空间方位来确定地面目标的空间位置^[3]。因此，对单像空间后方交会的解算进行研究具有典型的实际意义。空间后方交会是一个非线性模型，一般采用迭代方法进行解算^[4]，典型的解法包括欧拉角法^[1-2]、角锥体法^[5]及直接线性变换法^[6]，其中角锥体法对于外方位线元素初值依赖较强，直接线性变换法不适合高精度的测量任务^[7]。欧拉角法将成像姿态用欧拉角描述，对共线方程线性化，并利用最小二乘法迭代求解^[2]。在应用中，欧拉角法的限制主要表现为：①需要已知良好的外方位元素初值, 当初值与真值差距较大时，收敛速度慢，甚至不收敛；②对原始非线性模型用线性模型近似代替，有一定的模型误差，随着观测精度的提高，因线性近似所产生的模型误差与观测误差相当，有的甚至超过观测误差^[8]；③其迭代过程中需假定法方程系数矩阵满秩，否则将可能失效^[9]。虽然外方位元素初值可从全球定位系统(Global Positioning System，GPS)、惯性导航系统和轨迹数据中获取，但无人机等平台自带GPS及姿态量角系统本身的精确性不高，并存在记录滞后的现象，通常不能直接使用^[10]。

为减少初值依赖及线性近似的影响，学者们提出了许多有参考价值的方法，主要分为直接解法和迭代解法两类^[11]。直接解法无需外方位元素初值，但解算精度与鲁棒性较低，文献[12]利用单位四元数分别解算外方位线元素与角元素。代表性的迭代解法如文献[13-14]以线性解为初值进行高斯-牛顿全局优化；文献[15-16]采用对偶四元数描述坐标系之间的平移和旋转。迭代解法虽具有较高精度和鲁棒性，但计算复杂度高且对参数初值依然具有依赖性。文献[17]指出，在控制点间高差较小时，虽空间后方交会模型为良态，但其法方程系数矩阵也会存在严重的病态问题；文献[18]提出当像平面与控制点平面平行时，法方程系数矩阵的病态性将大大加强。对于病态方程，通常采用正则化方法求解^[19-22]。文献[23]采用谱修正迭代算法来改善方程的病态性，用以保持法方程在解算过程中的数值稳定性；文献[9]则利用Levenberg-Marquardt(LM)方法来避免法方程系数矩阵病态时造成的不收敛问题。

为减少线性化影响及加快收敛速度，考虑用牛顿法代替欧拉角法中的最速下降法，即对共线方程二阶泰勒近似，用二次函数模型近似原始非线性模型。非线性迭代求解的本质是在给定初值的前提下，计算梯度矩阵或海森矩阵决定下降方向以达到函数极值点，此时，梯度矩阵和海森矩阵的病态性强弱将决定法方程计算时是否考虑系数矩阵病态性的影响。为此，作为一种方法上的有效补充，本文提出一种求解单像空间后方交会的监督学习方法，在摄站或摄影基线及其附近一定范围内模拟生成样本集，并用样本集的整体下降方向来代替牛顿法中的牛顿方向，从而避免梯度矩阵和海森矩阵的直接求解，以消除牛顿法中因梯度矩阵和海森矩阵病态性对法方程求解的影响。此外，将函数凹凸性隐藏在学习过程中，同时达到增强初值选取鲁棒性的目的。

1 监督学习方法求解单像空间后方交会 1.1 非线性过程的二次表达

设$\hat \theta $为待估计参数，V为改正数，Y为观测值，则原始误差方程为：

$ \mathit{\boldsymbol{V}} = f\left( {\hat \theta } \right) - \mathit{\boldsymbol{Y}} $

(1)

式中，f为非线性函数，非线性最小二乘法就是在式(1)的条件下，求解以下目标函数^[24]：

$ \begin{array}{*{20}{c}} {\hat \theta = \arg \min {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{PV}} = }\\ {\arg \min \left( {{{\left( {f\left( {\hat \theta } \right) - Y} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {f\left( {\hat \theta } \right) - Y} \right)} \right)} \end{array} $

(2)

式中，P为权矩阵，当同精度观测时为单位矩阵，即P=I。令函数F(θ)=V^TV，则式(2)可化简为：

$ \begin{array}{*{20}{c}} {\hat \theta = \arg \min {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}} = \arg \min F\left( {\hat \theta } \right) = }\\ {\arg \min \left( {{{\left( {f\left( {\hat \theta } \right) - Y} \right)}^{\rm{T}}}\left( {f\left( {\hat \theta } \right) - Y} \right)} \right)} \end{array} $

(3)

牛顿法是解决式(3)中非线性最小二乘的常用工具，诸多其他方法均以其为基础进行改进。它假设在最小值的邻域内，函数F(θ)可近似为一个二次函数，即对式(3)按二阶泰勒展开：

$ \begin{array}{*{20}{c}} {F\left( \theta \right) \approx F\left( {\hat \theta } \right) + \nabla _F^{\rm{T}}\left( {\hat \theta } \right)\left( {\theta - \hat \theta } \right) + }\\ {\frac{1}{2}{{\left( {\theta - \hat \theta } \right)}^{\rm{T}}}{H_F}\left( {\hat \theta } \right)\left( {\theta - \hat \theta } \right)} \end{array} $

(4)

式中，${\nabla _F}\left( {\hat \theta } \right)$与${\mathit{\boldsymbol{H}}_F}\left( {\hat \theta } \right)$是函数F在估计值${\hat \theta }$处的梯度矩阵和海森矩阵。式(4)对(θ-${\hat \theta }$)求偏导并令其为0，若海森矩阵正定，则可得到非线性最小二乘牛顿法迭代解算更新公式：

$ {\hat \theta _{k + 1}} - {\hat \theta _k} = - H_F^{ - 1}\left( {{{\hat \theta }_k}} \right){\nabla _F}\left( {{{\hat \theta }_k}} \right) $

(5)

牛顿法只是局部收敛^[25]，并且在每次迭代中要计算函数在当前估计值处的梯度矩阵和海森矩阵。

结合式(3)、式(5)，将式(5)中的${\nabla _F}\left( {{{\hat \theta }_k}} \right)$按链式法则展开：

$ {\hat \theta _{k + 1}} - {\hat \theta _k} = - 2H_F^{ - 1}\left( {{{\hat \theta }_k}} \right){\nabla _f}\left( {{{\hat \theta }_k}} \right)\left( {f\left( {{{\hat \theta }_k}} \right) - Y} \right) $

(6)

为提高计算精度，借鉴最小二乘法多次迭代收敛的思想，可得到监督学习方法解决非线性最小二乘问题第k次迭代的更新公式：

$ {\hat \theta _{k + 1}} - {\hat \theta _k} = {R_k}\left( {f\left( {{{\hat \theta }_k}} \right) - Y} \right) $

(7)

式中，${\mathit{\boldsymbol{R}}_k} = - 2\mathit{\boldsymbol{H}}_F^{ - 1}({\hat \theta _k}){\nabla _f}({\hat \theta _k})$, 称为第k次迭代的整体下降方向，其通过最小化m个样本构成的训练集中所有样本第k+1次迭代的预测值和真值之间的差得到，即：

$ \begin{array}{l} \min \sum\limits_{i = 1}^m {{{\left\| {\theta _ * ^i - \hat \theta _{k + 1}^i} \right\|}^2}} = \\ \min \sum\limits_{i = 1}^m {{{\left\| {\theta _ * ^i - \left( {\hat \theta _k^i + {R_k}\left( {f\left( {\hat \theta _k^i} \right) - Y_ * ^i} \right)} \right.} \right\|}^2}} \end{array} $

(8)

式中，$\theta _*^i$为第i个样本的待估计参数真值(训练阶段已知)；$\hat \theta _k^i$为第i个样本第k次迭代待估计参数的估计值；$\mathit{\boldsymbol{Y}}_*^i = f\left( {\theta _*^i} \right)$。

由式(7)、式(8)可知，监督学习方法的解算过程如图 1所示，具体步骤为：①训练阶段，给定由多个样本组成的训练集，其中每个样本均已知待估计参数的真值与初值，监督学习方法从训练数据中学习一个使得所有样本都由初值下降到真值的整体下降方向R_k，即最小化式(8)。为提高预测精度，采用多次迭代，得到R_k，将其代入式(7)更新后得到的${\hat \theta _{k + 1}}$作为新的初值，再利用新初值学习一个新的整体下降方向R_k+1，经k次迭代后，得到整体下降方向组成的集合{R_k}。②测试阶段，给定待估计参数初值，依次按训练阶段得到的整体下降方向集合{R_k}，迭代式(7)进行待估计参数的优化。

监督学习方法与牛顿法解决非线性最小二乘的区别如图 2所示。图 2中，牛顿法每次的迭代方向由函数在当前参数估计值处的梯度矩阵和海森矩阵共同决定。监督学习方法是用训练阶段学习得到的整体下降方向进行优化，即是沿着初值到真值的理想方向进行，如图 2(b)所示，第0个样本第k+1次迭代中，初值到真值的理想方向和优化更新方向分别用虚线箭头和实线箭头表示。由图 2可知，实线箭头沿着虚线箭头方向进行优化。

由上述可知，监督学习方法较之牛顿法有以下优点：

1) 无需计算函数的梯度矩阵和海森矩阵，因此，对函数模型是否可导无要求。

2) 式(5)成立的前提是函数模型在当前参数估计值处的海森矩阵正定，但对于其海森矩阵非正定的情况，牛顿法无法直接通过式(5)迭代求解。而由式(7)、式(8)可知，监督学习方法在求解整体下降方向R_k以及预测新估计值的过程中只用到了当前估计值的函数值，所以求解过程不会受到梯度矩阵和海森矩阵病态性的影响。因此，对于牛顿法中因法方程系数矩阵病态而无法解决的问题，监督学习方法可用从训练数据中学习到的整体下降方向进行待估计参数的优化。

3) 无初值依赖性，由式(7)可知待估计参数当前值${\hat \theta _k}$与计算值${\hat \theta _{k + 1}}$由整体下降方向R_k决定。而式(8)中R_k的计算不依赖于函数的梯度矩阵与海森矩阵，即与函数的凹凸性质无关。因此，不会出现牛顿法中因初值选取不当而造成无法收敛的情况。

1.2 单像空间后方交会监督学习过程

当使用欧拉角描述摄站姿态时，设像片内方位元素已知，将待求解外方位元素采用列向量形式表示为$\mathit{\boldsymbol{p = }}{\left( {{X_S}{Y_S}{Z_S}\varphi \omega \kappa } \right)^{\rm{T}}}$，则共线方程表示为^[26]：

$ \left\{ \begin{array}{l} {h_x}\left( p \right) = x - {x_0} = - 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)}}\\ {h_y}\left( p \right) = y - {y_0} = - 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} \right. $

(9)

式中，$\left( {{x_0},{y_0},f} \right)$为像片内方位元素；(x, y)为像点坐标；(X, Y, Z)为地面坐标；(X_S, Y_S, Z_S)为摄站在地面坐标系中的坐标；${a_i},{b_i},{c_i}\left( {i = 1,2,3} \right)$为变换前后两坐标轴相应夹角的余弦。

令$\mathit{\boldsymbol{h}} = {({h_x}{h_y})^{\rm{T}}}$，$\mathit{\boldsymbol{l}} = {(x - {x_0}y - {y_0})^{\rm{T}}}$，则可得每个点的误差方程：

$ \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{h}}\left( {\mathit{\boldsymbol{\hat p}}} \right) - \mathit{\boldsymbol{l}} $

(10)

为提高外方位元素的求解精度，常有多余观测。若有n个控制点，则可依据式(10)建立n组误差方程$\mathit{\boldsymbol{V = }}{\left( {{v_1}{v_2} \cdots {v_n}} \right)^{\rm{T}}}$，因此，单像空间后方交会问题可描述为如下最小二乘过程，令$F\left( {\mathit{\boldsymbol{\hat p}}} \right) = {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}}$，则：

$ \begin{array}{l} \mathit{\boldsymbol{\hat p}} = \arg \min F\left( {\hat p} \right) = \\ \;\;\;\;\;\;\arg \min \left( {{{\left( {H\left( {\hat p} \right) - L} \right)}^{\rm{T}}}\left( {H\left( {\hat p} \right) - L} \right)} \right) \end{array} $

(11)

式中，$\mathit{\boldsymbol{H = }}{\left( {{\mathit{\boldsymbol{h}}_1}{\mathit{\boldsymbol{h}}_2} \cdots {\mathit{\boldsymbol{h}}_n}} \right)^{\rm{T}}}$; $\mathit{\boldsymbol{L = }}{\left( {{\mathit{\boldsymbol{l}}_1}{\mathit{\boldsymbol{l}}_2} \cdots {\mathit{\boldsymbol{l}}_n}} \right)^{\rm{T}}}$。依据1.1节非线性最小二乘的描述可知，若函数在当前外方位元素估计值${\mathit{\boldsymbol{\hat p}}_k}$处的海森矩阵${\mathit{\boldsymbol{H}}_F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)$正定，则可得单像空间后方交会的牛顿法迭代更新公式：

$ {\hat p_{k + 1}} - {\hat p_k} = - H_F^{ - 1}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right){\nabla _F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) $

(12)

式中，${\nabla _F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)$和${\mathit{\boldsymbol{H}}_F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) = \nabla \left( {{\nabla _F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)} \right)$是函数F在当前外方位元素估计值处的梯度矩阵和海森矩阵。将${{\nabla _F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)}$用链式法则展开：

$ \begin{array}{*{20}{c}} {{\nabla _F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) = 2{\nabla _V}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)\mathit{\boldsymbol{V}} = }\\ {2{\nabla _H}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)\left( {H\left( {\mathit{\boldsymbol{\hat p}}} \right) - L} \right)} \end{array} $

(13)

式中，

$ \begin{array}{l} {\nabla _H}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) = \left( {{\nabla _{{h_1}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right){\nabla _{{h_2}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) \cdots {\nabla _{{h_n}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)} \right) = \\ \;\;\;\;\;\;\left( {{\nabla _{{h_{1x}}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right){\nabla _{{h_{1y}}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) \cdots {\nabla _{{h_{nx}}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right){\nabla _{{h_{ny}}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)} \right) \end{array} $

(14)

其中,

$ {\nabla _{{h_{1x}}}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) = {\left( {\frac{{\partial {h_{1x}}\left( {\mathit{\boldsymbol{\hat p}}} \right)}}{{\partial {X_S}}}\frac{{\partial {h_{1x}}\left( {\mathit{\boldsymbol{\hat p}}} \right)}}{{\partial {Y_S}}}\frac{{\partial {h_{1x}}\left( {\mathit{\boldsymbol{\hat p}}} \right)}}{{\partial {Z_S}}}\frac{{\partial {h_{1x}}\left( {\mathit{\boldsymbol{\hat p}}} \right)}}{{\partial \varphi }}\frac{{\partial {h_{1x}}\left( {\mathit{\boldsymbol{\hat p}}} \right)}}{{\partial \omega }}\frac{{\partial {h_{1x}}\left( {\mathit{\boldsymbol{\hat p}}} \right)}}{{\partial \kappa }}} \right)^{\rm{T}}} $

(15)

将式(13)${{\nabla _F}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)}$的值代入式(12)，并令${\mathit{\boldsymbol{R}}_k} = - 2\mathit{\boldsymbol{H}}_F^{ - 1}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right){\nabla _H}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right)$，可得监督学习方法第k次迭代下单像空间后方交会解算的更新公式：

$ {\mathit{\boldsymbol{\hat p}}_{k + 1}} - {\mathit{\boldsymbol{\hat p}}_k} = {\mathit{\boldsymbol{R}}_k}\left( {\mathit{\boldsymbol{H}}\left( {{{\mathit{\boldsymbol{\hat p}}}_k}} \right) - L} \right) $

(16)

由1.1节监督学习方法解算过程可知，式(16)中的R_k可通过最小化m个样本所构成的训练集中，所有样本第k+1次迭代的预测值与真值之间的差得到，即：

$ \begin{array}{l} \min \sum\limits_{i = 1}^m {{{\left\| {\mathit{\boldsymbol{p}}_ * ^i - \mathit{\boldsymbol{\hat p}}_{k + 1}^i} \right\|}^2}} = \\ \;\;\;\;\;\min \sum\limits_{i = 1}^m {{{\left\| {\mathit{\boldsymbol{p}}_ * ^i - \left( {\mathit{\boldsymbol{\hat p}}_k^i + {\mathit{\boldsymbol{R}}_k}\left( {H\left( {\hat p_k^i} \right) - L_ * ^i} \right)} \right)} \right\|}^2}} \end{array} $

(17)

式中，p_*ⁱ表示第i个样本，即第i个外方位元素的真值；$\mathit{\boldsymbol{\hat p}}_k^i$为第i个样本第k次迭代外方位元素的估计值；$\mathit{\boldsymbol{L}}_*^i = H\left( {\mathit{\boldsymbol{p}}_*^i} \right)$表示控制点在外方位元素真值p_*ⁱ下的像点坐标。

由式(16)、式(17)可知，采用监督学习方法求解单像空间后方交会的两个条件为：

1) 训练阶段需要同一测区内多组已知外方位元素的像片组成训练集对外方位元素的整体下降方向进行学习；

2) 因不同测区地面控制点分布一般不具有相似性，为此，训练阶段与测试阶段为同一地区同一组控制点。

在给定一测区内一组地面控制点及像片内方位元素的前提下，考虑采用由该测区多组模拟的外方位元素作为样本组成训练集进行外方位元素整体下降方向的学习。其具体计算过程如下：

1) 样本数据准备阶段

① 原始数据，包括像片内方位元素、地面控制点在地面坐标系中的坐标M∈R^3×n，n表示地面控制点个数。

② 随机产生m个外方位元素${\mathit{\boldsymbol{P}}_*} = \left( {\mathit{\boldsymbol{p}}_*^1\mathit{\boldsymbol{p}}_*^2 \cdots \mathit{\boldsymbol{p}}_*^m} \right)$，P_*∈R^6×m。其中p_*^m∈R⁶为${\left( {{X_S}{Y_S}{Z_S}\varphi \omega \kappa } \right)^{\rm{T}}}$，表示第m个外方位元素，并通过共线方程分别计算M在m个外方位元素下的像点坐标${\mathit{\boldsymbol{U}}_*} = \left( {\mathit{\boldsymbol{L}}_*^1\mathit{\boldsymbol{L}}_*^2 \cdots \mathit{\boldsymbol{L}}_*^m} \right)$，U_*∈R^2n×m，其中$\mathit{\boldsymbol{L}}_*^m = \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{p}}_*^m) \in {\bf{R}^{2\mathit{\boldsymbol{n}}}}$为${({x_1}{y_1} \cdots {x_n}{y_n})^{\rm{T}}}$，表示外方位元素p_*^m下的像点坐标。

2) 训练阶段

① 初始化外方位元素${\mathit{\boldsymbol{P}}_k} = (\mathit{\boldsymbol{p}}_k^1\mathit{\boldsymbol{p}}_k^2 \cdots \mathit{\boldsymbol{p}}_k^m) = 0$，此时迭代次数k=0，P_k∈R^6×m，其中p_k^m∈R⁶为${\left( {{X_S}{Y_S}{Z_S}\varphi \omega \kappa } \right)^{\rm{T}}}$，表示第m个外方位元素第k次迭代的计算值。

② 通过式(9)计算地面控制点在当前外方位元素P_k下的像点坐标${\mathit{\boldsymbol{U}}_k} = (\mathit{\boldsymbol{L}}_k^1\mathit{\boldsymbol{L}}_k^2 \cdots \mathit{\boldsymbol{L}}_k^m)$，U_k∈R^2n×m，其中$\mathit{\boldsymbol{L}}_k^m = \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{p}}_k^m) \in {\bf{R}^{2n}}$为${({x_1}{y_1} \cdots {x_n}{y_n})^{\rm{T}}}$，表示外方位元素p_k^m下的像点坐标，计算U_*-U_k及P_*-P_k的值。通过解线性方程组${\mathit{\boldsymbol{P}}_*} - {\mathit{\boldsymbol{P}}_k} = {\mathit{\boldsymbol{R}}_k}\left( {{\mathit{\boldsymbol{U}}_*} - {\mathit{\boldsymbol{U}}_k}} \right)$得到整体下降方向R_k，其中R_k∈R^6×2n。

③ 用学习得到的R_k，代入式(16)计算P_k+1，k=k+1返回步骤②，直到满足指定的迭代次数或达到收敛条件为止。

④ 输出整体下降方向集合{R_k}。

3) 测试阶段

① 待估计外方位元素像片的像点观测值${\mathit{\boldsymbol{L}}_*} = {({x_1}{y_1} \cdots {x_n}{y_n})^{\rm{T}}}$。

② 初始化外方位元素p_k=0∈R⁶，此时迭代次数k=0。

③ 由式(9)计算地面控制点在当前外方位元素p_k下的像点${\mathit{\boldsymbol{L}}_k} = {(x_1^ky_1^k \cdots x_n^ky_n^k)^{\rm{T}}}$及L_k-L_*。代入式(16)，用{R_k}中R_k计算p_k+1，并将p_k+1作为外方位元素新的估计值。

④ 迭代次数k=k+1，返回步骤③，直至{R_k}遍历完为止。

⑤ 输出外方位元素估计值。

2 监督学习方法单像空间后方交会有效性试验

为验证本文监督学习方法解决单像空间后方交会的有效性，分别进行病态性试验与初值依赖性试验。

2.1 病态性试验

假设存在理想情况，如图 3所示。已知摄像机焦距f为35 mm，像平面坐标系的o'-xy平面与世界坐标系O-XYZ的O-XY平面平行。航高为500 m，摄影中心S在世界坐标系的坐标为(0，0，150)，单位为m，像片外方位角元素的真实值为(0，0，0)，单位为rad。假设地面控制点的坐标分别为A(-50，50，0)，B(-50，-50，0)，C(50，-50，0)，D(50，50，0)，单位为m，则根据三角形相似性，可得对应像点坐标a(-3.5，3.5)，b(-3.5，-3.5)，c(3.5，-3.5)，d(3.5，3.5)，单位为mm。根据地面控制点及其对应像点坐标，采用欧拉角法计算像片的外方位元素。由于像平面与控制点平面平行，导致欧拉角法中的法方程系数矩阵病态，因此，无法直接通过法方程解算此时的外方位元素^[17]。

以线元素均值为0、方差为0.1，以角元素均值为0、方差为0.05的高斯噪声，在外方位元素真实值(0，0，150，0，0，0)的基础上扰动产生多个外方位元素样本组成监督学习方法的训练集，训练整体下降方向并求解以上问题中的外方位元素，计算结果如表 1所示。

由表 1可知，当欧拉角法中因法方程系数矩阵病态而无法求解时，监督学习方法仍能准确解算，且能保持较高的计算精度。因为其通过训练集来学习获得整体下降方向，不受到梯度矩阵或海森矩阵病态的限制。

2.2 初值依赖性试验

给定地面控制点坐标，如表 2所示。

试验方案如图 4所示。已知像片参数为：镜头焦距f为24 mm，相对航高H为114 m，x₀=y₀=0，以表 2中第5个控制点(即图 4中A)为被摄地面的中心点，坐标为A(X_A, Y_A, 0)。基于计算方便及像片上控制点全部可见两点考虑，设摄影基线的方程为：Y=Y_A, Z=H。

为保证定位精度最佳，选取φ∈[-45°, 45°]。则图 4中B、C、D的外方位元素分别为B(X_A Y_A H 0 0 0)，C(X_A - H Y_A H φ_max 0 0)，D(X_A + H Y_A H φ_min 0 0)。以B为基准，沿摄影基线以Δx=±13 m为步长、C、D为终点产生如图 4所示的系列基准摄站，每个摄像机光轴对准地面中心点A。即B左边第i个摄站外方位元素为(X_A-i×Δx Y_A H、$\arctan \frac{{i \times \Delta x}}{H}$ 0 0)，右边第i个摄站外方位元素为(X_A+i×Δx Y_A H-$\arctan \frac{{i \times \Delta x}}{H}$ 0 0)。再以每个基准摄站为基准，线元素增量范围为$\left[ { - \frac{{\Delta x}}{2},\frac{{\Delta x}}{2}} \right]$，角元素增量范围为[-5°, 5°]，随机产生新的外方位元素。试验所有基准摄站共产生了1 000个外方位元素作为训练集进行训练。

测试阶段，为使训练得到的整体下降方向集合能准确解算测试阶段的外方位元素，假设测试阶段的外方位线元素与角元素的最大范围不超过训练阶段的训练样本。具体为：以Δx=±8 m产生系列基准摄站，在每个基准摄站的基础上，以线元素增量范围为$\left[ { - \frac{{\Delta x}}{2},\frac{{\Delta x}}{2}} \right]$，角元素增量范围为[-3°, 3°], 随机产生新的外方位元素作为测试样本。试验共产生301个外方位元素。

计算地面控制点在每个测试样本下的像点坐标，采用训练阶段得到的整体下降方向计算得到外方位元素的计算值，并与其随机产生的真实值对比，验证本文方法的正确性。考虑到计算效率，采用Eigen进行编程实现，训练过程花费时间为2 521 ms，平均计算每个测试样本的时间花费约为5 ms。为方便比较本文方法的准确性，将计算结果投影至二范数空间进行比较。为验证本文方法的初值依赖性，将训练阶段与测试阶段外方位元素的初值设为0，即$X_S^0 = Y_S^0 = Z_S^0 = 0$，${\varphi ^0} = {\omega ^0} = {\kappa ^0} = 0$。剔除计算错误的点后，所训练的整体下降方向准确解算测试样本数达283个，占全部测试样本的94%，计算结果如图 5所示。

由图 5可知，计算值曲线与真值曲线基本重合。经统计，测试集中摄站位置最大差值为0.016 m，摄影姿态最大差值为0.001°，满足测量精度要求^[15]，表明该方法的正确性，可用于求解单像空间后方交会。

为验证本文方法的初值依赖性，分别采用本文方法、欧拉角法和谱修正方法解算上述的283个测试样本，并将计算结果投影至二范数空间进行对比。由于试验采用交向摄影，部分倾角较大，无法直接采用欧拉角法进行计算，因此，采用直接线性变换方法提供初值，3种方法的计算结果误差曲线如图 6所示，迭代次数如图 7所示。

由图 6、图 7可知，在计算精度上，欧拉角法、监督学习方法与谱修正方法对于外方位元素的求解精度相当。但欧拉角法依赖于迭代初值，会因迭代初值的选取不当而无法收敛。因此，在不提供初值的前提下，较难解算大姿态下像片的外方位元素。在迭代次数上，监督学习方法总体优于其他两种算法，具有较快的收敛速度。虽然监督学习方法的迭代更新公式是在牛顿法的基础上推导的，但较之牛顿法却无较强的初值依赖且能保持较快的收敛速度。其原因在于，监督学习方法是通过大量的样本，在已知外方位元素初值和真值的前提下，通过每次迭代中最小化所有样本当前外方位元素初值到真值之差，最终得到整体下降方向集合{R_k}，并以此解算测试阶段待求像片的外方位元素。因此，只需测试样本与训练样本近似，且外方位元素的初值相同，则利用训练阶段得到的整体下降方向集合就可以准确解算待求的外方位元素。

从测试数据解算结果可知，监督学习方法具有解算精度较高、收敛速度较快及初值依赖性弱3个特点，但在测试阶段要求待求像片的所有像点均可见，且要求在同一航线上对测区进行不同倾角的重复观测，需采用多重交向摄影的方式进行试验。

3 结语

本文提出了一种求解单像空间后方交会的监督学习方法，并在牛顿法的基础上推导了迭代更新公式。对比试验验证了监督学习方法克服了法方程系数矩阵病态、快速收敛及初值依赖性弱的缺点，与此同时，能够保持较高的解算精度与效率。由于监督学习方法是对摄影基线上系列基准摄站扰动产生训练样本以获得整体下降方向，因此，目前只能将其应用于摄影基线附近一定范围内的像片外方位元素的解算。文中试验采用的摄影基线比较简单，较为复杂的摄影基线下外方位元素的解算有待进一步研究。此外，采用交向摄影构成的测量网与常规航空摄影测量不同，并且线阵CCD（charge coupled devices）比面阵CCD的姿态运动方程复杂，能否将本文方法分别应用于光束法平差和线阵CCD像片的外方位元素求解是今后的研究方向。