光束法区域网平差是线阵卫星遥感影像有地面控制条件下几何定位的重要方法。其基本思想是利用地面控制点坐标及其对应像点观测值等测量数据，同时求解成像传感器的外方位元素和地面点的三维坐标^[1]。当前，常规的光束法平差方法中均利用Gauss-Markov模型描述像点观测方程的函数模型，即假定函数模型已知、非随机，且观测值中仅含有随机误差。在此基础上引入外方位元素的虚拟观测方程以及控制点误差等方程，利用最小二乘法化迭代求解未知数^[2]。但是，在像点观测方程中，观测向量以及描述函数模型的系数矩阵均由观测值组成，即两者均包含随机误差。此时，Gauss-Markov模型并不严格成立，常规的光束法平差计算方法虽在应用中切实可行，但是理论上不够严密。

近些年来，国内外学者针对系数矩阵含有随机误差的函数模型进行了研究，引入了EIV(errors in variables)模型描述随机误差，并研究了基于该模型的总体最小二乘方法^{[3, 4, 5, 6, 7, 8]}，已在坐标转换^{[9, 10, 11]}、面阵影像空间后方交会^[12]以及影像去噪^[13]等方面得到了成功的应用。但是现有加权总体最小二乘难以直接应用于光束法平差计算，主要原因是光束法平差模型中仅有像点观测方程的系数矩阵包含随机误差，虚拟观测方程以及控制条件方程的设计矩阵均为单位阵，二者仍然满足Gauss-Markov条件。现存解决思路是将虚拟观测方程视为含约束条件的总体最小二乘^[14]，或者将EIV模型转化为方差含参数的Gauss-Markov模型利用自由极值法推导包含虚拟观测方程的总体最小二乘解^[16]，但是这两种解决思路均假定观测向量随机误差与系数矩阵随机误差为等精度观测量，且后者仅能处理虚拟观测为等精度观测的情形。能否合理给定各类观测值的权值，是线阵卫星遥感影像的光束法平差计算正确与否的重要前提，因而这两类思路难以直接应用于光束法平差计算。

本文拟引入EIV模型描述像点观测的函数模型，结合外方位元素虚拟观测方程与控制点误差方程，利用条件极值原理推导可处理多类虚拟观测值的总体最小二乘方法，基于该方法实现线阵卫星遥感影像的光束法区域网平差计算，之后分别使用模拟平差算例与两地区真实影像对本文方法与现有算法进行对比试验验证。

光束法平差的误差方程主要包括像点观测方程、外方位元素的虚拟观测方程以及控制点的误差方程。其数学模型可以表示为

点击浏览公式

常规光束法计算方法是直接利用最小二乘估计法化求解式(1)，但是式中的像点观测方程是由共线条件方程经线性化得到，且系数矩阵本身由外方位元素、地面控制点坐标和内方位元素共同求得，因此矩阵C、K_G和K_T中均包含随机误差，直接利用最小二乘进行光束法平差计算并不严密。此时，需要引入EIV模型对像点观测方程进行描述，并在此基础上构建基于总体最小二乘的光束法平差模型。

式(1)中的像点观测方程可利用Gauss-Markov模型表示为

点击浏览公式

点击浏览公式

式(1)中的外方位元素虚拟观测方程与控制条件方程虽然属两类不同方程，但是在形式上一致，可将两式联立并表示为关于X的方程

点击浏览公式

基于总体最小二乘的光束法平差模型中，共包括像点观测误差、虚拟观测误差和系数矩阵误差三类，根据极值条件原理可得如下目标优化函数

点击浏览公式

根据Kronecker积的性质可知^[15]

点击浏览公式

为使式(5)取得极小值，根据极值定理得到方程组

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

注意到式(17)左右两端均包含未知数向量，因而需要通过迭代求解$\hat X$。本文按照如下迭代计算步骤求解未知数：

(1) 给定权初值P₁⁽⁰⁾和P₂⁽⁰⁾，按最小二乘法计算X的初始值为

(2) 按方差分量估计法确定各类观测方程的方差^[17]，并确定P₂。其中，P₂为分块对角阵，由虚拟观测值类型确定其具体形式。

(3) 计算l⁽ⁱ⁾=L₁-A$\hat X$⁽ⁱ⁾，i为迭代次数；给定参数k；计算F⁽ⁱ⁾(k)。

(4) 计算

(5) 重复步骤(3)、步骤(4)，当‖$\hat X$⁽ⁱ⁺¹⁾-$\hat X$⁽ⁱ⁾‖＜ε时,停止迭代。

式(17)中，C^TP₂C是一个对角阵，对角线的非零元素是P₂的权值，因而C^TP₂C与广义岭估计矩阵等价，此时F(k)是一个增益系数，可用于调整P₂中各分量的大小以保证病态问题的稳定求解。

为验证本文方法的正确性与有效性，将利用模拟数据和真实线阵卫星遥感影像对本文方法进行验证。由于式(17)中参数需要确定参数k的取值，因而在试验中将文方法分成两种情况讨论：第1种按照先验信息确定参数，称为本文方法1；第2种从解算病态问题的角度出发，利用岭迹法确定参数k的值，称为本文方法2。

在模拟试验中，将分别对常规最小二乘算法(least squares algorithm,LS)^[18]、基于最小二乘的虚拟观测法(least squares based virtual observation method,LSVO)^[17]、文献^[16]中的基于总体最小二乘的虚拟观测法(total least squares based virtual observation method,TLSVO)方法和本文方法进行试验，用于对比几种方法的单次平差计算精度。

由于TLSVO方法仅能处理单类虚拟观测值的情况，且需要岭迹法确定参数，难以处理对权值依赖较高的光束法平差计算，因而在真实线阵卫星遥感影像试验中，仅将惯用光束法平差方法LSVO和本文方法1进行比较。

本试验用于验证本文方法在单次解算中的性能，试验算例仍然采用文献^[16]中的平差算例。该算例在给定系数矩阵、观测向量和未知数真值后，分别在系数矩阵和观测向量中加入随机误差，并进行平差计算，将平差解算得到的未知数估计量与未知数真值进行比较，以二者的偏差平方和作为计算精度的比较依据。原始系数矩阵A⁰、观测向量L⁰和真值X⁰分别为

系数矩阵和常数项中分别加入随机误差e₁～N(0,σ₀²I_n)和e_A～N(0,σ₀²I_mn)(σ₀²=0.1)后，得到的系数矩阵与观测向量为

矩阵N=A₂^TA的条件数cond(N)=2.08×10⁴，求解该方程属病态问题。分别采用最小二乘法、LSVO法、TLSVO法和本文方法对以上问题进行计算，并利用4种方法求得的估计量$\hat X$与真值X⁰的偏差平方和‖ΔX‖₂=‖$\hat X$-X⁰‖₂衡量解算结果与真实值的偏差。

注意到试验中观测向量和系数矩阵中加入的随机误差的方差大小相同，此时系数矩阵与观测向量应为等精度观测，因此在本方方法1中取系数矩阵观测量的权矩阵P_A=P₁，此时k=σ₀²/σ_p²。此外，为减少试验变量的差异，本文沿用文献^[16]中试验假设：实际观测向量和虚拟观测向量均为独立、等精度观测。

试验中利用岭迹法确定TLSVO算法中的准则参数λ，并令虚拟观测向量L₂=0。LSVO法与本文两类方法试验中，均利用验后方差分量估计方法确定实际观测量与虚拟观测量的方差，进而确定两类观测方程的权值。为减少主观因素对试验的影响，试验中TLSVO算法参数λ和本文方法2的参数k分别在各未知数变化趋于稳定的区间内取为

点击浏览公式

分别利用4种方法估计得到的及对应的偏差平方和的统计结果如表 1所示，本文方法2的岭迹如图 1如所示，图 2为‖ΔX‖₂随参数k变化的曲线图。

图 1 本文方法岭迹图 Fig. 1 Ridge trace of proposed algorithm

图选项

图 2 估计偏差随参数k变化曲线 Fig. 2 Change curve of bias with parameter k

图选项

由本算例的计算结果可知，引入虚拟观测方程后，LSVO算法的计算精度优于LS算法，主要因为虚拟观测方程的引入改善了原法方程的病态程度，从而提高了平差计算的稳定性与精度。本文方法2与TLSVO算法的平差精度相当，略高于本文方法1；3种方法的计算精度均优于LSVO算法，这是因为这3种方法均顾及了系数矩阵中存在的偶然误差，从而有更高的平差精度，这也体现出EIV模型在平差计算中的优势。本文方法1中F(k)虽由精度确定，但并非处理病态问题时的最优增益值，因而并不能确保处理病态问题时取得最高的估计精度。

为了进一步比较本文方法1、本文方法2与TLSVO算法之间的差异，在A₀⁰和L⁰中分别加入方差σ₀²=0.1,0.2,…,0.5的随机误差，共生成10组数据，每个方差对应两组数据，分别利用LSVO算法、TLSVO算法、本文方法1和本文方法2进行平差计算，所得偏差统计如表 2所示，待定参数同样按照式(18)确定。分别用LSVO算法、TLSVO算法和本文方法1计算得到的偏差平方和减去本文方法2得到的偏差平方和，并将10组试验所得到的偏差较差值绘制于图 3。图中任意离散点p(x,y)到x轴的距离表示按该点对应方法进行平差后，计算得到的偏差平方和比本文方法2的偏差平方和高y。若p(x,y)在x轴的上方，则该方法的求解精度低于本文方法2。图中，方块所对应曲线d₁=‖ΔX_LSVO‖-‖ΔX₂‖，三角对应曲线d₂=‖ΔX_TLSVO‖-‖ΔX₂‖，圆对应的曲线d₃=‖ΔX₁‖-‖ΔX₂‖，‖ΔX₁‖为本文方法1偏差，‖ΔX₂‖为本文方法2偏差。

图 3 4种方法在10组试验中的偏差差值 Fig. 3 The difference of biases from 4 algorithms in 10 experiments

图选项

从表 2与图 3的结果可知，在第3、第6、第7和第10组试验中，本文方法1和本文方法2的计算结果基本相当；本文方法1在大多数情况下精度优于TLSVO算法，部分情况下略低于TLSVO的解算精度；TLSVO算法部分情况下也能达到本文方法2的解算精度，但是结果并不稳定，在第7组试验中解算精度甚至低于LSVO算法；本文方法2在所有试验中均得到最高的参数估计精度，相比于本文方法1，该方法从处理病态问题的角度确定了最优的k值，在处理病态问题时有一定优势。

然而在实际光束法平差的迭代计算中，需快速选择待定参数，减少主观因素对结果的影响，如采用岭迹法确定各次迭代的待定参数，则会造成运算效率的显著下降。实际光束法平差中应依据先验知识确定参数k的取值，即使用本文方法1进行光束法平差计算，适当放宽精度要求，提高计算的效率与自动化水平。

本试验用于验证本文方法在真实线阵卫星遥感影像光束法平差中的计算精度，将分别采用河南某地(图 4(a))与中国西部某地(图 4(b))的SPOT-5 HRS影像进行验证。河南地区影像以平原地为主，包括部分山区，范围内共包含45个均匀分布的地面控制点。西部某地影像所摄地区以丘陵地形为主，由于地物特征不明显，该影像范围内仅包含10个地面控制点。两幅影像内控制点均由外业测量并通过内业人工量测得到，量测精度在1个像元左右。

图 4 SPOT-5 HRS影像控制点分布示意图 Fig. 4 Distribution of ground control points for SPOT-5 HRS images

图选项

SPOT-5卫星配备了高精度的定轨和定姿设备，因而可利用卫星辅助数据将成像时刻的瞬时外方位元素X_St、Y_St、Z_St、φ_t、ω_t和κ_t表示为^{[19, 20]}

点击浏览公式

在光束法平差迭代计算中，通常取P₁=I，并将单位权方差取为像点观测值的方差估计值₀^[17]，即$\hat \sigma _0^2 = \hat \sigma _p^2$。对于外方位元素的虚拟观测值和控制点观测值，通常依据参数的精度给定权初值进行预平差，并将第i类观测值在第n次迭代中的权值计算为^[21]

点击浏览公式

将两组影像数据分别按照式(1)列出像点观测方程、外方位元素虚拟观测方程和控制点的误差方程。由于缺少系数矩阵观测量精度的先验信息，试验中取k=1。分别利用惯用方法(LSVO)和本文方法进行光束法平差运算，所得试验结果如表 3和表 4所示，影像1中的3种试验方案中控制点分布依次为：4个角点、4个角点和1个中心点以及9个标准点位；影像2中控制点数量有限，因此仅比较控制点数为4个、5个和6个3种情况下两种方法光束法平差结果的差异。

由真实影像试验结果可知：

(1) 在河南地区影像试验中，按本文方法1进行光束法平差计算后，平面精度和高程精度均优于常规平差LSVO算法1 m左右；在西部地区影像试验中，本文方法1的平面精度高于LSVO算法0.5 m左右，高程精度优于LSVO算法1 m以上。这说明了本文所提出方法有更高的平差计算精度，与模拟数据试验的结论一致。

(2) 河南地区影像试验中，当4个角控制点参与平差时，两者平差结果相差最大，随着控制点数量的增加，两种方法平差结果的差距逐渐减小；在西部地区影像试验中，也有相同的结论。这表明控制点数量的增加对常规方法的精度提升更为显著，当控制点数量较少时，本文方法1的优势更为明显。

(3) 在3种控制点布设方案下，河南地区影像的光束法平差计算精度提升并不显著；当增设中心点后，检查点的平面精度略有下降，这表明该试验中的部分控制点精度较低，主要因为对控制点进行量测时，部分地面点难以辨识。尽管误差方程中引入了控制点的误差方程，但是试验未采用选权迭代等方法降低质量较差的控制点对试验结果的影响。如要进一步提高光束法平差计算的精度，需要引入抗差估计法对平差模型进行改进，并引入系统误差参数进行自检校光束法平差，这是需要进一步研究的内容。

本文顾及光束法平差中像点观测方程的误差特性，引入EIV模型描述该方程中的观测误差和系数矩阵误差，利用条件极值原理推导了基于总体最小二乘的线阵卫星遥感影像光束法平差计算方法。模拟试验结果表明，相比于传统方法，本文方法在单次平差计算中即可获得较高的平差计算精度。在真实影像试验中，本文方法的光束法平差计算精度优于常规算法，是线阵影像光束法平差计算的一条新的途径。然而，本文方法仍有不足之处需要完善，比如实际光束法平差计算中，参数k的确定存在一定的主观性，而该参数取值直接影响到平差结果，当从系数矩阵观测精度的角度或者改善法方程病态状况的角度出发，对k值的选取进行优化。此外，本文方法未考虑控制点中存在粗差的情况，且未引入系统参数进行自检校光束法平差试验。如要将本文方法进行实际应用，还需要提高算法的计算效率。这些都是今后研究的重点内容。