基于光电融合的目标三维结构反演

邢笑宇,霍超颖,殷红成,王 静

(1.北京环境特性研究所电磁散射重点实验室,北京 100854;2.北京环境特性研究所光学辐射重点实验室,北京 100854)

0 引 言

随着雷达组网技术的蓬勃发展,多源信息融合研究受到了学者们的广泛关注,并被应用于目标特征提取与识别领域。目标特征反演本身属于不适定问题,而单部传感器获取的目标信息往往具有局限性、模糊性及非连续性,这将严重影响反演的准确度。相较之下,多源信息融合通过利用不同数据源之间的信息互补,能够为反演提供更多的有效信息,从而大大降低目标特征提取的不适定性。

可见光传感器与雷达传感器在目标信息获取方面具有天然互补性:可见光传感器成像通常较为清晰但易受气候、光照、探测距离等因素的影响；雷达传感器探测具有远距离、全天候、全天时的特点,但是成像效果往往不及可见光图像;可见光图像的成像面垂直于传感器视线,而雷达图像的成像面平行于传感器视线,两者恰好呈垂直关系;可见光图像对目标轮廓的捕捉较为精确,而通过雷达目标散射中心参数估计可以获得目标的结构细节信息[1-2]。因此,若将可见光图像与雷达图像加以融合,必定能对目标的特征挖掘与识别产生积极影响。

目前多源图像融合可分为3个层级[3]:像素级融合、特征级融合及决策级融合。像素级融合直接对原始图像进行配准和融合,其方法主要包括加权平均法[4]、塔式分解法[5]、小波变换法[6]等。像素级融合属于低层级的数据融合,对目标信息的损失最少,但其对单源图像的一致性要求很高,需要各图像的成像范围、质量及成像角度大致相同才能进行融合。特征级融合先提取多源图像特征再对所得特征进行综合利用,属于中间层融合,其对单源图像的一致性要求较低,信息损失也较为适中,常用的特征级融合方法有D-S推理法[7]、聚类分析法[8]、神经网络法[9]等。决策级融合先针对单源图像分别给出判别结果,再基于置信度对所得结果进行融合,这类融合对待融合图像的一致性要求最低,但是信息损失也最大,决策级融合的代表性算法包括贝叶斯推理[10]、模糊集合理论[11]、证据理论[12]等。

在上述3个融合层级中,像素级融合的相关研究最多,但是若要将其应用于可见光图像和雷达图像融合,则需要雷达与可见光传感器相对目标的视线为大体垂直的关系。对于空天目标观测来说,这个条件显然是难以满足的。决策级融合目前主要应用于大型识别系统,灵活性高,但是信息损失较大。相对而言,特征级融合对图像一致性要求及信息损失都较为适中,更适合于光电特征融合应用。

然而，当前利用特征级融合来处理光电数据的相关研究较少,且多集中在基于一致性较高图像的地物特征提取上[13-14]。本文针对空天目标光电协同探测问题,提出了一种在同传感器视角下融合逆合成孔径雷达(inverse synthetic aperture radar,ISAR)图像散射中心特征及光学图像轮廓角点特征的三维光电特征融合提取方法。首先,利用正交匹配追踪[15-16](orthogonal matching pursuit,OMP)和Freeman链码法[17-18]分别提取出序列ISAR图像中的二维散射中心特征以及光学图像中的角点特征。然后，引入图论中的二分图[19]思想实现序列二维图像中的散射中心匹配及角点匹配，并基于逆投影法[20]进行二维散射中心及角点特征的三维重构及融合。最后，综合利用两种特征进行目标的结构解译,相对于单源图像而言,该方法能够提取出更为完整精确的目标特征。

1 基于ISAR图像的目标散射中心特征提取

根据几何绕射理论[21](geometrical theory of diffraction,GTD),光学区目标的总散射场可以被近似认为是由若干个局部散射中心的散射贡献所合成的。因此,对于ISAR目标,可以采用二维GTD(two dimensioal GTD,2D-GTD)[1]模型来近似表示其高频散射响应。

1.1 2D-GTD模型

2D-GTD模型的表达式为

exp(-j4πf(xmcosφ+ymsinφ)/c)

(1)

式中：S为回波信号；f为雷达观测频率；fc为中心频率；φ为观测角度；c为光速；M为散射中心个数；Am为第m个散射中心的散射强度；xm和ym分别为目标散射中心在ISAR二维成像面上投影的横向位置及径向位置；αm为散射中心的类型参数,其取值为半整数形式且不同数值代表了不同的局部散射结构,如表1中所示。

表1 类型参数取值与局部典型结构之间的对应关系Table 1 Corresponding relationship between type parameter value and local typical structure

从式(1)可以看出,对2D-GTD模型进行参数估计,不仅可以得到目标主要散射体的位置信息,还可以反演其典型结构,这是从光学图像中难以获取的。

1.2 基于OMP方法的2D-GTD模型参数估计

为了保证2D-GTD模型的参数估计精度,本文基于散射中心稀疏分布的特点,引入稀疏思想来使散射中心的位置参数提取具有高分辨力。按照稀疏理论,2D-GTD模型可以被写成：

S=Φw

(2)

式中：测量向量S被定义为

S=[S(f1,φ1),S(f2,φ2),…,S(fN,φN)]T

(3)

S中的每个元素对应一次观测信号，N为总观测次数，Φ是字典矩阵:

Φ=[φ1，φ2，…,φP]

(4)

其中,φp是对应于第p个散射成分的原子(p=1,2,…,P),P为散射成分总数。

φp=(jf/fc)αpexp(-j4πf(xpcosφ+ypsinφ)/c)

(5)

式中：αp、xp和yp分别为第p个散射成分对应的类型、横向位置和径向位置参数。系数向量w为

w=[w1,w2,…,wP]T

(6)

由于2D-GTD模型的位置参数在连续区间内取值,稀疏字典维数较大,为了保证参数估计速度,这里采用OMP方法来进行稀疏参数估计。OMP[15-16]是一种较为经典的贪婪算法,能够基于局部最优准则实现快速稀疏求解,其主要步骤如下。

步骤 1输入:S,Φ,信号稀疏度K;

步骤 2初始化:

残差r=S,支撑集I=∅,迭代次数t=1;

步骤 3寻找局部最优解:

在每一次循环迭代过程中寻找满足目标函数的索引popt:

(7)

步骤 4更新支撑集和迭代次数:

I=I∪popt;t=t+1

步骤 5计算系数向量:

(8)

式中：ΦI为由Φ的第I列组成的矩阵。

步骤 6更新残差:

(9)

步骤 7收敛条件判断:

否则,重复步骤3～步骤7,直到满足收敛条件。

2 基于可见光图像的目标轮廓角点特征提取

2.1 Canny算子轮廓检测

相较于ISAR图像[22],可见光图像对目标轮廓的显示更为清晰,因此文中采用可见光图像来提取目标轮廓特征。

Canny算子[23-25]是一种常用的边缘轮廓检测方法,具有较好的鲁棒性及检测精度,其算法流程如图1所示。主要包括4个步骤:高斯滤波、梯度计算、非极大值抑制、双阈值检测与边缘连接。

图1 轮廓检测Fig.1 Contour detection

首先，采用高斯滤波[23]来对图像进行平滑,平滑后的图像用于梯度计算,其计算公式如下:

(10)

式中：M(x,y)为梯度幅值,反映了边缘的强度；θ(x,y)为梯度方向角,反映了边缘的方向；kx(x,y)和ky(x,y)分别为平滑后图像沿x和y方向的偏导。

在计算出图像中各像素点的梯度幅值及方向后,可以采用非极大值抑制的方法来找出图像中所有可能的边缘点。其基本思想为:将梯度方向角划分为0、1、2、3共4个扇区,如图2所示。在各像素点的3×3邻域内,沿同扇区方向判断该像素点的梯度幅值M(x,y)是否为极大值,若是,则该点为候选边缘点,梯度幅值保持不变;否则为非边缘点,梯度幅值设为0。

最后,采用双阈值法[23]从候选点中找出真正的边缘点并进行边缘连通,确定目标轮廓。

图2 梯度方向的扇区划分Fig.2 Sector division of gradient direction

2.2 基于Freeman链码分析的角点检测

相对于目标轮廓而言,轮廓角点是一种更为稳健的特征,因此本文采用轮廓角点特征进行光电融合，角点提取采用Freeman[17-18]链码法。

Freeman链码法通过单位向量序列来描述二值化轮廓的边界,其中较为常用的为8邻接链码法。该方法将两个相邻接轮廓像素点之间的方向向量归为8个链码值,如图3所示。

图3 Freeman 链码Fig.3 Freeman chain code

假设用li(li=0,1,…,7)表示由第i个轮廓像素点到第i+1个轮廓像素点的方向链码,这些链码li将组成一个有序而封闭的Freeman链码序列,以逆时针方向搜索链码,并将轮廓边缘差别码[17]定义为

(11)

由式(11)可知,差别码的可能取值包括0、±1、±2、±3和±4,接下来,以下准则可用于进行角点判断:

准则 1若di=0或di=4,则第i个轮廓点不是角点;

准则 2若|di|=3,则第i个轮廓点是角点;

准则 3若|di|=1或di=2,则第i个轮廓点是候选角点,需要通过曲率计算[26]来进一步进行判断。

3 基于二分图的特征点匹配算法

在提取了一系列的二维图像特征点后,即可对上述特征点进行匹配关联。现有的特征点匹配方法大多需要利用多幅目标连续运动图像来实现关联,如一维-高维散射映射图法[27]、Kalman滤波法[28-29]等。但是对于空天目标,受到信噪比、高速运动等因素的影响,连续获得多幅高清晰度的目标图像较为困难。为此,本文将图论思想引入到目标特征点匹配当中,采用二分图法[19]来对两幅或少量连续图像进行特征匹配。

首先,将特征点匹配问题转化为二分图问题,这里以两幅图像的特征点匹配为例进行介绍。将二分图记为G=(V1,V2,E),如图4所示。

图4 特征点匹配问题二分图Fig.4 Bipartite graph for feature point matching

其中，V1和V2为顶点集,分别对应第一幅图像中的N1个特征点和第二幅图像中的N2个特征点；E为二分图的边集,E={eij}(i=1,2,…,N1;j=1,2,…,N2),eij代表V1第i个顶点和V2第j个顶点之间的连线边,其取值为

eij=max(R)-rij

(12)

式中：rij为第一幅图像第i个特征点和第二幅图像第j个特征点之间的欧式距离；R为由全部rij组成的集合。

在特征点匹配二分图建立后,采用Kuhn-Munkres算法[30-31]求解二分图的最大权值作为最优匹配，算法流程如图5所示。

图5 Kuhn-Munkres算法流程图Fig.5 Flow chart of Kuhn-Munkres algorithm

算法具体步骤如下:

步骤 1构建赋值顶点集;

步骤 2构建相等子图;

步骤 3初始化匹配集；

步骤 4采用交替宽度优先搜索树法在相等子图中搜索可扩路;

步骤 5匹配扩展;

步骤 6判断当前匹配是否为最大匹配,若是则输出匹配结果,否则修改顶标值,重复步骤步骤2～步骤5,直到满足输出条件。

4 特征点三维重构与融合

最后通过特征点的三维位置重构来实现光电特征融合。图6给出了三维目标对雷达成像面及光学成像面的投影关系。假设o-xyz为目标本体坐标系,目标绕原点o旋转,其垂直于雷达视线方向的角速度分量单位向量为Ω,雷达视线单位向量为s,则垂直于Ω且包含s的平面为雷达成像面,记作α;垂直于s且包含Ω的平面为光学成像面,记作β。

图6 雷达成像及光学成像的投影关系Fig.6 Projection relationship between radar imaging and optical imaging

在平面α上建立ISAR成像坐标系uov,v轴与雷达中心视线s方向一致,u轴垂直于v轴。并定义光学成像坐标系uow,u轴方向与ISAR成像坐标系中相同,w轴沿Ω方向。假设P为目标上一点,位置为P,其在ISAR图像及可见光图像上的二维坐标位置可通过投影定理获得,对于ISAR图像有:

(13)

可见光图像中的投影位置为

(14)

若将P作为未知数,则公式皆为欠定方程,因此要想获得特征点的三维位置解,不管对于可见光图像还是ISAR图像,都需要联合至少两幅俯仰角不同的图像来求解。

以利用两幅图像进行目标三维特征点重构为例。目标在俯仰向上的转动可以等效为目标不动而传感器视线的俯仰角发生变化。假设两幅图像对应的传感器视线中心向量分别为s1和s2,对应的角速度向量分别为Ω1和Ω2,则根据式(13),散射中心的三维位置Ps可表示为

(15)

式中：vs1和vs2分别为两幅ISAR像中散射中心的径向位置；us1为第一幅ISAR像中散射中心的横向位置。

同理,基于两幅光学图像进行角点的三维位置重构,有

(16)

式中:Pc为角点三维位置；wc1和wc2分别为两幅光学图像中角点的纵向位置；uc1为第一幅光学图像中角点的横向位置。

在利用式(15)和式(16)完成特征点重建后,三维光学角点与三维散射中心都被归于同一目标坐标系下,从而实现了光电特征的融合。

5 仿真分析

为了验证提出算法的有效性,本文选取图7中的组合体作为研究对象,基于其电磁及可见光仿真数据来进行光电特征提取及融合试验。如图7所示,组合体由两个圆柱、两个平板以及两个连接杆组成,其坐标原点位于大圆柱的几何中心上,连接杆与大圆柱的交点处于圆柱的母线中点。表2中列出了组合体各部件的详细尺寸。

图7 组合体Fig.7 Combination

表2 组合体尺寸Table 2 Size of the combination

在组合体电磁仿真数据上叠加信噪比为15 dB的高斯白噪声作为回波信号,并分别基于雷达视线俯仰角θ1=85°和θ2=90°观测下的组合体回波信号进行二维散射中心特征提取,其中传感器视线的方位中心角保持在0°。图8中黑色点表示提取散射中心所在位置,表3列出了各散射中心参数的估计值。可见组合体在雷达视线上主要有4个散射中心,中间的两个散射中心类型参数为0.5,对应柱面反射机理,两端的散射中心类型参数为1,对应平面反射机理。

图8 组合体二维散射中心提取结果Fig.8 Extraction results of two-dimensional scattering center for combination

表3 组合体散射中心特征Table 3 Scattering centers features of combination

在与雷达同传感器视线的可见光图像中加入信噪比为15 dB的高斯噪声,形成测试图像,其角点特征提取结果如图9所示。该视角下组合体共有14个角点,图9中对角点进行了编号,相应的角点坐标如表4所示。

图9 组合体可见光图像(左上)及轮廓角点提取结果Fig.9 Combination visible image (on the top left)and contour corner extraction results

表4 组合体角点特征Table 4 Corner feature of the combination

接下来对散射中心和角点特征进行三维重构并融合,结果如图10和表5所示。

图10 组合体(灰色)及其光电特征Fig.10 Combination (gray)and its electro-optical features

表5 组合体光电特征Table 5 Electro-optical features of combination

图10中蓝色点代表散射中心,红色点代表角点,根据各特征点与组合体之间的位置关系可以看出,散射中心和角点的三维位置重构结果大体正确。

由表5可知,组合体的散射中心共有4个,根据类型参数可以推断,这4个散射中心分别对应两个柱面和两个平板,其几何中心位置可由位置参数确定。组合体角点共14个,分别对应着柱面和平板的边缘点,通过角点链接关系可计算出组合体的主要尺寸,如表6所示,其中边1～2代表角点1和角点2之间的边,其他边的表示方法与此类似。通过对比可知,表6中的尺寸计算值与目标的真实尺寸大致相同,尺寸反演结果正确。

表6 基于角点特征的组合体尺寸计算Table 6 Size calculation of combination based on corner feature

综上所述,根据表5中的散射中心参数可以确定组合体的主要结构,并能定位各结构的几何中心位置,但是无法确定各结构的轮廓尺寸。角点特征给出了组合体主要部件的轮廓尺寸信息,但缺少对各部件所属结构的认知。经光电融合后,如图10及表5所示,目标的特征点更为丰富,能够同时实现对目标的结构类型解译及轮廓尺寸反演,基本补全了目标朝向传感器视线一侧的几何信息,对特征的挖掘更为全面。

6 结 论

本文提出了一种综合利用同传感器视线下的雷达及光学数据进行目标特征反演融合的新方法。雷达方面采用散射中心特征实现目标基本结构的位置及类型判断;光学方面采用Canny算法及Freeman链码法提取目标的轮廓角点所在位置;在光电融合方面引入图论思想进行少量图像的特征点关联并利用逆投影原理实现特征点的三维重构,从而达到对散射中心与角点特征进行融合应用的目的。最后，基于组合体仿真数据进行了方法验证,仿真试验结果表明,该方法能够同时反演出目标的典型结构及轮廓尺寸,从而提供出较单一雷达特征或可见光特征更为完整丰富的目标信息。