利用遥感影像特别是卫星影像可快速、高效地获取地理信息，但常规的遥感影像定位方法对控制点的依赖程度较高，控制点的数量与分布都影响着影像的定位精度。中国海岛礁数量众多、分布零散，控制点的获取效率低、成本高、风险大，因此，对于海岛礁影像，在无地面控制点(简称无控)条件下进行定位是首选。

文献[1]在严格成像几何模型的基础上，通过轨道外推技术提高了海岛礁卫星影像的无控定位精度；文献[2-3]针对航空影像，利用惯性测量单元(inertial measurement unit, IMU)/差分全球定位系统(differential global positioning system, DGPS)以及引入潮汐模型等对海岛礁区域在无控或稀少控制条件下的定位进行了尝试，得到一些有益的成果。

天绘一号和资源三号卫星是中国第一代高分辨率传输型立体测绘卫星，其摄影测量目标是无控测制1:50 000比例尺地形图，但由于其装备的星敏感器精度较低，无控条件下不经过平差的直接定位精度不能达到上述目标的要求^[4]。在这种情况下，区域网平差^[4-6]和利用已有DEM、正射影像等数据进行辅助^[7-9]是提高影像无控定位精度的有效方法，但前者计算过程通常较复杂，往往还受区域大小等的限制，后者则受制于辅助数据。

航天飞机雷达地形测绘任务(shuttle radar topography mission, SRTM) DEM等公开DEM具有较高的精度和全球陆地覆盖率^[10]，可为上述基于辅助数据的无控定位方案提供完备的基准数据。文献[9]将公开DEM用于辅助国产卫星的影像定位，实验表明，利用SRTM DEM辅助定位时，单景立体影像的无控定位精度可以满足1:50 000比例尺地形图的测制要求，为国产卫星影像的无控定位提供了较好的解决方案。该方法需要从影像中提取DEM，但海岛礁影像中存在或大或小的海域，其具有的纹理重复等特点不利于DEM的提取，当海域比例较大时，还会引起陆域DEM出现大量粗差甚至造成DEM提取失败，目前针对海岛礁区域的DEM提取方法也较少；此外，该方法也仅考虑了影像绝对误差的补偿。

本文针对海岛礁影像的无控定位问题进行研究时，将SRTM DEM作为辅助数据，在文献[9]所提方法的基础上，为克服海域影像纹理重复等对DEM提取的不利影响，将SRTM DEM也用于辅助DEM提取，并同时对多视影像中的相对误差进行补偿，提出一种SRTM DEM辅助的海岛礁卫星影像无控定位方法。

1 SRTM DEM

SRTM DEM是目前应用最广泛的全球公开DEM数据之一，其最新版本SRTM DEM V4.1由国际热带农业中心发布(http://srtm.csi.cgiar.org)，标称的绝对平面精度和绝对高程精度分别为±20 m、±16 m，置信度均为90%，空间分辨率为3″(≈90 m)，水平基准和高程基准均分别为1984世界大地测量坐标系(world geodetic system 1984, WGS84)和1996地球重力模型(earth gravitational model 1996, EGM96)。该DEM自发布以来就受到了研究者的广泛关注，目前，官方公布的和研究者的评价结果都表明其具有比先进星载热发射和反射辐射仪全球DEM (advanced spaceborne thermal emission and reflection radiometer global DEM, ASTER GDEM)更高的精度^[10]，故本文将该DEM的最新版本数据作为辅助定位的基准数据。

SRTM DEM具有较高的高程精度，且海域的高程统一为指定值，在利用其辅助提取DEM时，不仅可将其高精度的高程作为提取DEM的初始高程，还可用于大致区分开影像覆盖范围内陆域和海域对应的DEM格网点，提取DEM时则只需考虑确定陆域格网点的高程，从而能有效克服海域影像纹理重复等特点对影像匹配的不利影响。由于SRTM DEM的高程基准为EGM96模型，因此在对其加以利用前需将其换算为大地高，后文所提及的高程均为大地高。

2 有理函数模型

有理函数模型(rational function model, RFM)形式简单，与传感器无关，具有很好的通用性、高效性，且能很好地隐藏传感器的核心信息^[11]，得到了广泛的应用。RFM的数学模型为：

$ \left\{\begin{aligned} s &=F(\varphi, \lambda, h)=\frac{P_{1}(U, V, W)}{P_{2}(U, V, W)} \\ l &=G(\varphi, \lambda, h)=\frac{P_{3}(U, V, W)}{P_{4}(U, V, W)} \end{aligned}\right. $

(1)

式中，(φ, λ, h)是地面点坐标；(U, V, W)和(s, l)分别是正则化的地面坐标和像坐标；P_i(i=1, 2, 3, 4)为一般多项式，且多项式中3个分量各自的幂及幂之和均不超过3次^[12]。

由于RFM参数没有实际的物理意义，采用像方补偿方案对其进行改正时采用偏移补偿的方式^[11]：

$ \left\{\begin{array}{l}{S+\Delta S=S_{S} \cdot F(\varphi, \lambda, h)+S_{0}} \\ {L+\Delta L=L_{S} \cdot G(\varphi, \lambda, h)+L_{0}}\end{array}\right. $

(2)

式中，(S, L)为像坐标；(ΔS, ΔL)为像坐标系统误差的补偿值；(S_S, L_S)和(S₀, L₀)分别为像坐标的正则化缩放系数和平移参数。

采用仿射变换模型可对RFM中的系统误差进行有效的消除，本文就利用该模型对影像定位结果中的误差进行补偿，其形式为：

$ \left\{\begin{array}{l}{\Delta S=e_{0}+e_{1} \cdot S+e_{2} \cdot L} \\ {\Delta L=f_{0}+f_{1} \cdot S+f_{2} \cdot L}\end{array}\right. $

(3)

式中，e_i、f_i(i=1, 2, 3)为待求的RFM改正数。

3 基本原理

本文方法包括概略定位和精定位两个阶段，两个阶段均包括DEM提取和DEM匹配等步骤。为保证计算效率，概略定位阶段，首先选取少量样本点确定提取DEM的高程搜索范围，然后提取DEM与SRTM DEM进行匹配，从而将影像的定位精度提高到较高水平。在此过程中，若多视影像的交会条件较差，则进行自由网平差初步改善交会条件。精定位阶段，首先利用连接点进行自由网平差，对影像中的相对误差进行补偿，提高各视影像的交会条件；然后提取更精细的DEM进行匹配，进一步补偿影像中的绝对误差，实现影像的精确定位。本文方法的主要流程如图 1所示，下面分别对两个阶段DEM的提取及DEM匹配等主要步骤进行介绍。

3.1 概略定位DEM的提取

多视匹配提取DEM主要有两种策略：一种是逐像对匹配与匹配结果一致性检验相结合，另一种是通过物方高程变化形成的投影轨迹搜索对应匹配点^[13]。前者一般以核线立体像对为基础，不需要通过地面投影寻找对应像点；后者可靠性较高，尽管遮挡区域会形成一定的干扰，但当遮挡只存在于部分影像中时，仍可能得到正确的匹配结果^[14]。

本文采用文献[14]提出的利用地面铅垂线进行物方空间约束的多视铅垂线轨迹(multi-view vertical line locus, MVLL)法提取DEM，基本原理如图 2所示。首先在以近似高程Z₀为中心、容差ΔZ为半径确定的高程搜索范围[Z₀-ΔZ, Z₀+ΔZ]内根据高程间隔dZ可得到一系列的高程；然后分别投影至各视影像上并计算影像窗口之间的相关性，若相关性峰值满足阈值要求，则将峰值对应的高程赋给该格网点；反之，该点则未能确定有效的高程。下面对DEM提取的主要步骤进行介绍。

1) DEM范围及有效陆域点的确定。RFM参数中包括将像坐标和地面坐标进行正则化的参数，通过正则化可避免计算过程中参数数值量级差别过大而引入舍入误差。因此，可根据地面坐标的正则化参数大致确定影像的地面覆盖范围，此处将各视影像覆盖范围并集的外接矩形作为待提取DEM的范围。

设定DEM的分辨率R后，即可确定各格网点的平面坐标，进而可从SRTM DEM中内插各格网点的初始高程Z_ORI。但上述确定的矩形DEM范围内的格网点不全在影像的覆盖范围内，如4个角附近，这些点称为域外点；此外，还有分布在海域内的点，称为海域点。因此，需将这些格网点进行标记，不参与后续的高程值确定及DEM匹配等过程，后续处理只针对影像覆盖范围内的陆域格网点(称为有效陆域点)。根据内插所得高程即可较容易地区分出海域点；根据陆域点的平面坐标、初始高程以及下视影像RFM参数计算像坐标，在下视影像边界外的点则判断为域外点。

2) 高程搜索范围的确定。根据MVLL法的基本原理可知，高程搜索范围和高程间隔共同决定了DEM提取的效率和精度：高程搜索范围较大会使得地面铅垂线上的高程数量较多，影响计算效率，较小时则不一定能涵盖正确的高程值；当高程间隔较大时，最优相似测度确定的高程可能与实际高程相差较大，导致DEM的高程精度较低，较小时相邻高程值体现在像坐标上的差异太小，对DEM精度的提高没有太大的实际意义，还会降低计算效率。尽管此处将具有较高精度的SRTM DEM的高程作为各点的初始高程Z_ORI，但不同影像的直接定位精度差异较大，将初始高程Z_ORI直接作为近似高程Z₀时容差ΔZ不能太小。

为兼顾DEM提取的精度和效率，首先选取少量分布均匀的有效陆域点作为样本点，通过分析其高程误差的分布，为整个区域DEM的提取提供一个较小的更准确的搜索范围，具体步骤如下：

1) 在初始高程Z_ORI两侧较大的搜索范围[Z_ORI—ΔZ₁, Z_ORI+ΔZ₁]内，以稍大的高程间隔dZ₁通过MVLL法确定各样本点的高程Z_MVLL；

2) 若大部分样本点都通过MVLL法确定了高程，且Z_MVLL与Z_ORI之差dZ_MVLL-ORI分布集中，则确定上述高差集中分布的范围，该范围的中值和半径分别设置为高程改正值Z_GZ和容差ΔZ₂，整个区域DEM提取的高程搜索范围则为[Z_ORI+Z_GZ-ΔZ₂, Z_ORI+Z_GZ+ΔZ₂]；

3) 若上述条件不满足，则需首先通过自由网平差提高多视影像的交会条件，再按上述方法确定搜索范围。这是因为当多视影像的交会条件较差时，样本点处的铅垂线投影至各视影像上的投影线分布在不同位置，不会交于某点，因此难以正确提取各点的高程；而确定了高程的点也大多是粗差点，各点的高差dZ_MVLL-ORI分布也相对分散。

3.2 多视影像的自由网平差

在进行自由网平差时，首先通过影像匹配和选权迭代等处理确定w个分布均匀的连接点(文中w取300)，然后按照前视、下视、后视的顺序依次对各连接点列出如式(4)所示的误差方程：

$ \boldsymbol{V}=\boldsymbol{A} t+\boldsymbol{B} \boldsymbol{x}-\boldsymbol{l} $

(4)

式中，V是像点坐标的残差向量；t和x分别是仿射参数的改正数和连接点物方坐标的改正数；A和B分别为对应的系数矩阵；l是像坐标观测值和计算值之差。

对应的法方程为：

$ \left[\begin{array}{cc}{\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P A}} & {\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{B}} \\ {\boldsymbol{B}^{\mathrm{T}} \boldsymbol{P A}} & {\boldsymbol{B}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{B}}\end{array}\right]\left[\begin{array}{l}{\boldsymbol{t}} \\ {\boldsymbol{x}}\end{array}\right]=\left[\begin{array}{c}{\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P l}} \\ {\boldsymbol{B}^{\mathrm{T}} \boldsymbol{P l}}\end{array}\right] $

(5)

式中，P为权矩阵。

式(5)可简化为：

$ \left[\begin{array}{cc}{\boldsymbol{U}} & {\boldsymbol{W}} \\ {\boldsymbol{W}^{\mathrm{T}}} & {\boldsymbol{V}}\end{array}\right]\left[\begin{array}{l}{\boldsymbol{t}} \\ {\boldsymbol{x}}\end{array}\right]=\left[\begin{array}{l}{\boldsymbol{l}_{\boldsymbol{t}}} \\ {\boldsymbol{l}_{\boldsymbol{x}}}\end{array}\right] $

(6)

由式(6)可知，U和V都是块状对角阵，且t的参数数量远大于x的参数数量。此外，自由网平差缺少物方的控制，影像在空间中呈现为一种自由状态，解算时表现为法方程系数矩阵不满秩，即该平差是秩亏区域网平差^[15]。因此，此处采用块状高斯消元法^[16]，首先将式(6)中待求参数x消去，然后再依次对参数t和x进行求解：

$ \left\{\begin{array}{l}{\left(\boldsymbol{U}-\boldsymbol{W} \boldsymbol{V}^{-1} \boldsymbol{W}^{\mathrm{T}}\right) \boldsymbol{t}=\boldsymbol{l}_{\boldsymbol{t}}-\boldsymbol{W} \boldsymbol{V}^{-1} \boldsymbol{l}_{\boldsymbol{x}}} \\ {\boldsymbol{V} \boldsymbol{x}=\boldsymbol{l}_{\boldsymbol{x}}-\boldsymbol{W}^{\mathrm{T}} \boldsymbol{t}}\end{array}\right. $

(7)

这样可有效减小参数求解时法方程的阶数，保证计算效率。

3.3 DEM匹配及RFM改正数的求解

最小高差(least z-difference, LZD)法是基于DEM的2.5维结构提出的一种DEM匹配方法，对于规则格网数据的匹配，比其他方法具有更高的计算效率和匹配精度，此处采用文献[9]提出的能自适应探测及剔除DEM间的差异的引入截尾最小二乘(least trimmed square, LTS)估计的LZD法进行DEM匹配。

在确定各有效陆域点的高程时，由于噪声、云雾等的影响，部分点可能无法确定高程。通过分析可知，LZD法在确定对应点时从待匹配DEM格网点出发，在基准DEM中内插确定对应点高程，即当基准DEM是规则格网数据时，即可保证较高的计算效率，而对待匹配DEM的数据分布无特定要求。另一方面，在已有数据的基础上进行内插，对数据中的有效信息并无实质性增加，对匹配精度的提高也无实质意义，相反还会因内插以及参与匹配的点数增加导致计算效率降低。故此处直接将确定了高程的陆域格网点用于匹配，而不再进行内插得到无空洞的规则格网DEM。完成DEM匹配后，从两个DEM中选取分布均匀且高差较小的对应点求解或更新误差补偿模型的参数。

3.4 精定位DEM的提取

经过概略定位阶段对误差的补偿，影像的定位精度会有较大幅度的改善。本阶段则对上述误差补偿过程进行精化，以进一步提高影像的定位精度。同样地，先从有效陆域点中选取一定数量的分布均匀的点作为连接点，通过匹配得到各点在各视影像上的准确位置，然后一方面用于自由网平差进一步改善多视影像的交会条件；另一方面，各连接点前方交会确定物方坐标，再通过分析该高程与从SRTM DEM中内插得到的初始高程的差值Z_MVLL-ORI确定提取DEM的高程搜索范围。

与前一阶段确定高程搜索范围的原理不同，此处各格网点的高程改正值从上述各连接点的高差Z_MVLL-ORI中内插得到，容差设置为高差的3倍标准差，由此确定的各点处高程改正值各不相同，能较好地体现出定位模型在不同区域的高程误差，故所确定的高程搜索范围也更具针对性，有利于提取DEM时不同地形区域都能得到较准确的高程值。同时，利用较小的高程间隔dZ₂进行搜索。

4 实验结果与分析 4.1 实验数据

本文采用6景浙江舟山地区的天绘一号立体影像数据对本文方法的性能进行验证与分析，该测区以丘陵地形为主，高程范围约为0~800 m，每景立体影像包括1B级三线阵全色影像及对应的RFM参数，影像分辨率为5 m，数据获取时间为2014年1月至2015年2月，各景立体影像的生产方式和覆盖范围内陆域比例如表 1所示，其中生产方式包括先进产品生产系统(advanced products generation system, APGS)和等效框幅影像(equivalent frame photo, EFP)两种。实验区内共54个已知点，平面精度优于0.1 m，高程精度优于0.2 m，像点坐标采用人工转刺，精度优于1个像元。各景立体影像的覆盖范围和已知点的分布如图 3所示。

4.2 实验结果及分析

概略定位阶段的样本点数取500，高程搜索范围ΔZ₁取200 m，高程间隔dZ₁取10 m；精定位阶段样本点数取1 000，高程间隔dZ₂取5 m；从影像中提取的待匹配DEM(由于实验影像均为天绘一号影像，故后文简称TH-DEM)的分辨率R取0.000 7°；后续实验除R的取值有所差异，其他参数的设置均相同。本文实验均在Intel i7-4720HQ 2.60 GHz、16 GB RAM和MATLAB R2014a的软硬件环境下进行。

为验证本方法的可行性和有效性，首先对概略定位阶段的结果进行分析。各影像的直接定位结果如表 2所示，各影像直接定位的高程精度与概略定位阶段确定的高程搜索范围之间的关系如图 4所示。需要说明的是，为更清楚地展示两者之间的数值关系，图中高程精度区分正负，其正负根据检查点的误差分布确定。

由图 4可知，概略定位阶段利用少量样本点确定的高程初始改正值Z_GZ与影像直接定位的高程精度比较接近，高程搜索范围有效涵盖了高程精度。影像863-145和864-146对应的高程搜索范围较大，前者的平面误差较大，平面方向上的偏差将导致种子点高差分布相对分散；尽管后者的平面精度较高，但其覆盖区域的地形起伏相对较大，较小的平面偏差也可能导致较大的高差。另几景立体影像的高程搜索范围则相对较小。

这表明在概略定位阶段, 通过一定数量的样本点确定高程搜索范围是可行的，并且该高程搜索范围可根据影像覆盖范围内地形的起伏情况自适应确定，从而有效保证计算效率。

对本方法两个阶段的定位精度进行分析，并与有控定位的结果进行比较。对于有控定位方案，选择靠近影像角点且分布较均匀的3个已知点作为控制点，其余点作为检查点。各影像的定位结果如图 5(a)、5(b)所示，计算时间与陆域比例之间的关系如图 5(c)所示。

由图 5(a)、5(b)可知，相对于直接定位的精度，概略定位对精度有大幅度提高，可将影像的定位精度提高至较高的水平，具有不同陆域比例的影像均能得到较高的定位精度，表明提取DEM时利用SRTM DEM进行辅助可有效避免落水区域影像的不利影响。精定位能进一步提高定位精度，所得精度与有少量控制点参与定位的精度相当。采用本方法各景立体影像的平面和高程精度均分别优于6.1 m和5.2 m，能较好满足1:50 000比例尺地形图的精度要求。

经过两个阶段的纠正，各影像中检查点的误差分布盒状图如图 6所示，对于每一个盒子，红线表示各检查点误差的中位数，上下边缘分别表示上四分位数(q₃)和下四分位数(q₁)。将误差大于q₃+w(q₃-q₁)或小于q₁-w(q₃-q₁)的检查点视为粗差点，文中该线条的最大长度w设置为1.5，即将误差绝对值大于2.7σ的数据视为粗差。

由图 6可知，各检查点的误差分布比较集中，只存在个别的粗差点，即使最大误差的数值也较小，这表明本方法所得的定位结果不仅具有统计意义上的较高精度，影像覆盖范围内不同区域的精度也相对一致。

实验表明，本方法切实可行，能有效解决前文提及的具有大面积落水区域的海岛礁卫星影像的无控定位问题。

由图 5(c)可知，计算效率与陆域比例有很强的线性相关性，陆域比例越大，所需的计算时间越长。其中影像863-145的计算时间与陆域比例之比与其他影像有一定的差异，这是因为该影像的交会条件较差，概略定位阶段进行了自由网平差，而其陆域比例仅约4%，DEM提取及匹配等过程所需时间较少，自由网平差对总时间的影响较明显。而对于同样进行了自由网平差的影像864-147，因其陆域面积大，自由网平差对总时间的影响较小，故与另几景立体影像的计算效率与陆域比例之比较一致。对于上述两景由APGS生产方式得到的立体影像，概略定位阶段均需进行自由网平差，说明APGS方式所得RFM参数的精度低于由EFP所得的参数，这与实际情况相符。

为进一步分析待匹配DEM的分辨率对定位精度和计算效率的影响，分别从各影像中提取不同分辨率的TH-DEM用于匹配、定位，所得的定位精度和计算效率如图 7所示。

由图 7可知，从影像中提取不同分辨率的TH-DEM用于匹配以补偿误差时，不同陆域比例、不同生产方式的影像均能得到较稳定且较高的定位精度，受TH-DEM分辨率的影响较小，特别是高程精度基本保持不变，这表明通过两个阶段的误差补偿，影像中存在的相对误差和绝对误差都得到了充分的补偿。而文献[9]仅通过DEM匹配进行绝对定向，以对影像中的绝对误差进行补偿，并未考虑多视影像之间存在的相对误差，当多视影像的交会条件较差时，采用该方法进行定位所得的精度不稳定。影像863-145的平面精度随TH-DEM分辨率的变化有小幅度的波动，主要原因是该景立体影像的陆域比例太小，另几景则相对稳定。

对于计算效率，TH-DEM的分辨率越低，计算效率越高；当分辨率较高时，随着分辨率的降低，计算效率提高较快，且陆域比例越大，计算效率提高越明显；随着分辨率继续降低，计算效率提高速度逐渐变慢，计算时间趋于稳定。这是因为整个定位过程中涉及到DEM提取、匹配等过程，计算量较大，整体计算时间不可能太短。上述针对待匹配DEM分辨率与计算效率之间关系的分析，可为实际应用中的参数设置提供一定指导。

5 结语

本文提出一种SRTM DEM辅助的海岛礁卫星影像无地面控制点定位方法，利用该方法可对影像中的相对误差和绝对误差进行充分有效的补偿，对具有不同陆域比例、不同生产方式的天绘一号海岛礁影像均能得到较稳定且较高的定位精度，且基本不受待匹配DEM分辨率的影响，平面和高程精度分别优于6.2 m、5.2 m，能较好地满足1:50 000比例尺地形图的精度要求。

提取DEM时能根据影像覆盖范围内地形起伏情况自适应地确定较准确的高程搜索范围，从而保证较高的计算效率和DEM的提取精度。计算效率取决于陆域比例和待匹配DEM的分辨率，为保证计算效率，当待定位的海岛礁影像的陆域比例较大时，可提取比SRTM DEM分辨率略低的待匹配DEM；当陆域比例很小时(如不足10%)，待匹配DEM的分辨率应适当高于SRTM DEM或与之相当。