基于光電融合的目標三維結構反演

邢笑宇,霍超穎,殷紅成,王 靜

(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鏈碼法提取目標的輪廓角點所在位置;在光電融合方面引入圖論思想進行少量圖像的特征點關聯并利用逆投影原理實現特征點的三維重構,從而達到對散射中心與角點特征進行融合應用的目的。最后，基于組合體仿真數據進行了方法驗證,仿真試驗結果表明,該方法能夠同時反演出目標的典型結構及輪廓尺寸,從而提供出較單一雷達特征或可見光特征更為完整豐富的目標信息。