同心球透镜与微相机拼接阵列复合结构的成像系统通过集成多个光电探测器并配合多路电子学处理单元，可以实现几亿至几十亿像素的成像，它能够解决传统光电系统大视场与高分辨率的相互制约的问题，并可以有效地降低了极大信息量光学系统的实现成本，可以满足航天遥感、航空侦查、边海防预警、城市安防等领域对于大视场、远距离的高像素成像要求，具有极其重要的战略需求和广阔的市场应用前景^{[1, 2, 3, 4, 5]}。目前，国内尚未见报道有科研机构开展十亿像素成像系统研制工作，而国外已经研制出多款亿级像素到数十亿像素不等的光电成像系统。其中，美国DUKE大学Marks等人^[6]研制出大视场超高像素成像系统(AWARE 10)。该系统由一个同心球透镜和98个微相机阵列拼接组成，最大输出像素为10亿，水平视场为120°，垂直视场为50°，焦距为35 mm，可以清晰分辨1 km外的人脸。此外，Brian Leininger等人^[7]在美国国防部先进项目研究局(DARPA)资助下研发了搭载于无人侦察机上的实时地面广域监视成像系统(ARGUS-IS)。该系统由4个完全相同的镜头和368个500万像素的探测器阵列拼接组成，可位于6 000 m高空拍摄单张覆盖40平方公里的照片，最大输出像素为18亿，视场角为60°，可以精确监控地面上只有15 cm大小的物件，智能探测和定位跟踪65个动态目标。

由于高分辨率微镜头拼接过程中涉及到巨大了数据处理量、并且拼接计算过程复杂度极高，如何保证图像拼接算法的准确性和实时性成为系统研制所要解决的核心问题之一。采用传统的串行图像拼接算法难以满足工程应用的实时性要求。近些年，随着多核CPU与图像处理单元(GPU)硬件技术的快速发展，图像拼接的并行加速算法已成为信息光学与机器视觉领域的一个重要研究方向。国内学者ZHOU等人^[8]提出了GPU加速SIFT算法具体实现方法。国内学者徐晶等人^[9]提出一种基于GPU的快速图片检测方案，采用SURF算法和SVM算法对图像进行特征提出和分类，使检测速度提升近10倍。国外学者Vincent Garcia等人^[10]全面阐述了基于GPU的KNN高维空间特征匹配算法，并给出了算法性能比对。国内学者郭一汉等人^[11]采用NCC匹配方法进行区域相似性度量，通过RANSAC算法估计图像身影变换模型，并行化实现该算法。

本文依据已设计完成的基于同心球透镜与微相机阵列复合结构的十亿像素瞬态成像系统，利用GPU强大的并行计算、浮点计算能力以及在内存资源分配管理方面的优势，提出了一种基于统一计算设备架构(CUDA)与先验信息相结合的自适应图像拼接并行加速算法。首先，利用高精度的四维标定平台在事先标定的成像重叠区域内，采用基于统一计算设备架构(Compute Unified Device Architecture，CUDA)的快速鲁棒特征(Speed-up Robust Feature,SURF)方法提取图像的特征点。然后，利用基本线性代数运算子程序(CUDA basic linear algebra subroutines,CUBLAS)加速基于随机KD-Tree索引的近似最近邻搜索算法，得到拼接成像重叠区域的初始的匹配点对。最后，提出一种改进的并行渐近式抽样一致性(IPROSAC)算法用于剔除误匹配点和空间变换矩阵参数估计。实验结果表明，该算法可以满足工程项目应用对图像处理关于实时性好、抗干扰能力强的要求。

成像系统主要由微相机阵列光电系统和高速图像处理软硬件集成平台两部分构成。

超高像素瞬态成像光电系统的第一级成像系统是一个同心球透镜光学系统，由它来将大视场的探测目标完全一致地成像在一个球形焦平面上。再由数百个微相机构成的第二级成像系统将球形第一焦平面二次成像到相应的CMOS探测器上，同时保持每个探测器所成图像都能有足够的重叠区域以进行计算拼接，而微相机阵列排布方式则类似于富勒烯(C60)上的碳原子分布^[12]，如图1所示。

图  1  微镜头阵列成像系统构成

Figure 1.  Micro-lens array imaging system

高速图像处理系统由CMOS微相机阵列成像分系统、图像实时监控分系统、图像高速存储分系统、图像实时拼接及控制分系统等五个部分构成，各分系统之间的关系如图2所示。

图  2  高速图像处理系统构成

Figure 2.  High speed image processing system

本文采用基于计算统一设备架构(Compute Unified Device Architecture，CUDA)与先验信息相结合的自适应图像拼接并行加速算法完成成像系统图像拼接过程。CUDA是一种由NVIDIA推出的通用并行计算架构^[13]，如图3所示。

图  3  计算统一设备架构原理图

Figure 3.  Schematic diagram of the CUDA

该架构使GPU能够解决复杂的计算问题，它包含了CUDA指令集架构(ISA)以及GPU内部的并行计算引擎。基于CUDA的图像拼接算法执行的过程分为Host/CPU端和Device/GPU端两部分，其中CPU主要负责算法逻辑与控制，GPU主要负责并行计算与浮点运算等。其中，CUDA并行计算过程包括两个层次的并行：一是线程网络(Grid)中的线程块(Block)间并行，二是线程块中的线程(Thread)间并行，这种高效的多线程并行数据处理模式使得面对计算量巨大、复杂度极高的任务时可以更加高效地利用GPU强大的并行计算资源。

通过高精度四维标定平台对微相机阵列拼接成像重叠区域的预标定，确定图像特征检测的空间范围进行特征检测提取，采用基于CUDA的SURF方法提出图像特征点，该过程包括构造尺度空间、检测提取特征点、确定特征点主方向、生成特征描述向量四个步骤^[14]，如图4所示。

图  4  基于CUDA的图像特征检测与描述

Figure 4.  Feature detection and description based CUDA

SURF特征检测是基于尺度空间理论以保证图像尺寸不变性，首先构建高斯金字塔尺度空间。这一过程与SIFT方法有很大的不同，SIFT采用高斯差分方法，而SURF首先采用Hessian矩阵行列式求得近似值图像，正因为存在这些不同才加快了SURF特征检测的速度。

式中，D_xx、D_yy、D_xy为方框滤波器模板与图像进行卷积运算后的值，0.9是Bay等人给出的一个经验值。通过矩阵行列式和特征值来筛选极值点，若Fast-Hessian矩阵的行列式为正，且两个特征值为一正一负，则保留该极值点。该过程只得到了一张Fast-Hessian矩阵行列式图，但是在金字塔图像中分为很多层，利用上述方法将构造金字塔尺寸空间过程并行化，如图5所示。其中，octave代表金字塔的顺序，这一逼近过程被分为k个线程进行并行计算，Hessian矩阵行列式的值在所有并行线程执行完成后被计算。

图  5  尺度空间并行构造过程

Figure 5.  Process of the scale space parallel construction

基于尺度空间的SURF特征检测并行加速过程如下所示^[15]：

Step1：采用非极大值抑制(non-maximum suppression，NMS)方法初步确定特征点，利用CUDA平台对尺度空间所有像素进行并行处理，每个线程负责尺度空间内的像素点是否为特征点的判断，根据Fast-Hessian矩阵行列式求得尺度图像的极值点，与其三维邻域26个点进行比较，如果它是这26个点中的最大值或最小值，则保留下来。

Step2：为了能够对候选特征点进行亚像素精确定位，采用三维线性插值法，同时筛掉那些小于一定阈值的点，得到精确定位的特征点。

(1)确定特征点主方向

特征点主方向是通过计算特征领域内的Harr小波响应值得到的。

Step1：构建一个以特征点为圆心、6 s(s为特征点的尺度)为半径的邻域，每个线程负责计算邻域内的点在水平方向和垂直方向Harr小波(Harr小波边长取4 s)响应值，并给这些响应值赋予不同的高斯权重系数。

Step2：标定60度扇形区域的区域间隔进行旋转，每个线程负责将扇形区域水平方向和垂直方向Harr小波响应值相加求和，遍历整个圆形区域。

Step3：将相加求和得到的最大值扇形区域的方向作为该特征点的主方向。

(2)生成特征描述子向量

由式(2)可知，SURF特征描述向量的维数是16×4＝64维，特征描述向量的构成如图6所示。

图  6  特征描述向量构成示意图

Figure 6.  Constitution sketch diagram of the feature description vector

SURF特征描述子属高维空间特征向量，特征匹配过程相当于高维空间中的最近邻搜索问题。本文采用由NVIDIA推出的基于CUDA平台的基础数值线性代数运算库BLAS的GPU实现-CUBLAS来加速基于随机KD-Tree索引的近似最近邻搜索算法，得到参考图像与待配准图像的初始的匹配点对。

基于随机KD-Tree的近似最近邻搜索算法原理如下^[16]：

x=(x₁,x₂,…,x₆₄)^T和y=(y₁,y₂,…,y₆₄)^T为待匹配的两个特征点SURF特征向量。

Step2：根据距离计算结果建立一个64维的二叉树索引数据结构。

Step3：针对每个特征点遍历整个64维二叉树，搜索近似最近邻结点和次近邻结点，如果最近邻距离和次近邻距离的比值小于预先设定的阈值(取值范围为0.2~0.4)，则认为最近邻结点即为该特征点的匹配点对；否则，舍弃该点，返回没有找到。

式中，N_R代表矩阵中第i行元素等于||r_i||²，第j列元素都等于||q_j||²。dist²(R,Q)在式(5)中的表达可以应用于CUBLAS。为了对内存资源合理配置，我们采用以下方法：首先，分别用m维和n维的向量存储N_R和N_Q。然后，式(5)中的加减运算采用传统CUDA内核来处理。

基于CUDA和CUBLAS的ANN搜索的并行计算如下：

Step1：利用CUDA计算向量N_R和N_Q；

Step2：利用CUBLAS计算m×n维矩阵A=-2R^TQ；

Step3：将矩阵A第i行的每个元素与向量N_R的第i个元素相加得到矩阵B；

Step4：采用并行插值排序法对矩阵B的各列进行排序，获得的矩阵称为C；

Step5：给矩阵C第j列的前k个元素与向量N_Q的第j个元素相加得到矩阵D；

Step6：计算矩阵D前k个元素的平方根，得到了k个最小的距离，获得的矩阵记作E；

Step7：提取矩阵E中最前面的k×n子矩阵，得到的矩阵就是所要求的k个最临近每个待匹配点处的距离矩阵。

(1)确定空间变换模型

(2)变换矩阵参数估计

本文将一种改进的渐进抽样一致性(PROSAC)算法^[18]并行化用于完成对空间变换模型的参数估计。传统PROSAC算法将匹配点对按照相似性程度进行了排序，并且按照排序结果进行采样，所以比随机抽样一致性(RANSAC)算法更快完成空间变换模型参数估计。但由于匹配点对距离较近，如果严格按照排序结果采样将导致参数估计陷入局部极值而过早地结束算法而没有找到最优的参数^[19]。因此采用了一种穿插方法对当前内点集合进行采样，避免上述问题的出现。改进的PROSAC算法并行化步骤如下：

Step1:初始化匹配点对集合S，利用并行排序算法对匹配点对欧式距离dist排序。

Step2：确定样本采用次数N，内点集合大小为n。为充分利用GPU资源，提高算法效率，采样次数N设定为64的整数倍。

Step3：每个Thread负责从集合S中随机按顺序抽取3个数据和第n个数据组成初始采样样本。

Step4：利用式(4)并行计算每个匹配点经矩阵M变换后到对应匹配点的欧式距离d，通过与门限阈值T比较，若d＜T，则将该点作为内点。

Step5：选取包含内点数目最多的一个内点集合(数目相同时，选择标准差较小的点集)再次并行计算变换矩阵M的参数。

其中，样本采样满足如下约束条件：

(1)必须保证采样样本中的4个匹配点对都在内点集合中的概率足够高，一般取值95%；

(2)3个内点不能共线，并行计算变换矩阵M的参数。

实验测试采用的操作系统为Windows 7 64位SP1，处理器为Intel Core i7-4790K@4.00GHz四核，内存为16 GB(DDR3 1 333 MHz)，显卡为NVIDIA GeForce GTX 980，开发工具为Visual Studio 2010 SP1和NVIDIA CUDA 6.5，程序设计语言C++。

随机选取相邻排布的微相机从中采集三幅图像进行拼接算法性能测试，该过程由图像特征检测提取、图像特征点匹配、空间变换矩阵参数估计组成。

(1)特征检测与提取

分别对无穷远距离具有相同场景、不同像素值(包括2M、5M、8M和13M)进行特征检测与提取实验测试，特征检测样本集如图7所示。

图  7  特征检测样本集

Figure 7.  Test samples of the image mosaic

采用本文提出的SURF并行特征检测算法提取图7中(a)、(b)两幅图像的特征点，特征检测结果如图8所示。

图  8  特征检测结果描述

Figure 8.  Description of the feature detection result

采用特征检测串行算法重复该实验过程，测试计算时间20次取平均值，对比结果见表1。

(2)图像特征匹配

采用本文提出的并行特征匹配算法对图8中的两幅图像特征进行匹配，匹配结果如图9所示。

图  9  特征匹配结果

Figure 9.  Result of the feature matching

采用特征匹配串行算法重复该实验过程，测试计算时间20次取平均值，对比结果见表2。

(3)空间变换矩阵参数估计

首先，采用本文提出的并行参数估计算法对图9特征匹配点对进行空间变换矩阵参数估计。

然后，采用与并行算法理论相同的参数估计串行算法重复该实验过程，测试计算时间20次取平均值，对比结果见表3。

采用本文提出的基于统一设备架构CUDA与先验信息相结合的自适应图像拼接并行加速算法，用于完成对图7微相机阵列相邻排布三幅图像的拼合，生成最终拼接结果如图10所示。

图  10  图像拼接结果示意图

Figure 10.  Result of the image mosaic

通过实验测试数据结果得知，随着图像分辨率的提高，采用本文提出的并行加速算法将表现出更好的性能，图像拼接计算时间与加速比见表4。

针对利用微相机阵列实现超分辨成像的图像拼接问题，提出一种基于统一设备架构(CUDA)与先验信息相结合的并行加速算法。首先，依据事先标定的成像重叠区域，采用基于CUDA的快速鲁棒特征(SURF)方法提取候选特征点集。然后，利用CUDA基本线性代数运算子程序(CUBLAS)加速基于随机KD-Tree索引的近似最近邻搜索算法，得到初始的匹配点对。最后，提出一种改进的并行渐近式抽样一致性(IPROSAC)算法对空间变换矩阵进行参数估计，得到空间几何变换关系。实验结果表明，该算法的图像拼接时间为287 ms，与单独采用CPU相比速度提高近30倍，满足工程应用对图像拼接实时性要求。