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

(1.北京環(huán)境特性研究所電磁散射重點(diǎn)實(shí)驗(yàn)室,北京 100854;2.北京環(huán)境特性研究所光學(xué)輻射重點(diǎn)實(shí)驗(yàn)室,北京 100854)

0 引 言

隨著雷達(dá)組網(wǎng)技術(shù)的蓬勃發(fā)展,多源信息融合研究受到了學(xué)者們的廣泛關(guān)注,并被應(yīng)用于目標(biāo)特征提取與識別領(lǐng)域。目標(biāo)特征反演本身屬于不適定問題,而單部傳感器獲取的目標(biāo)信息往往具有局限性、模糊性及非連續(xù)性,這將嚴(yán)重影響反演的準(zhǔn)確度。相較之下,多源信息融合通過利用不同數(shù)據(jù)源之間的信息互補(bǔ),能夠?yàn)榉囱萏峁└嗟挠行畔?從而大大降低目標(biāo)特征提取的不適定性。

可見光傳感器與雷達(dá)傳感器在目標(biāo)信息獲取方面具有天然互補(bǔ)性:可見光傳感器成像通常較為清晰但易受氣候、光照、探測距離等因素的影響；雷達(dá)傳感器探測具有遠(yuǎn)距離、全天候、全天時的特點(diǎn),但是成像效果往往不及可見光圖像;可見光圖像的成像面垂直于傳感器視線,而雷達(dá)圖像的成像面平行于傳感器視線,兩者恰好呈垂直關(guān)系;可見光圖像對目標(biāo)輪廓的捕捉較為精確,而通過雷達(dá)目標(biāo)散射中心參數(shù)估計可以獲得目標(biāo)的結(jié)構(gòu)細(xì)節(jié)信息[1-2]。因此,若將可見光圖像與雷達(dá)圖像加以融合,必定能對目標(biāo)的特征挖掘與識別產(chǎn)生積極影響。

目前多源圖像融合可分為3個層級[3]:像素級融合、特征級融合及決策級融合。像素級融合直接對原始圖像進(jìn)行配準(zhǔn)和融合,其方法主要包括加權(quán)平均法[4]、塔式分解法[5]、小波變換法[6]等。像素級融合屬于低層級的數(shù)據(jù)融合,對目標(biāo)信息的損失最少,但其對單源圖像的一致性要求很高,需要各圖像的成像范圍、質(zhì)量及成像角度大致相同才能進(jìn)行融合。特征級融合先提取多源圖像特征再對所得特征進(jìn)行綜合利用,屬于中間層融合,其對單源圖像的一致性要求較低,信息損失也較為適中,常用的特征級融合方法有D-S推理法[7]、聚類分析法[8]、神經(jīng)網(wǎng)絡(luò)法[9]等。決策級融合先針對單源圖像分別給出判別結(jié)果,再基于置信度對所得結(jié)果進(jìn)行融合,這類融合對待融合圖像的一致性要求最低,但是信息損失也最大,決策級融合的代表性算法包括貝葉斯推理[10]、模糊集合理論[11]、證據(jù)理論[12]等。

在上述3個融合層級中,像素級融合的相關(guān)研究最多,但是若要將其應(yīng)用于可見光圖像和雷達(dá)圖像融合,則需要雷達(dá)與可見光傳感器相對目標(biāo)的視線為大體垂直的關(guān)系。對于空天目標(biāo)觀測來說,這個條件顯然是難以滿足的。決策級融合目前主要應(yīng)用于大型識別系統(tǒng),靈活性高,但是信息損失較大。相對而言,特征級融合對圖像一致性要求及信息損失都較為適中,更適合于光電特征融合應(yīng)用。

然而，當(dāng)前利用特征級融合來處理光電數(shù)據(jù)的相關(guān)研究較少,且多集中在基于一致性較高圖像的地物特征提取上[13-14]。本文針對空天目標(biāo)光電協(xié)同探測問題,提出了一種在同傳感器視角下融合逆合成孔徑雷達(dá)(inverse synthetic aperture radar,ISAR)圖像散射中心特征及光學(xué)圖像輪廓角點(diǎn)特征的三維光電特征融合提取方法。首先,利用正交匹配追蹤[15-16](orthogonal matching pursuit,OMP)和Freeman鏈碼法[17-18]分別提取出序列ISAR圖像中的二維散射中心特征以及光學(xué)圖像中的角點(diǎn)特征。然后，引入圖論中的二分圖[19]思想實(shí)現(xiàn)序列二維圖像中的散射中心匹配及角點(diǎn)匹配，并基于逆投影法[20]進(jìn)行二維散射中心及角點(diǎn)特征的三維重構(gòu)及融合。最后，綜合利用兩種特征進(jìn)行目標(biāo)的結(jié)構(gòu)解譯,相對于單源圖像而言,該方法能夠提取出更為完整精確的目標(biāo)特征。

1 基于ISAR圖像的目標(biāo)散射中心特征提取

根據(jù)幾何繞射理論[21](geometrical theory of diffraction,GTD),光學(xué)區(qū)目標(biāo)的總散射場可以被近似認(rèn)為是由若干個局部散射中心的散射貢獻(xiàn)所合成的。因此,對于ISAR目標(biāo),可以采用二維GTD(two dimensioal GTD,2D-GTD)[1]模型來近似表示其高頻散射響應(yīng)。

1.1 2D-GTD模型

2D-GTD模型的表達(dá)式為

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

(1)

式中：S為回波信號；f為雷達(dá)觀測頻率；fc為中心頻率；φ為觀測角度；c為光速；M為散射中心個數(shù)；Am為第m個散射中心的散射強(qiáng)度；xm和ym分別為目標(biāo)散射中心在ISAR二維成像面上投影的橫向位置及徑向位置；αm為散射中心的類型參數(shù),其取值為半整數(shù)形式且不同數(shù)值代表了不同的局部散射結(jié)構(gòu),如表1中所示。

表1 類型參數(shù)取值與局部典型結(jié)構(gòu)之間的對應(yīng)關(guān)系Table 1 Corresponding relationship between type parameter value and local typical structure

從式(1)可以看出,對2D-GTD模型進(jìn)行參數(shù)估計,不僅可以得到目標(biāo)主要散射體的位置信息,還可以反演其典型結(jié)構(gòu),這是從光學(xué)圖像中難以獲取的。

1.2 基于OMP方法的2D-GTD模型參數(shù)估計

為了保證2D-GTD模型的參數(shù)估計精度,本文基于散射中心稀疏分布的特點(diǎn),引入稀疏思想來使散射中心的位置參數(shù)提取具有高分辨力。按照稀疏理論,2D-GTD模型可以被寫成：

S=Φw

(2)

式中：測量向量S被定義為

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

(3)

S中的每個元素對應(yīng)一次觀測信號，N為總觀測次數(shù)，Φ是字典矩陣:

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

(4)

其中,φp是對應(yīng)于第p個散射成分的原子(p=1,2,…,P),P為散射成分總數(shù)。

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

(5)

式中：αp、xp和yp分別為第p個散射成分對應(yīng)的類型、橫向位置和徑向位置參數(shù)。系數(shù)向量w為

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

(6)

由于2D-GTD模型的位置參數(shù)在連續(xù)區(qū)間內(nèi)取值,稀疏字典維數(shù)較大,為了保證參數(shù)估計速度,這里采用OMP方法來進(jìn)行稀疏參數(shù)估計。OMP[15-16]是一種較為經(jīng)典的貪婪算法,能夠基于局部最優(yōu)準(zhǔn)則實(shí)現(xiàn)快速稀疏求解,其主要步驟如下。

步驟 1輸入:S,Φ,信號稀疏度K;

步驟 2初始化:

殘差r=S,支撐集I=?,迭代次數(shù)t=1;

步驟 3尋找局部最優(yōu)解:

在每一次循環(huán)迭代過程中尋找滿足目標(biāo)函數(shù)的索引popt:

(7)

步驟 4更新支撐集和迭代次數(shù):

I=I∪popt;t=t+1

步驟 5計算系數(shù)向量:

(8)

式中：ΦI為由Φ的第I列組成的矩陣。

步驟 6更新殘差:

(9)

步驟 7收斂條件判斷:

否則,重復(fù)步驟3～步驟7,直到滿足收斂條件。

2 基于可見光圖像的目標(biāo)輪廓角點(diǎn)特征提取

2.1 Canny算子輪廓檢測

相較于ISAR圖像[22],可見光圖像對目標(biāo)輪廓的顯示更為清晰,因此文中采用可見光圖像來提取目標(biāo)輪廓特征。

Canny算子[23-25]是一種常用的邊緣輪廓檢測方法,具有較好的魯棒性及檢測精度,其算法流程如圖1所示。主要包括4個步驟:高斯濾波、梯度計算、非極大值抑制、雙閾值檢測與邊緣連接。

圖1 輪廓檢測Fig.1 Contour detection

首先，采用高斯濾波[23]來對圖像進(jìn)行平滑,平滑后的圖像用于梯度計算,其計算公式如下:

(10)

式中：M(x,y)為梯度幅值,反映了邊緣的強(qiáng)度；θ(x,y)為梯度方向角,反映了邊緣的方向；kx(x,y)和ky(x,y)分別為平滑后圖像沿x和y方向的偏導(dǎo)。

在計算出圖像中各像素點(diǎn)的梯度幅值及方向后,可以采用非極大值抑制的方法來找出圖像中所有可能的邊緣點(diǎn)。其基本思想為:將梯度方向角劃分為0、1、2、3共4個扇區(qū),如圖2所示。在各像素點(diǎn)的3×3鄰域內(nèi),沿同扇區(qū)方向判斷該像素點(diǎn)的梯度幅值M(x,y)是否為極大值,若是,則該點(diǎn)為候選邊緣點(diǎn),梯度幅值保持不變;否則為非邊緣點(diǎn),梯度幅值設(shè)為0。

最后,采用雙閾值法[23]從候選點(diǎn)中找出真正的邊緣點(diǎn)并進(jìn)行邊緣連通,確定目標(biāo)輪廓。

圖2 梯度方向的扇區(qū)劃分Fig.2 Sector division of gradient direction

2.2 基于Freeman鏈碼分析的角點(diǎn)檢測

相對于目標(biāo)輪廓而言,輪廓角點(diǎn)是一種更為穩(wěn)健的特征,因此本文采用輪廓角點(diǎn)特征進(jìn)行光電融合，角點(diǎn)提取采用Freeman[17-18]鏈碼法。

Freeman鏈碼法通過單位向量序列來描述二值化輪廓的邊界,其中較為常用的為8鄰接鏈碼法。該方法將兩個相鄰接輪廓像素點(diǎn)之間的方向向量歸為8個鏈碼值,如圖3所示。

圖3 Freeman 鏈碼Fig.3 Freeman chain code

假設(shè)用li(li=0,1,…,7)表示由第i個輪廓像素點(diǎn)到第i+1個輪廓像素點(diǎn)的方向鏈碼,這些鏈碼li將組成一個有序而封閉的Freeman鏈碼序列,以逆時針方向搜索鏈碼,并將輪廓邊緣差別碼[17]定義為

(11)

由式(11)可知,差別碼的可能取值包括0、±1、±2、±3和±4,接下來,以下準(zhǔn)則可用于進(jìn)行角點(diǎn)判斷:

準(zhǔn)則 1若di=0或di=4,則第i個輪廓點(diǎn)不是角點(diǎn);

準(zhǔn)則 2若|di|=3,則第i個輪廓點(diǎn)是角點(diǎn);

準(zhǔn)則 3若|di|=1或di=2,則第i個輪廓點(diǎn)是候選角點(diǎn),需要通過曲率計算[26]來進(jìn)一步進(jìn)行判斷。

3 基于二分圖的特征點(diǎn)匹配算法

在提取了一系列的二維圖像特征點(diǎn)后,即可對上述特征點(diǎn)進(jìn)行匹配關(guān)聯(lián)?，F(xiàn)有的特征點(diǎn)匹配方法大多需要利用多幅目標(biāo)連續(xù)運(yùn)動圖像來實(shí)現(xiàn)關(guān)聯(lián),如一維-高維散射映射圖法[27]、Kalman濾波法[28-29]等。但是對于空天目標(biāo),受到信噪比、高速運(yùn)動等因素的影響,連續(xù)獲得多幅高清晰度的目標(biāo)圖像較為困難。為此,本文將圖論思想引入到目標(biāo)特征點(diǎn)匹配當(dāng)中,采用二分圖法[19]來對兩幅或少量連續(xù)圖像進(jìn)行特征匹配。

首先,將特征點(diǎn)匹配問題轉(zhuǎn)化為二分圖問題,這里以兩幅圖像的特征點(diǎn)匹配為例進(jìn)行介紹。將二分圖記為G=(V1,V2,E),如圖4所示。

圖4 特征點(diǎn)匹配問題二分圖Fig.4 Bipartite graph for feature point matching

其中，V1和V2為頂點(diǎn)集,分別對應(yīng)第一幅圖像中的N1個特征點(diǎn)和第二幅圖像中的N2個特征點(diǎn)；E為二分圖的邊集,E={eij}(i=1,2,…,N1;j=1,2,…,N2),eij代表V1第i個頂點(diǎn)和V2第j個頂點(diǎn)之間的連線邊,其取值為

eij=max(R)-rij

(12)

式中：rij為第一幅圖像第i個特征點(diǎn)和第二幅圖像第j個特征點(diǎn)之間的歐式距離；R為由全部rij組成的集合。

在特征點(diǎn)匹配二分圖建立后,采用Kuhn-Munkres算法[30-31]求解二分圖的最大權(quán)值作為最優(yōu)匹配，算法流程如圖5所示。

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

算法具體步驟如下:

步驟 1構(gòu)建賦值頂點(diǎn)集;

步驟 2構(gòu)建相等子圖;

步驟 3初始化匹配集；

步驟 4采用交替寬度優(yōu)先搜索樹法在相等子圖中搜索可擴(kuò)路;

步驟 5匹配擴(kuò)展;

步驟 6判斷當(dāng)前匹配是否為最大匹配,若是則輸出匹配結(jié)果,否則修改頂標(biāo)值,重復(fù)步驟步驟2～步驟5,直到滿足輸出條件。

4 特征點(diǎn)三維重構(gòu)與融合

最后通過特征點(diǎn)的三維位置重構(gòu)來實(shí)現(xiàn)光電特征融合。圖6給出了三維目標(biāo)對雷達(dá)成像面及光學(xué)成像面的投影關(guān)系。假設(shè)o-xyz為目標(biāo)本體坐標(biāo)系,目標(biāo)繞原點(diǎn)o旋轉(zhuǎn),其垂直于雷達(dá)視線方向的角速度分量單位向量為Ω,雷達(dá)視線單位向量為s,則垂直于Ω且包含s的平面為雷達(dá)成像面,記作α;垂直于s且包含Ω的平面為光學(xué)成像面,記作β。

圖6 雷達(dá)成像及光學(xué)成像的投影關(guān)系Fig.6 Projection relationship between radar imaging and optical imaging

在平面α上建立ISAR成像坐標(biāo)系uov,v軸與雷達(dá)中心視線s方向一致,u軸垂直于v軸。并定義光學(xué)成像坐標(biāo)系uow,u軸方向與ISAR成像坐標(biāo)系中相同,w軸沿Ω方向。假設(shè)P為目標(biāo)上一點(diǎn),位置為P,其在ISAR圖像及可見光圖像上的二維坐標(biāo)位置可通過投影定理獲得,對于ISAR圖像有:

(13)

可見光圖像中的投影位置為

(14)

若將P作為未知數(shù),則公式皆為欠定方程,因此要想獲得特征點(diǎn)的三維位置解,不管對于可見光圖像還是ISAR圖像,都需要聯(lián)合至少兩幅俯仰角不同的圖像來求解。

以利用兩幅圖像進(jìn)行目標(biāo)三維特征點(diǎn)重構(gòu)為例。目標(biāo)在俯仰向上的轉(zhuǎn)動可以等效為目標(biāo)不動而傳感器視線的俯仰角發(fā)生變化。假設(shè)兩幅圖像對應(yīng)的傳感器視線中心向量分別為s1和s2,對應(yīng)的角速度向量分別為Ω1和Ω2,則根據(jù)式(13),散射中心的三維位置Ps可表示為

(15)

式中：vs1和vs2分別為兩幅ISAR像中散射中心的徑向位置；us1為第一幅ISAR像中散射中心的橫向位置。

同理,基于兩幅光學(xué)圖像進(jìn)行角點(diǎn)的三維位置重構(gòu),有

(16)

式中:Pc為角點(diǎn)三維位置；wc1和wc2分別為兩幅光學(xué)圖像中角點(diǎn)的縱向位置；uc1為第一幅光學(xué)圖像中角點(diǎn)的橫向位置。

在利用式(15)和式(16)完成特征點(diǎn)重建后,三維光學(xué)角點(diǎn)與三維散射中心都被歸于同一目標(biāo)坐標(biāo)系下,從而實(shí)現(xiàn)了光電特征的融合。

5 仿真分析

為了驗(yàn)證提出算法的有效性,本文選取圖7中的組合體作為研究對象,基于其電磁及可見光仿真數(shù)據(jù)來進(jìn)行光電特征提取及融合試驗(yàn)。如圖7所示,組合體由兩個圓柱、兩個平板以及兩個連接桿組成,其坐標(biāo)原點(diǎn)位于大圓柱的幾何中心上,連接桿與大圓柱的交點(diǎn)處于圓柱的母線中點(diǎn)。表2中列出了組合體各部件的詳細(xì)尺寸。

圖7 組合體Fig.7 Combination

表2 組合體尺寸Table 2 Size of the combination

在組合體電磁仿真數(shù)據(jù)上疊加信噪比為15 dB的高斯白噪聲作為回波信號,并分別基于雷達(dá)視線俯仰角θ1=85°和θ2=90°觀測下的組合體回波信號進(jìn)行二維散射中心特征提取,其中傳感器視線的方位中心角保持在0°。圖8中黑色點(diǎn)表示提取散射中心所在位置,表3列出了各散射中心參數(shù)的估計值?？梢娊M合體在雷達(dá)視線上主要有4個散射中心,中間的兩個散射中心類型參數(shù)為0.5,對應(yīng)柱面反射機(jī)理,兩端的散射中心類型參數(shù)為1,對應(yīng)平面反射機(jī)理。

圖8 組合體二維散射中心提取結(jié)果Fig.8 Extraction results of two-dimensional scattering center for combination

表3 組合體散射中心特征Table 3 Scattering centers features of combination

在與雷達(dá)同傳感器視線的可見光圖像中加入信噪比為15 dB的高斯噪聲,形成測試圖像,其角點(diǎn)特征提取結(jié)果如圖9所示。該視角下組合體共有14個角點(diǎn),圖9中對角點(diǎn)進(jìn)行了編號,相應(yīng)的角點(diǎn)坐標(biāo)如表4所示。

圖9 組合體可見光圖像(左上)及輪廓角點(diǎn)提取結(jié)果Fig.9 Combination visible image (on the top left)and contour corner extraction results

表4 組合體角點(diǎn)特征Table 4 Corner feature of the combination

接下來對散射中心和角點(diǎn)特征進(jìn)行三維重構(gòu)并融合,結(jié)果如圖10和表5所示。

圖10 組合體(灰色)及其光電特征Fig.10 Combination (gray)and its electro-optical features

表5 組合體光電特征Table 5 Electro-optical features of combination

圖10中藍(lán)色點(diǎn)代表散射中心,紅色點(diǎn)代表角點(diǎn),根據(jù)各特征點(diǎn)與組合體之間的位置關(guān)系可以看出,散射中心和角點(diǎn)的三維位置重構(gòu)結(jié)果大體正確。

由表5可知,組合體的散射中心共有4個,根據(jù)類型參數(shù)可以推斷,這4個散射中心分別對應(yīng)兩個柱面和兩個平板,其幾何中心位置可由位置參數(shù)確定。組合體角點(diǎn)共14個,分別對應(yīng)著柱面和平板的邊緣點(diǎn),通過角點(diǎn)鏈接關(guān)系可計算出組合體的主要尺寸,如表6所示,其中邊1～2代表角點(diǎn)1和角點(diǎn)2之間的邊,其他邊的表示方法與此類似。通過對比可知,表6中的尺寸計算值與目標(biāo)的真實(shí)尺寸大致相同,尺寸反演結(jié)果正確。

表6 基于角點(diǎn)特征的組合體尺寸計算Table 6 Size calculation of combination based on corner feature

綜上所述,根據(jù)表5中的散射中心參數(shù)可以確定組合體的主要結(jié)構(gòu),并能定位各結(jié)構(gòu)的幾何中心位置,但是無法確定各結(jié)構(gòu)的輪廓尺寸。角點(diǎn)特征給出了組合體主要部件的輪廓尺寸信息,但缺少對各部件所屬結(jié)構(gòu)的認(rèn)知。經(jīng)光電融合后,如圖10及表5所示,目標(biāo)的特征點(diǎn)更為豐富,能夠同時實(shí)現(xiàn)對目標(biāo)的結(jié)構(gòu)類型解譯及輪廓尺寸反演,基本補(bǔ)全了目標(biāo)朝向傳感器視線一側(cè)的幾何信息,對特征的挖掘更為全面。

6 結(jié) 論

本文提出了一種綜合利用同傳感器視線下的雷達(dá)及光學(xué)數(shù)據(jù)進(jìn)行目標(biāo)特征反演融合的新方法。雷達(dá)方面采用散射中心特征實(shí)現(xiàn)目標(biāo)基本結(jié)構(gòu)的位置及類型判斷;光學(xué)方面采用Canny算法及Freeman鏈碼法提取目標(biāo)的輪廓角點(diǎn)所在位置;在光電融合方面引入圖論思想進(jìn)行少量圖像的特征點(diǎn)關(guān)聯(lián)并利用逆投影原理實(shí)現(xiàn)特征點(diǎn)的三維重構(gòu),從而達(dá)到對散射中心與角點(diǎn)特征進(jìn)行融合應(yīng)用的目的。最后，基于組合體仿真數(shù)據(jù)進(jìn)行了方法驗(yàn)證,仿真試驗(yàn)結(jié)果表明,該方法能夠同時反演出目標(biāo)的典型結(jié)構(gòu)及輪廓尺寸,從而提供出較單一雷達(dá)特征或可見光特征更為完整豐富的目標(biāo)信息。