摄影测量有着悠久的历史，目前已广泛用于不同比例尺的地形测量及工业、建筑等非地形测量^[1-3]。尽管其对象与任务千差万别，但主要方法都是从二维影像重建三维模型^[4-5]，在重建的三维模型上提取所需的各种信息，其中重建三维模型的基础就是根据影像确定目标点的空间位置。作为一种常用定位方法，立体像对的空间前方交会是根据立体像对左右两影像的内、外方位元素和同名像点的影像坐标量测值来确定该点的物方空间坐标^[6]。这种方法易于理解，运算简便，却无法充分利用更多像片的观测数据。一旦立体像对中的某张影像量测数据包含较大误差，则必然导致错误的解算结果，即其抗差性能较差。随着近景、航空、航天等摄影测量技术的发展^[7-10]，包含同一目标点的影像在数量上已经相当可观，充分利用这些像片信息，并尽可能降低误差对解算结果的影响，最终实现地面点定位，具有重要意义。文献[11]曾提出基于多张像片前方交会的人工标志点匹配方法，并取得较好的匹配效果，然而该方法仍是通过依次实现双片空间前方交会的方法进行物方点筛选。在地面点定位过程中，若能直接利用所有的观测信息，不仅可以减少依次筛选与逐个剔除而引起的时间耗费，也能有效避免因个别像片观测条件差而引起的定位偏差，最终解算出更精确的地面点位置。总体最小二乘方法^[12-18]顾及了所有的观测信息，被广泛用于解决测量学中观测数据含有误差的问题，均取得较好成效。考虑到空间前方交会的主要思路是“同名射线必在物方空间交于一点”，而实际作业中通常会有测量误差的存在，同名射线不会严格相交，因此，可利用总体最小二乘方法，根据“目标点到多条同名射线的距离平方和最短”这一原则构建目标函数来实现地面点定位^[19-21]。最后，通过引入稳健估计理论和借助选权迭代法以规避粗差干扰，并为精度高的观测信息赋予较高的权重^[22-27]，给出可直接用于多张像片的空间前方交会公式。

1 多像片前方交会定义

如图 1所示，在多个摄影站点S₁、S₂、…、S_n对地面进行摄影，可获取n张像片。地面点A在这n张像片的构象分别为a₁、a₂、…、a_n。已知n张像片的内、外方位元素，将像片按内、外方位元素置于摄影时的位置，在理想情况下，同名射线S₁a₁、S₂a₂、…、S_na_n必将交于一点，即地面点A。这种利用多张像片的内、外方位元素及同名像点的像空间坐标，通过交会解算相应模型点坐标的过程，称为多张像片的空间前方交会。

2 多像片前方交会过程 2.1 构建目标函数

考虑到像点的量测误差等因素，这些同名射线并不一定交于同一点甚至不一定相交，在综合考虑所有观测信息的条件下，可通过构建目标函数，求到同名射线距离的平方和最小的点作为A点

式中，p_i为权重；是点A到射线S_ia_i的距离。记点A、各像片摄影中心S_i在地面摄影测量坐标系中的坐标分别为X∈R^3×1、X_Si∈R^3×1，像点a_i的像空间坐标为x_i∈R^3×1，各像片的旋转矩阵为R_i∈R^3×3，则根据坐标系间的旋转关系，可得射线S_ia_i的方向向量在地面摄影测量坐标系中可表示为R_ix_i。因此，式(1) 可等价为

考虑到矢量叉乘与点乘之间的关系，则式(2) 可等价变换为

式中，各像片权重p_i，旋转矩阵R_i，摄影中心X_Si，像点坐标x_i均为已知，仅地面点A的坐标X为未知。通过解算“到同名射线距离的平方和最小”这一条件而得出的X，即为考虑了所有观测点信息的地面点坐标。

2.2 总体最小二乘估计

为使得地面点A的坐标X满足目标函数φ=min这一条件，即，因此可对式(3) 逐项展开并对X求一阶导数，如下

由于式(4) 中X^TR_ix_i是一个标量，满足其转置等于自身，有X^TR_ix_i=(X^TR_ix_i)^T=x_i^TR_i^TX，将式(4) 进一步化简，可得

式中, I为单位阵，当n=1时，I-R_ix_ix_i^TR_i^T/|x_i|²的行列式为0，即系数矩阵不可逆，无法进一步求出地面点A的坐标，这是因为单个像点及其投影中心只能确定地面点的空间方向，而无法确定其位置。当n≥2，若n条射线S_ia_i相互平行时，此时每张像片上I-R_ix_ix_i^TR_i^T/|x_i|²项相同，系数矩阵可合并为(I-R₁x₁x₁^TR₁^T/|x₁|²)，其行列式为0，不可逆，因此也无法求出地面点A的坐标，这是因为相互平行的n条射线在空间中并不相交，满足到它们距离平方和最小的点位于一条直线上，而无法确定其具体位置。

由于式(5) 系数矩阵的每一项中I-R_ix_ix_i^TR_i^T/|x_i|²具有特征值0、1、1，故对应二次型χ^T(I-R_ix_ix_i^TR_i^T/|x_i|²)χ(χ≠0) 非负。在其他情况下，即存在n≥2条不平行的同名射线时，χ^T(I-R_ix_ix_i^TR_i^T/|x_i|²)χ不可能同时为0。因此，可得式(5) 系数矩阵的二次型

故式(5) 系数矩阵正定，可逆。因此，必然可确定满足目标函数最小的地面点A，为

由于目标函数对X的一阶导数为0的点有且仅有一个，且其二阶导数即为式(5) 的系数矩阵，根据上述判定，可知其正定，说明该驻点确实是极小值点，故式(7) 即为地面点A的加权总体最小二乘估计。

2.3 选权迭代法

由于在实际测量工作中客观条件的限制，很难完全避免粗差的存在或做到完全同等精度量测。在平差过程中，通常引入权作为比较观测值之间相对精度高低的指标，并为精度较高的观测数据赋予较高的权重。但是，在各像片量测误差分布未知的情况下，难以对各像片进行定权。随着解算地面点坐标时可利用像片数目的增多，为使参数的估值尽可能充分利用观测数据中的有效信息，规避有害信息的干扰，经常会引入稳健估计的方法。本文选用IGG选权迭代法来为观测像片赋权。主要步骤如下：

(1) 令目标函数中p_i初值均为1或其他先验权，利用式(7)，得地面点A坐标的第一次估值。

(2) 分别求地面点A到各射线S_ia_i的距离d_i，及标准差，如下

(3) 根据IGG方案^[27]，定义等价权因子w_i，如下

式中，(淘汰点)。

(4) 令代替式(7) 中的p_i，重新估计新的地面点A坐标。

(5) 重复(2)~(4) 步，进行迭代运算，直至两次估值结果差异小于限差为止。

3 算例分析

为验证本文方法的可靠性，并分析误差对于该方法的影响，以及其对粗差的探测效果，假定地面点A的地面摄影测量坐标真值X₀=[200 100 50]^T，单位为m以及3组已知内、外方位元素的像片(见表 1—表 3)，其中f为焦距；X_S=(X_S Y_S Z_S)^T为摄影中心的地面摄影测量坐标；φ、ω、κ为像片旋转角。

3.1 算法可靠性验证

3.1.1 算法正确性验证

为验证算法的正确性，利用严格共线方程可以推算出该地面点在不同像片上对应像点的像空间坐标x_i。之后，将这3组像片分别进行空间前方交会，可得出地面点A的新坐标X，最后将其与真值进行比较。

该算例中已知每组6张像片的内、外方位元素，地面点坐标真值按照严格共线方程模拟出每张像片上对应的像点坐标，此时可认为各像片具有同等精度条件，即各像片的权重均为1，利用本文方法解算出的地面点坐标为X=[200 100 50]^T，单位为m与真值相同，该算例验证了本文算法的正确性。

3.1.2 算法稳健性验证

为验证本文方法的稳健性，以第1组数据为例，在各像片像点坐标上加入0.5个像素(像素大小为4 μm)的高斯噪声，并在像片1中的像点横坐标加入2 mm的粗差，再利用本文方法进行多张像片的空间前方交会，ds为交会结果与真值的位移偏差。初次运算时为每张像片赋予权重为1，每次运算的权和地面点估计结果，列于表 4。

而若采用传统立体像对空间前方交会方法，由于它仅考虑两张像片的信息且无法剔除粗差，故将产生C₆²=15种计算结果，具体数值见表 5。为对比不同方法的解算结果，将地面点真值及不同方法解算值的空间分布情况绘于图 2。

由表 4可以看出，经过多次迭代，含有粗差的像片权重逐渐减小到0，使其在平差中不起作用。本文方法可通过在平差过程中为观测值重新赋权来规避粗差，具有较好的抗差性。而传统立体像对空间前方交会方法，由于难以利用多余的观测信息，受单张像片观测精度的影响较大，抗差性能较差。从图 2、表 5中可以看出，由于粗差的存在，立体像对空间前方交会结果中至少有5个估值与地面点真值相比具有较大偏差。

此外，本文方法除了可剔除粗差的影响外，对于随机误差也有较好的平衡作用。由图 2、表 4、表 5可以看出，采用加权总体最小二乘的多像空间前方交会结果精度优于立体像对的解算精度，能取得更稳定的交会结果，该算例验证了本文方法的稳健性。

3.2 不同误差对算法的影响分析

为系统分析该算法的特点，需分别研究像点观测值误差及外方位线元素、角元素误差对解算结果的影响，以及它们的综合影响。以下将从误差对于交会精度的影响以及目标点交会偏差两方面进行分析。

3.2.1 误差对交会精度的影响分析

经试验，3种误差分别对本文算法交会精度的影响具有类似的特性。因此，为节省篇幅，文中将重点针对像点坐标误差对于交会精度的影响进行分析，同时直接给出外方位线元素、角元素分别对于算法交会精度的影响结果，而对其过程不再赘述。

为分析本文方法交会精度受像点坐标误差的影响情况，重复进行10 000次试验后，将99.9%置信度对应的精度值作为其交会精度。图 3(a)—图 3(f)分别为不同误差水平下(像点观测误差由0.5~3个像素(像素大小为4 μm)递增)，10 000次试验中小航高像片的定位精度频数直方图。

由图 3可以看出，当像点坐标的误差水平不同，本文方法具有不同的交会精度。像点坐标含有的误差越小，交会精度越高；像点坐标误差越大，交会结果的精度越低。其中，在像点坐标含有1.5个像素误差水平下，可实现1 dm的交会精度。考虑到本文交会方法可以综合考虑多张像片的信息，具有更好的稳健性，因此它对于观测误差的敏感程度较低。为说明这一特性，并与传统立体像对空间前方交会方法进行对比，在置信度为99.9%的情况下，可绘制其交会精度随着像点坐标误差(0~100像素误差，间隔0.5像素进行试验)的变化曲线，如图 4(a)所示。

同理，可绘制出另两组像片组合(大航高像片组合(0~0.5像素误差，间隔0.01像素进行试验)及混合航高像片组合(0~15像素误差，间隔0.5像素进行试验))的交会精度随着像点坐标误差的变化曲线，分别如图 4(b)、图 4(c)所示。

由图 4可以看出，传统立体像对与本文多像空间前方交会方法的交会精度均随着像点坐标误差的增大而近似线性下降。且同等误差水平条件下，各组像片采用本文方法具有明显较高的交会精度。

为给出本文算法交会精度随着像点坐标误差的变化率，可利用Origin的Linear Fit功能拟合出直线的斜率。类似地，由于交会精度随着外方位线元素(试验数据：小航高像片组合外方位线元素误差0~5 m，间隔0.1 m进行试验；大航高像片组合外方位线元素误差0~0.5 m，间隔0.01 m进行试验；混合航高像片组合0~5 m外方位线元素误差，间隔0.1 m进行试验)、角元素(试验数据：小航高像片组合外方位角元素误差0~600″，间隔10″进行试验；大航高像片组合外方位角元素误差0~5″，间隔0.01″进行试验；混合航高像片组合0~90″外方位角元素误差，间隔1″进行试验)误差也呈线性递减趋势，可同样给出直线斜率，列于表 6，其中像点坐标、线元素、角元素误差的横坐标单位分别为像素、米及角秒。

由表 6可以得出，采用本文多像空间前方交会方法，当小航高像片上像点坐标在1个像素误差水平时，可实现约7 cm的交会精度，而当误差为100个像素时，仍可实现7 m的交会精度；大航高像片上像点坐标在0.5个像素误差水平条件下，具有约18 m的交会精度，而混合航高像片的像点坐标在1个像素误差水平时，具有1 m左右的交会精度。小航高像片外方位线元素在1 m误差水平条件下，可实现约2.3 m的交会精度；大航高像片外方位线元素在0.5 m误差水平条件下，具有约23 m的交会精度，而混合航高像片外方位线元素在1 m误差水平条件下，具有约3 m的交会精度。小航高像片外方位角元素在1″误差水平条件下，具有约1 cm的交会精度；大航高像片外方位角元素在1″误差水平下，具有5 m左右的交会精度，而混合航高像片外方位线元素在1″误差水平条件下，具有1.5 dm左右的交会精度。

结合图 4及表 6可以得出，同等误差水平条件下，本文方法前方交会精度远远优于传统立体像对方法。利用本文方法进行空间前方交会的过程中，小航高像片的交会精度受误差影响最小，具有最高的交会精度，混合航高像片次之，而大航高像片对误差最为敏感。因此，在处理大航高像片数据时，可适当增加小航高像片信息，可有效降低其对误差的敏感性，进而提高交会精度。

3.2.2 误差对交会点位置偏差的影响分析

为充分分析本文算法受误差的影响特点，需研究地面点交会结果随着误差的偏移情况。参考文献[28]，设计以下仿真试验。

(1) 地面点坐标真值：在一定区域范围内选取均匀分布的100个地面点，其坐标已知。

(2) 误差大小：分别在像点坐标、外方位线元素、角元素上引入10像素、100 m、30″的误差。

(3) 地面点坐标解算值：引入步骤(2) 中的误差，利用本文方法进行多像空间前方交会，获得地面点位置的解算值。

(4) 交会结果分析：将地面点解算值与真值进行比较，绘制地面点解算值各分量的偏差示意图，如图 5—图 7所示。最后，分析不同误差对于交会结果的影响及它们的综合影响。

由图 5可以看出，小航高像片在单独加载像点坐标、外方位线元素、角元素误差时，其X、Y、Z 3分量之间的偏差曲线幅度较为相似，没有明显差别，但整体来说，线元素上的100 m误差将使得交会结果与真值之间产生最大的偏差。在3种噪声的综合影响下，X、Y、Z 3分量之间的偏差幅度相近，整体幅度水平接近于单独加载100 m线元素误差的情形。

由图 6可以看出，大航高像片在单独加载像点坐标、外方位线元素、角元素误差时，其X、Y分量的偏差幅度相近，明显较Z分量更为稳定，受误差的影响较小。整体来说，线元素上的100 m误差将使得交会结果与真值之间产生最大的偏差。在3种噪声的综合影响下，整体偏差幅度与单独加载100 m线元素误差的情形相近。

由图 7可以看出，相对于大航高像片组合，在混合航高像片中, Z分量对于各种误差的敏感程度大大降低，甚至在单独加载像点坐标和角元素误差的情况下，较X及Y分量偏差幅度更小。而在单独加载线元素误差时, Z分量表现的稍为敏感，偏差幅度较X及Y分量略大。

综合图 5—图 7可得，同等噪声水平分别加载到不同航高的像片组合中，会引起交会结果的偏差程度不同，其中小航高像片交会偏差最小，混合航高像片次之，大航高像片交会结果最差。在利用大航高进行多像空间前方交会时，目标点的Z分量对于误差较为敏感，以致交会结果产生较大偏差。此时，可以适当利用小航高像片信息，有效降低Z分量对误差的敏感程度，进而实现更高的交会精度。

3.3 对粗差的剔除效果分析

为系统分析本文算法的特性，需研究该方法对于粗差的剔除能力。经验证，当粗差分别存在于像点坐标或外方位元素中时，本文算法对于粗差的探测能力具有一些通性，此部分仅以粗差存在像点坐标中予以说明。同时，考虑到粗差出现的概率低且文中各组像片的数目较少，因此文中仅考虑每组像片数据中仅有一个粗差的情况。

首先，分别在小航高、大航高像片组合中第1张像片的横坐标上加入不同大小的粗差，同时使像点坐标误差逐渐增加，可绘制不同粗差的剔除率随着像点坐标误差变化的曲线，如图 8、图 9所示。

由图 8—图 9可以看出，同等粗差情况下，粗差剔除率随着像点坐标误差的增大而逐渐减小，即噪声越小，更容易发现粗差。在同等像点坐标误差水平时，大粗差更容易被剔除。小航高的像片在具有小于17个像素像点坐标误差时，可以完全剔除大于1 mm的粗差。在像点坐标误差小于100个像素时，小航高的像片可以不小于98%的概率剔除大于4 mm的粗差。大航高的像片在含有0.1个像素的像点坐标误差时，可以以95%的概率剔除1像素的粗差，在像点坐标0.5个像素以内的误差时，可以99%的概率剔除6个像素的粗差。

为了分析混合航高像片组合的粗差剔除能力，本文在混合航高像片组合中的低航高与高航高像片上分别加以同等大小粗差，并绘制粗差剔除率随着噪声的变化曲线，如图 10(a)—图 10(b)所示。

对比图 10(a)、图 10(b)可以看出，在混合航高像片组合中，当像点坐标误差为5个像素时，可以98%的概率剔除大航高像片上包含的0.1 mm的粗差；当像点坐标误差为15个像素时，可以接近100%的概率剔除大航高像片上包含的0.4 mm以上的粗差。而在同等像点坐标误差水平条件下，包含在小航高像片中的同等大小的粗差则难以被剔除。即混合航高像片组合对于大航高像片上的粗差更敏感，同等误差水平条件下，出现在大航高像片上的粗差更容易被发现。这是因为在同等大小粗差存在时，大航高像片上的粗差会引起它所决定的同名射线产生较大偏差，因此在选权迭代过程中会被赋予较小的权重。若小航高像片上的粗差继续增大超过一定数值时，即由它决定的同名射线发生更大偏差时，可使其很快被剔除。

4 结语

参考传统立体像对的空间前方交会原理，根据“目标点到同名射线距离的加权平方和最小”这一原则建立目标函数，实现多像空间前方交会，可得出如下结论：

(1) 根据总体最小二乘方法建立目标函数，对待求参数求一阶和二阶导数，当至少存在两条不平行的同名光线时，即可利用本文方法进行空间前方交会来确定地面点坐标。

(2) 采用选权迭代的方法，可以有效剔除观测数据中包含的粗差，规避粗差对于解算结果的干扰，同时为精度较高的观测数据赋予更高的权重，更利于得到精确的交会结果。

(3) 相对于立体像对，多张像片的空间前方交会能利用更多的观测信息，提高观测数据的利用率，可有效降低物点估计值Z分量对于大航高像片所含误差的敏感性，使估计结果具有更好的稳健性，一定程度上丰富了摄影测量理论。