线阵传感器是当前遥感影像获取的主要手段。航空领域以ADS40为代表的线阵相机成功应用于各种比例尺航测成图^{[1, 2, 3]}。航天领域SPOT 5、IKONOS、QuickBird、WorldView、GeoEye、Pleiades以及我国的天绘、资源三号高分辨率卫星均为线阵传感器^{[4, 5, 6]}。当前，深空探测也是各航天国家的研究热点，欧美国家已经对月球、火星的地形进行了多次测绘，我国成功发射了嫦娥一号、嫦娥二号卫星，深空测绘相机也多采用线阵传感器^{[7, 8, 9]}。一方面是海量航空、航天遥感影像数据量呈现爆炸式增长，另一方面是海量原始遥感影像的数据处理能力明显不足，这就造成了大量遥感影像数据得不到及时处理，不能转换为可以应用的空间信息。如何有效、自动地对获取的大量遥感影像数据进行快速处理，已经成为亟待解决的问题^[10]。

线阵传感器的几何定位模型主要有严密几何模型、仿射变换几何模型以及有理多项式函数模型，其中，严密几何模型基于共线条件方程，理论最为严密，精度最高^{[10, 11, 12]}。近年来，有理多项式函数模型在IKONOS、QuickBird等商业卫星中得到了广泛应用，但是其RPC(rational polynomial coefficients)参数的求解过程仍然需要严密几何模型^[13]。另外，航空、航天线阵影像的几何检校也是基于严密几何模型^{[14, 15, 16]}。基于严密几何定位模型对线阵影像进行微分纠正是生成各级影像产品的基础，线阵影像微分纠正的难点是地面点反投影，核心问题是确定地面点对应的最佳扫描行。

线阵影像的地面点反投影基于严密几何模型，即共线条件方程建立的地面点、相机投影中心、像点3点共线关系，公式如下

点击浏览公式

下面先介绍传统的基于像方搜索的最佳扫描行定位算法，如图 1所示，给定地面点P(X,Y,Z)，P的反投影像点为p，ab表示线阵列CCD，传统的基于像方的搜索策略是利用反投影像点与线阵列CCD的距离d作为判断依据，即将最佳扫描行i的外方位元素代入共线方程求得的像点坐标p应该与线阵列CCD像元的检校坐标一致。显然，基于像方的最佳扫描行搜索计算过程涉及多次共线方程迭代，效率较低。

针对线阵影像最佳扫描行搜索策略，文献^[17]提出了像方二分法迭代算法，基于搜索窗口二分的思想快速缩小最佳扫描行搜索范围，该算法仍然基于像方共线方程迭代，计算量较大。文献^[18]分析了嫦娥一号激光测高数据与影像数据的不一致性，并提出了一种基于双线性变换快速确定搜索窗口的激光点反投影算法。文献^[19]提出了一种基于物方几何约束的最佳扫描行搜索方法，减少了共线方程迭代次数，计算效率与实用性均优于基于像方的搜索方法。分析现有算法可知，利用物方投影面几何约束是解决最佳扫描行搜索问题的有效途径。本文基于线阵影像严密几何模型提出了一种利用物方投影面几何约束的反投影快速算法，通过判断地面点与扫描行投影面的空间距离来定位最佳扫描行，这一点与文献^[19]的基本思想是一致的。由于线阵CCD各个像元并不严格位于一条直线上，线阵CCD的分段处理以及分段投影面的快速确定是影响基于物方几何约束最佳扫描行搜索效率的关键，也是本文算法与文献^[19]的主要区别。

图 1 线阵影像地面点反投影算法基本原理 Fig. 1 Basic principle of back projection for linear pushbroom imagery

图选项

本文利用线阵影像各扫描行投影面的几何约束关系快速定位最佳扫描行，算法思路如下：如图 2所示，线阵影像每一扫描行与投影中心可构成物方空间上的一个投影面，ab为第i条扫描行线阵列CCD首、末像元，S_i为该扫描行的投影中心，像方ab与投影中心S_i以及地面点AB可构成第i条扫描行的投影面。同样可计算像方cd与投影中心S_i+1以及地面点CD构成的第i+1条扫描行投影面。地面点P位于第i条扫描行投影面与第i+1条扫描行投影面之间，利用物方各扫描行投影面间的几何关系，通过判断地面点P与各扫描行投影面的距离，可快速定位与地面点P最邻近的扫描行i与i+1，再利用地面点P与最邻近扫描行间的距离精确内插出最佳扫描行。

基于物方投影面几何约束的地面点反投影算法也需要迭代进行，但是迭代计算是简单的空间解析几何运算，相比共线条件方程计算量大大减少，且迭代计算时可以使用各扫描行投影面间的距离进行快速估计，相比传统像方多次的共线方程迭代计算，本文算法仅需少量几次空间解析几何运算，因此最佳扫描行搜索效率得到大幅度提升。

图 2 基于投影面几何约束的最佳扫描行快速搜索 Fig. 2 Best scan line searching strategy based on geometric constraints of central perspective plane

图选项

基于物方投影面几何约束的线阵影像最佳扫描行搜索算法计算步骤如下：

步骤1 获取线阵影像严密几何定位模型所需的各扫描行时间、外方位元素、相机焦距、线阵列CCD探元位置等参数。

步骤2 以图 2中S_iAB构成的投影面为例，给定物方平均高程可由ab的像点坐标计算出AB的地面点坐标，S_i为第i条扫描行的外方位线元素，是已知值，设S_i、A、B 3点的地面点坐标分别表示为(x₁,y₁,z₁)、(x₂,y₂,z₂)与(x₃,y₃,z₃)，则由空间3点可确定一个空间平面，计算公式为

点击浏览公式

点击浏览公式

点击浏览公式

步骤3 根据相机投影中心以及各扫描行的投影面方程计算各扫描行投影面间的距离D_i、D₁表示第1条扫描行投影面与第2条扫描行投影面间的距离，以此类推。D_i用于对地面点P的最邻近扫描行进行快速估计。

步骤4 先判断地面点P与第1条扫描行投影面的距离d_i(i=1)，若d_i大于两个扫描行投影面间的距离D_i(i=1)，则利用下式对扫描行增量n进行估计，n=d_i/n_i。

步骤5 迭代判断P与第i条扫描行投影面的距离，直到d_ii并转至步骤6。

步骤6 利用与P最邻近的扫描行精确内插出最佳扫描行

点击浏览公式

受安置误差、镜头畸变等各种因素影响，线阵列CCD像元实际并不严格位于一条直线上，由线阵传感器相机检校文件得出的线阵列CCD像元一般有几个像素至几十个像素的偏移。图 3是ADS40相机前视、下视以及后视3条线阵的实际像元排列示意图，其中前、后视线阵直线偏差较大，而下视CCD直线性相对较好。

图 3 ADS40线阵CCD排列示意图 Fig. 3 CCD Arrangement of ADS40

图选项

可见，CCD线阵列像平面坐标通常表现为曲线形式，对应的线阵列CCD扫描行的物方投影面实际上是一个空间曲面。计算空间一点到空间曲面的距离比较复杂，为便于计算，对线阵列CCD进行分段处理，即由多个直线段代替曲线，这样曲面形式的物方投影面可由多个分段投影面代替，确定地面点对应的分段投影面后即可利用式(4)计算出点到扫描行投影面的空间距离。由于CCD线阵列整体仍然接近于直线，分段数目并不会太大。

如图 4所示，线阵列CCD分段处理步骤如下:

步骤1 设线阵CCD像元数量为n，各CCD像元P_i(i=1,2,…,N)的像点坐标可由相机检校文件得出。确定线阵列CCD首、末探元MN对应的直线方程

点击浏览公式

点击浏览公式

步骤2 求线阵列CCD各探元P_i到直线MN的距离d_i，P_i(x_i,y_i)到直线MN的距离为

点击浏览公式

步骤3 求距离d_i的最大值d_h，即d_h=max(d₁,d₂,…,d_n)。

步骤4 比较d_h与分段拟合阈值ε(例如0.5个像素)的大小，并计算开关量Q。

点击浏览公式

图 4 线阵列CCD分段示意图 Fig. 4 Line array CCD partitioned into several line segments

图选项

步骤5 根据Q决定是否提取中间特征点

(1) 如果Q=0，则可直接使用直线MN(M、N为特征点)代替曲线MN，转步骤6。

(2) 如果Q=1，则将d_h所对应的点P_j(x_j,y_j)抽出，作为中间特征点；然后连接新弦线MP_j，转第一步，并以MP_j代替MN，继续进行计算和判断。

步骤6 将所有提取出的中间特征点从起点M开始，顺序排列至终点N，得出分段拟合的特征点，根据这些特征点计算每一扫描行的各个分段投影面。

线阵CCD分段后，对应于某一影像扫描行的分段投影面数量为m，判断地面点至投影面距离时需要确定该点对应的分段投影面。如图 5所示，线阵列CCD实际像元排列为曲线abcde，曲线abcde可利用4条线段拟合，即ab、bc、cd、de，各个分段对应的地面点坐标分别为AB、BC、CD、DE，相机投影中心为S，此时的分段扫描行投影面为SAB、SBC、SCD、SDE，地面点P对应的像点为p，投影光线为pSP，像方p点位于bc段，对应于物方SBC投影面，SBC投影面的两条边界投影光线为SB、SC，P点在SBC的垂点为P’，由于P’与SBC共面，可以将SBC与P’所在的空间平面转换为二维平面，并计算出SB与SC在二维平面的直线方程

点击浏览公式

点击浏览公式

点击浏览公式

由图 5中可以看出，通过上述的判断准则，只有投影面SBC满足Q=1的条件，因此，可以确定分段投影面。

由于扫描行投影面方程仅需计算一次，本文的分段投影面确定方法主要计算量是式(11)中的d₁与d₂，即简单的空间几何计算，计算效率仍然是比较高的。为进一步提高计算速度，反投影迭代初始时可直接使用首、末CCD探元构成的近似扫描行投影面进行距离计算，当地面点至扫描行投影面的距离较小时，再使用分段投影面进行精确计算。从试验结果来看，利用近似扫描行投影面迭代3~4次，即可接近最佳扫描行，再利用分段投影面迭代一至两次即可定位与地面点P最邻近的扫描行投影面，整个迭代过程都是基于物方的简单的解析几何运算，且无需任何初始值，相比像方基于共线条件方程迭代的计算量小很多，从而实现最佳扫描行的快速、精确定位。

图 5 判断地面点对应的线阵列CCD分段投影面 Fig. 5 Determination of CCD line segment of central perspective plane

图选项

线阵影像的各级产品如核线影像与正射影像生成均基于微分纠正算法。微分纠正按照方法可分为正解法与反解法，实际多使用反解法^[20]。反解法微分纠正原理如图 6所示，由纠正影像上的地面点坐标反求其对应的原始影像像点坐标，并将该像点坐标的灰度值赋给纠正影像，由于求得的原始影像像点坐标一般不位于像素中心，一般采用双线性内插方法进行重采样^[21]。框幅式影像一幅图像内所有像素的外方位元素是一样的，由地面点计算像点坐标可直接利用共线方程计算得出，因此框幅式影像微分纠正效率相对较高，而线阵影像微分纠正的效率要比框幅式影像低很多，主要原因就是线阵影像地面点反投影需要迭代进行。

图 6 反解法微分纠正示意图 Fig. 6 Basic principle of differential rectification

图选项

当前线阵影像数据量普遍较大，ADS40相机一个条带单个通道影像数据量一般为几个GB，很多卫星影像的数据量也在GB级。可见，对线阵影像的最佳扫描行定位算法进行优化，对于提高线阵影像微分纠正算法效率是十分有必要的。

以正射影像生成为例，基于最佳扫描行快速搜索策略的线阵影像反解法微分纠正算法流程如下:

步骤1 构建线阵影像的严密几何构型，对CCD线阵列进行分段处理，计算各扫描行的近似投影面以及分段投影面。

步骤2 利用原始线阵影像以及DEM数据计算纠正影像的坐标范围。根据纠正影像坐标范围以及设定的像素间距大小s，计算纠正影像宽高。

步骤3 逐像素计算纠正影像上每一个像点对应的原始影像像点坐标，这个过程需要先进行一系列坐标转换，将纠正影像的坐标转换至摄影测量物方坐标系，然后是地面点反投影计算。

步骤4 对求出的原始图像像点坐标进行灰度重采样。

步骤5 逐像素进行线阵影像微分纠正计算。

线阵影像数据量较大，一次读取所有数据会超出内存上限，针对线阵影像的这一特点，算法实现时对影像进行分块处理。为减少磁盘IO和提高算法性能，本文从目标正射影像出发，将待生成的正射影像进行分块，并计算分块后的正射影像对应的原始图像范围，分块纠正时仅读取该范围内的原始图像数据，且各个分块影像可采用多线程方式进行并行处理，进一步提高处理速度。

选取航空ADS40以及火星快车HRSC线阵影像进行试验，其中，ADS40线阵分段数较多，而HRSC线阵直线性较好，分1段即可，试验数据基本情况见表 1。

试验环境基于Microsoft XP操作系统，采用Visual Studio 2008进行程序开发，计算机CPU配置为Intel酷睿i7系列(8×3.4GHz)，内存为3.24GB。随机选取100万个像点，给定地面点高程计算出每一个像点对应的地面点坐标，分别利用本文的最佳扫描行快速搜索算法和传统的基于像方迭代算法进行地面点反投影实验，并与已知点进行对比，结果见表 2。

由表 2数据可知，本文算法与基于像方的扫描行搜索算法定位精度相当，均达到0.001个像素的高精度定位，采用两种算法进行微分纠正的精度也是一样的，因此微分纠正实验不再比较精度，仅比较效率。表 3中T₁为基于像方迭代搜索最佳扫描行的微分纠正处理时间，T₂为利用本文算法的微分纠正处理时间，T₃为利用本文算法进行多线程并行处理所需时间，并发线程数量为8。

分析试验结果，得出以下几点结论：

(1) 在扫描行定位精度方面，本文算法与传统的基于像方的迭代搜索算法基本一致，均小于0.001个像素，即给定地面点可实现最佳扫描行的精确定位，这是保证线阵影像几何纠正精度的基础。

(2) 线阵CCD分段阈值为0.5个像素大小，表 3中ADS40前、后视影像分段数目较多，分别为20段与13段，下视影像直线性较好，分为两段，而火星快车HRSC直线性最好，分为1段即可。分析ADS40影像分段数量与处理时间可知，当线阵CCD需要分段处理时，分段数量虽然对处理效率有一定影响，但是由于本文算法初始迭代时使用了近似投影面减少计算量，分段数量对计算效率的影响并不十分明显。

(3) 分析最佳扫描行计算效率可知，由于基于像方的搜索算法一般需要十几次的共线方程迭代运算，而本文算法仅需要少量几次的物方空间解析几何运算，且在CCD分段处理以及分段投影面确定方面进行了优化，因此本文算法计算速度快，火星快车HRSC影像效率提升达到20倍，而ADS40影像效率提升达7~9倍。分析几何纠正效率可知，基于像方搜索算法进行几何纠正效率较低，不具备实际应用价值，而基于本文算法的微分纠正处理效率得到大幅度提升，ADS40影像效率提升6倍以上，火星快车HRSC影像效率提升13倍以上。

(4) 由于线阵影像数据量较大，加上线阵影像微分纠正算法本身的复杂性，若仅使用单线程方式进行几何纠正处理，即使使用本文的优化算法，ADS40一条全色波段影像纠正耗时也要1h左右。可见，大数据量线阵影像的几何纠正处理需要采用并行方式进行加速。经过并行处理后，在单机环境下，2GB大小的ADS40影像微分纠正仅需十几分钟即可完成，利用本文算法进行微分纠正较好地解决了线阵影像大数据量几何纠正的性能瓶颈问题。由于磁盘读写、数据传输、计算资源调度等开销，本文试验使用8个线程进行并行处理仅获得了3~4倍的加速比提升，若对并行算法进行优化，计算效率还有较大提升空间。

基于严密几何模型进行线阵影像的几何纠正是当前海量遥感影像处理的性能瓶颈，其核心问题是地面点对应最佳扫描行的快速定位。本文基于物方扫描行投影面几何约束提出了一种最佳扫描行快速搜索算法，无需要任何初始值，仅需要少量的空间解析几何运算即可快速定位最佳扫描行，具有精度高且计算速度快的特点。利用本文算法可较好地解决大数据量线阵影像的微分纠正问题，具有较强的实用价值。本文试验所用的并行算法仅在单机上采用多个线程，若采用集群或者GPU进行并行加速，则大数据量线阵影像的几何纠正效率将进一步提高，利用普通配置的计算机即可在几分钟内完成GB级数据量的几何纠正，线阵影像几何纠正的GPU与集群并行处理也是笔者下一步的研究重点。