張建民，艾偉，張學(xué)軍，陳路路，杜楚，員建廈，周云，杜丹

（1.中國(guó)電子科技集團(tuán)公司第五十四研究所，河北 石家莊 050081；2.河北省智能化信息感知與處理重點(diǎn)實(shí)驗(yàn)室，河北 石家莊 050081；3.陸裝駐石家莊地區(qū)第一軍代室，河北 石家莊 050081）

0 引言

探月雷達(dá)是我國(guó)嫦娥探月工程中對(duì)月表進(jìn)行實(shí)地登錄探測(cè)的重要手段[1-2]，其原理是通過(guò)天線發(fā)射的高頻電磁波脈沖對(duì)隱伏目標(biāo)體進(jìn)行探測(cè)。其中嫦娥三號(hào)探月雷達(dá)系統(tǒng)具有兩個(gè)通道：第一通道中心工作頻率為60 MHz，厚度上分辨率為米級(jí)，用于淺層月殼結(jié)構(gòu)探測(cè)；第二通道中心工作頻率為500 MHz，厚度上分辨率小于30 cm，用于探測(cè)月壤厚度和結(jié)構(gòu)[3-4]。

然而由于天線耦合、著陸器相關(guān)金屬構(gòu)件等工作條件的影響，使得月面下目標(biāo)的反射信號(hào)受到強(qiáng)烈干擾。為抑制雜波噪聲干擾，獲取更高分辨率的月表層位成像，許多信號(hào)處理方法被用于探月雷達(dá)信號(hào)處理，如反褶積方法[5]、帶通濾波[6]、振幅補(bǔ)償[7]、二維經(jīng)驗(yàn)?zāi)Ｊ椒纸鉃V波[8]、背景噪聲去除方法[9-10]、shearlet 變換[11]等，并取得了良好的處理效果，但這些算法核心構(gòu)成略顯復(fù)雜或存在噪聲殘留問(wèn)題，為此考慮選取算法構(gòu)成簡(jiǎn)單、性能優(yōu)越、不需頻域變換的形態(tài)學(xué)濾波方法來(lái)壓制探月雷達(dá)信號(hào)中的噪聲，提高數(shù)據(jù)的信噪比。

數(shù)學(xué)形態(tài)方法是一種基于積分幾何、隨機(jī)集合論等數(shù)學(xué)理論發(fā)展而來(lái)的非線性信號(hào)處理方法，算法簡(jiǎn)單，容易實(shí)現(xiàn)，且在圖像識(shí)別、噪聲壓制等方面得到了廣泛應(yīng)用。李杰等在2012 年采用數(shù)學(xué)形態(tài)學(xué)細(xì)化算法對(duì)圖像邊緣進(jìn)行細(xì)化處理，并取得了較好的結(jié)果[12]。王大瑩等利用形態(tài)學(xué)方法對(duì)激光雷達(dá)影像進(jìn)行建筑物邊緣提取，獲取了更連續(xù)的邊緣[13]。石玉敏采用基于數(shù)學(xué)形態(tài)學(xué)的船舶圖像質(zhì)量提升方法來(lái)獲得信噪比更高的船舶圖像[14]。胡愛軍等研究了開閉、閉開等組合形態(tài)算子對(duì)于濾波降噪、提高測(cè)試信號(hào)信噪比的有效性[15]。方文等提出了一種基于數(shù)學(xué)形態(tài)學(xué)的密集假目標(biāo)干擾抑制算法，實(shí)現(xiàn)了干擾環(huán)境下的目標(biāo)檢測(cè)[16]。柴曉飛等提出了一種基于數(shù)學(xué)形態(tài)學(xué)的氣象雜波抑制方法，較好地消除了氣象雜波對(duì)雷達(dá)系統(tǒng)的影響[17]。

因而，本文將數(shù)學(xué)形態(tài)學(xué)優(yōu)越的信息提取和噪聲壓制能力應(yīng)用到探月雷達(dá)信號(hào)的處理中，以提高探月雷達(dá)信號(hào)的信噪比和成像分辨率。首先，本文介紹了數(shù)學(xué)形態(tài)學(xué)濾波的基本原理；其次，針對(duì)探月雷達(dá)噪聲信號(hào)特點(diǎn)，提出一種基于自適應(yīng)數(shù)學(xué)形態(tài)學(xué)的信噪分離方法，該方法能夠基于群體的隨機(jī)搜索全局優(yōu)化方法，自動(dòng)選取最佳形態(tài)學(xué)結(jié)構(gòu)元素的長(zhǎng)度和振幅值；最后，利用該形態(tài)學(xué)濾波對(duì)嫦娥三號(hào)探月雷達(dá)第二通道數(shù)據(jù)進(jìn)行逐道處理，結(jié)果顯示，該方法能夠更有效地去除探月雷達(dá)信號(hào)中的噪聲及提高探月雷達(dá)成像的分辨率。

1 嫦娥三號(hào)探月雷達(dá)數(shù)據(jù)特征

嫦娥三號(hào)探月雷達(dá)具有CH-2A 和CH-2B 兩個(gè)通道，其中，第二通道收發(fā)天線類型為領(lǐng)結(jié)天線，安裝在玉兔巡視器底部，離地面約30 cm，工作中心頻率為500 MHz，頻率范圍為250—750 MHz[1,5]，數(shù)據(jù)采樣間隔為0.312 5 ns。搭載探月雷達(dá)的玉兔巡視器共行駛了114.8 m，玉兔巡視器路線圖如圖1 所示，本文選取其中CD 段的Channel-2B 數(shù)據(jù)進(jìn)行分析。

圖1 玉兔巡視器路線圖

由于數(shù)據(jù)采集時(shí)存在數(shù)據(jù)重復(fù)采集、時(shí)間延遲、振幅衰減、直達(dá)波干擾等情況[5]，在對(duì)探月雷達(dá)數(shù)據(jù)進(jìn)行進(jìn)一步分析之前，需要對(duì)探月雷達(dá)數(shù)據(jù)進(jìn)行預(yù)處理[11]。預(yù)處理后的第二通道CD 段數(shù)據(jù)如圖2所示，可以看出，預(yù)處理后的探月雷達(dá)信號(hào)中具有明顯的噪聲，以至于信號(hào)中的層位結(jié)構(gòu)不能被很好地識(shí)別出來(lái)，尤其是70 ns 以后的層位結(jié)構(gòu)非常不清晰（圖2（a）），嚴(yán)重影響了中深部數(shù)據(jù)的解釋。對(duì)預(yù)處理結(jié)果進(jìn)行時(shí)頻分析，可以看出該預(yù)處理數(shù)據(jù)中具有明顯的偏離主頻信號(hào)的低頻和高頻干擾，主要集中在100 MHz 和1 500 MHz，這些噪聲很可能來(lái)源于接收天線儀器本身和周圍采集環(huán)境的共同干擾。為壓制其中低頻和高頻噪聲的干擾，本文引入數(shù)學(xué)形態(tài)學(xué)濾波方法來(lái)抑制噪聲。

圖2 嫦娥三號(hào)探月雷達(dá)第二通道數(shù)據(jù)預(yù)處理結(jié)果

2 方法原理

2.1 數(shù)學(xué)形態(tài)學(xué)濾波基礎(chǔ)

數(shù)學(xué)形態(tài)學(xué)的基本思想是通過(guò)集合論方法來(lái)定量地描述目標(biāo)信號(hào)的幾何結(jié)構(gòu)，即利用預(yù)先定義好的結(jié)構(gòu)元素與信號(hào)的幾何特征進(jìn)行局部匹配或修正，同時(shí)保留目標(biāo)信號(hào)主要的形狀特征，以實(shí)現(xiàn)提取有用信息的目的。其中，結(jié)構(gòu)元素是數(shù)學(xué)形態(tài)學(xué)的基本要素，具有任意形狀和尺寸。

數(shù)學(xué)形態(tài)學(xué)的基本運(yùn)算包括腐蝕、膨脹、開運(yùn)算和閉運(yùn)算。設(shè)原始信號(hào)f(n)為定義在P={0,1,…,N-1}上的離散函數(shù)，結(jié)構(gòu)元素g(n)為定義在Q={0,1,…,M-1}上的離散函數(shù)，且N＞M，那么f(n)關(guān)于g(n) 膨脹和腐蝕運(yùn)算分別定義為：

式中，符號(hào)⊕和Θ 分別表示膨脹和腐蝕運(yùn)算。形態(tài)學(xué)上膨脹和腐蝕運(yùn)算實(shí)質(zhì)上是離散函數(shù)在結(jié)構(gòu)元素定義域內(nèi)的最大值和最小值濾波。其中，膨脹運(yùn)算表示一個(gè)擴(kuò)張過(guò)程，可以填平邊界不平滑的凹陷部分；腐蝕運(yùn)算表示一個(gè)收縮的過(guò)程，能夠剔除邊界不平滑的凸起部分。形態(tài)學(xué)閉、開運(yùn)算是在膨脹和腐蝕運(yùn)算基礎(chǔ)上組合而來(lái)的，分別定義如下：

式中，符號(hào)·和?分別表示形態(tài)閉和形態(tài)開運(yùn)算。其中形態(tài)學(xué)閉運(yùn)算是對(duì)同一結(jié)構(gòu)元素進(jìn)行先膨脹后腐蝕，開運(yùn)算是對(duì)同一結(jié)構(gòu)元素進(jìn)行先腐蝕后膨脹。開運(yùn)算可以剔除信號(hào)中的尖峰，抑制正脈沖噪聲，閉運(yùn)算可以補(bǔ)償谷底，抑制負(fù)脈沖的噪聲。相比膨脹和腐蝕運(yùn)算結(jié)果，開閉運(yùn)算可以更好地保留原始信號(hào)的主要特征。

2.2 廣義數(shù)學(xué)形態(tài)學(xué)濾波器的構(gòu)造

為了更好地濾除目標(biāo)信號(hào)中的正、負(fù)脈沖噪聲，在閉、開運(yùn)算的基礎(chǔ)上，采用相同的結(jié)構(gòu)元素，可構(gòu)建出形態(tài)學(xué)開-閉和形態(tài)學(xué)閉-開復(fù)合濾波器如下：

但由于開運(yùn)算的收縮性導(dǎo)致開-閉濾波器的輸出偏小，閉運(yùn)算的擴(kuò)張性導(dǎo)致閉-開濾波器的輸出偏大，所以為了取得較好的濾波效果，通過(guò)二者的平均值，構(gòu)造出新的濾波器如下：

式中，Mgf表示該形態(tài)學(xué)濾波器在結(jié)構(gòu)元素為g時(shí)對(duì)信號(hào)f的濾波結(jié)果。

形態(tài)濾波的質(zhì)量不僅取決于所選擇的形態(tài)變換，而且和選取的結(jié)構(gòu)元素也有較大的關(guān)系。結(jié)構(gòu)元素在形態(tài)學(xué)運(yùn)算中的作用類似于一般信號(hào)處理時(shí)的濾波窗口或參考模板，其形狀和尺寸都對(duì)形態(tài)學(xué)的運(yùn)算結(jié)果產(chǎn)生影響。常見的結(jié)構(gòu)元素形狀有直線形、正方形、拋物線形、三角形、正弦形等。本文選擇的結(jié)構(gòu)元素形狀為拋物線形，定義如下：

從式(8) 可以看出，拋物線形結(jié)構(gòu)元素結(jié)構(gòu)簡(jiǎn)單，其中參數(shù)A和L控制結(jié)構(gòu)元素的幅值和長(zhǎng)度。一般來(lái)講，利用大尺寸結(jié)構(gòu)元素的形態(tài)濾波器可以提取目標(biāo)信號(hào)的大致輪廓信息；小尺寸結(jié)構(gòu)元素的形態(tài)濾波器可以提取目標(biāo)信號(hào)中相對(duì)細(xì)節(jié)的信息。然而，對(duì)于探月雷信號(hào)既有低頻噪聲又含高頻噪聲的特點(diǎn)，采用單一結(jié)構(gòu)元素的形態(tài)濾波器類似于低通濾波器，不能呈現(xiàn)出信號(hào)中間某一段頻率范圍的信息。為了增強(qiáng)形態(tài)學(xué)濾波器的適用性，提取信號(hào)中有用的信息，通過(guò)如下思路對(duì)信號(hào)進(jìn)行形態(tài)學(xué)處理：

假設(shè)目標(biāo)信號(hào)f被結(jié)構(gòu)元素為gu的形態(tài)學(xué)濾波器處理得fu，那么根據(jù)式(7)可有：

如果信號(hào)fu被尺寸大于gu結(jié)構(gòu)元素gv的形態(tài)濾波器處理，可獲得fuv和RMMF，表示如下：

其中，RMMF代表提取信號(hào)中某一頻段范圍的信息。

2.3 形態(tài)學(xué)結(jié)構(gòu)元素的全局優(yōu)化方法

利用提出的數(shù)學(xué)形態(tài)學(xué)濾波方法獲取輸入信號(hào)中的有用頻段信息時(shí)，形態(tài)學(xué)結(jié)構(gòu)元素A和L的選取對(duì)處理結(jié)果有較大影響，因此如何確定結(jié)構(gòu)元素尺寸的參數(shù)，是應(yīng)用數(shù)學(xué)形態(tài)學(xué)濾波方法實(shí)現(xiàn)有效壓制信號(hào)中噪聲的關(guān)鍵。本文利用基于群體的隨機(jī)搜索全局優(yōu)化方法求解形態(tài)學(xué)濾波方法的最優(yōu)A和L值，以構(gòu)建自適應(yīng)形態(tài)學(xué)濾波方法。

粒子群隨機(jī)搜索全局優(yōu)化方法具有收斂速度快、精度高、操作性強(qiáng)等優(yōu)點(diǎn)，在許多優(yōu)化問(wèn)題中得到成功應(yīng)用[18]。粒子群優(yōu)化算法首先在給定的空間內(nèi)隨機(jī)對(duì)粒子群進(jìn)行初始化，并且每個(gè)粒子都有一個(gè)由適應(yīng)度函數(shù)決定的適應(yīng)值。每個(gè)粒子根據(jù)個(gè)體經(jīng)驗(yàn)和鄰居經(jīng)驗(yàn)來(lái)決定飛行的距離和方向，從而粒子群就隨當(dāng)前的個(gè)體極值和全局極值，一代代在解空間中搜索出最優(yōu)解。

設(shè)在一個(gè)D維空間中，由M個(gè)粒子組成的種群為X=(X1,X2,...,XM)，第i個(gè)粒子的位置為。將Xi代入適應(yīng)度函數(shù)中計(jì)算出適應(yīng)度值及粒子位置。第i個(gè)粒子的速度為，其中粒子i的個(gè)體極值為,種群的全局極值為。粒子通過(guò)兩個(gè)極值迭代更新自身的速度和位置，如下：

式中，i=1,2,…,M，d=1,2,…,D，M和D分別為種群的大小和待優(yōu)化問(wèn)題的維數(shù)；w是慣性權(quán)重；k是當(dāng)前迭代次數(shù)；c1、c2是加速常數(shù)；r1、r2是隨機(jī)數(shù)，范圍為(0,1)。

本文采用包絡(luò)熵Ep作為適應(yīng)度函數(shù)[19]，來(lái)評(píng)價(jià)粒子位置的好壞程度。其原因?yàn)椋喊j(luò)熵Ep可定量表示原始信號(hào)的稀疏性，Ep值越大，則信號(hào)稀疏性越弱，即含有越多噪聲；反之，Ep值越小，則信號(hào)稀疏性越強(qiáng)，即含有噪聲越少，當(dāng)噪聲最少時(shí)，粒子所處位置最好，即對(duì)應(yīng)結(jié)構(gòu)元素的A和L最優(yōu)。信號(hào)S 通過(guò)Hilbert 解調(diào)后，所得包絡(luò)信號(hào)sig(i)(i=1,2,…,N)的包絡(luò)熵Ep計(jì)算公式如下：

式中，pi為包絡(luò)信號(hào)sig(i)歸一化形式。

形態(tài)學(xué)濾波方法最佳結(jié)構(gòu)元素振幅和長(zhǎng)度的優(yōu)化搜索過(guò)程如下：

①對(duì)粒子群優(yōu)化方法參數(shù)進(jìn)行初始設(shè)定，將相應(yīng)的包絡(luò)熵Ep確定為適應(yīng)度函數(shù)；

②隨機(jī)產(chǎn)生粒子種群中一些粒子的位置及其移動(dòng)速度；

③對(duì)信號(hào)在不同位置下進(jìn)行形態(tài)學(xué)濾波處理，計(jì)算每個(gè)粒子位置相應(yīng)的Ep值；

④對(duì)比Ep值大小，利用式(12)、(13)更新粒子的移動(dòng)速度和位置；

⑤循環(huán)步驟③～④，當(dāng)?shù)螖?shù)達(dá)到閾值后結(jié)束循環(huán)，并輸出最佳粒子位置。

3 實(shí)測(cè)數(shù)據(jù)測(cè)試與性能分析

為獲得更加清晰的探月雷達(dá)成像分辨率，利用以上提出的自適應(yīng)形態(tài)學(xué)濾波方法，對(duì)圖2（a）信號(hào)中噪聲部分進(jìn)行壓制，提取中間頻率范圍有用信息。其中，為加快運(yùn)算速度以及保證結(jié)構(gòu)元素幅值作用的一致性，采用結(jié)構(gòu)元素長(zhǎng)度不同、幅值相同的兩個(gè)結(jié)構(gòu)元素，來(lái)構(gòu)造式(11)中的形態(tài)學(xué)濾波器。對(duì)兩個(gè)結(jié)構(gòu)元素尺寸自動(dòng)尋優(yōu)的粒子群算法參數(shù)設(shè)置如表1 所示，尋優(yōu)過(guò)程中極小包絡(luò)熵值隨種群進(jìn)化代數(shù)的變化如圖3 所示，極小包絡(luò)熵最小時(shí)，完成優(yōu)化過(guò)程，對(duì)應(yīng)的兩個(gè)結(jié)構(gòu)元素的長(zhǎng)度分別為5 和6，幅值為3。

表1 粒子群算法各項(xiàng)參數(shù)

圖3 極小包絡(luò)熵值隨種群進(jìn)化代數(shù)的變化

利用帶通濾波、EMD（Empirical Mode Decomposition，經(jīng)驗(yàn)?zāi)B(tài)分解）方法以及自適應(yīng)形態(tài)學(xué)濾波方法，分別對(duì)預(yù)處理后的探月雷達(dá)數(shù)據(jù)進(jìn)行處理，結(jié)果如圖4、圖5、圖6 所示。利用帶通濾波處理后的雷達(dá)成像結(jié)果中有明顯的條帶狀噪聲以及70 ns 后的層位結(jié)構(gòu)難以識(shí)別（圖4（a）），時(shí)頻譜特征顯示，帶通濾波能夠很好地去除信號(hào)中的高頻噪聲以及部分低頻噪聲，但80 ns 后仍殘留較多的低頻噪聲（圖4（b））。EMD 方法處理后的結(jié)果顯示，EMD 方法雖然能夠較好地去除信號(hào)中的低頻和高頻干擾，但雷達(dá)數(shù)據(jù)成像結(jié)果中仍存在比較明顯的條帶狀噪聲（圖5），影響數(shù)據(jù)的解釋和分析。圖6 顯示了利用自適應(yīng)形態(tài)學(xué)方法處理的結(jié)果，可以看出自適應(yīng)形態(tài)學(xué)方法相比帶通濾波以及EMD方法具有更強(qiáng)的噪聲壓制能力，能夠很好地去除信號(hào)中低頻和高頻噪聲，保留預(yù)處理后數(shù)據(jù)中的有用信息，使反射的回波信號(hào)更加連續(xù)清晰，同相軸特征明顯。同時(shí)，利用式(15)計(jì)算了嫦娥三號(hào)探月雷達(dá)預(yù)處理數(shù)據(jù)以及利用帶通濾波、EMD 方法和自適應(yīng)形態(tài)學(xué)方法處理結(jié)果的包絡(luò)熵，結(jié)果分別為6.041 0、6.036 7、6.014 5、5.873 3，可以看出，利用自適應(yīng)形態(tài)學(xué)濾波方法獲得的包絡(luò)熵結(jié)果相比帶通濾波和EMD 方法獲得的包絡(luò)熵結(jié)果更小，即剩余噪聲越少，這進(jìn)一步說(shuō)明了自適應(yīng)形態(tài)學(xué)濾波具備更強(qiáng)的噪聲壓制能力。

圖4 圖2（a）數(shù)據(jù)經(jīng)帶通濾波處理后的結(jié)果

圖5 圖2（a）數(shù)據(jù)經(jīng)EMD方法處理后的結(jié)果

圖6 圖2（a）數(shù)據(jù)經(jīng)自適應(yīng)形態(tài)學(xué)方法處理后的結(jié)果

4 結(jié)束語(yǔ)

本文提出了一種基于自適應(yīng)形態(tài)學(xué)的探月雷達(dá)噪聲壓制方法，該方法利用兩個(gè)不同尺寸的拋物線形結(jié)構(gòu)元素，構(gòu)造出能夠提取信號(hào)中間尺度范圍的形態(tài)學(xué)濾波器，并將包絡(luò)熵作為目標(biāo)函數(shù)，通過(guò)粒子群隨機(jī)搜索全局優(yōu)化方法對(duì)形態(tài)學(xué)結(jié)構(gòu)元素的幅值和長(zhǎng)度進(jìn)行優(yōu)選。通過(guò)對(duì)嫦娥三號(hào)探月雷達(dá)第二通道數(shù)據(jù)測(cè)試結(jié)果表明，相比常用帶通濾波方法、EMD 方法，該方法具有更強(qiáng)的噪聲壓制能力，能夠更有效地去除探月雷達(dá)信號(hào)中的干擾信息，大幅度提高了雷達(dá)信號(hào)成像的分辨率，從而為月表結(jié)構(gòu)的合理劃分提供了較好的支撐依據(jù)。同時(shí)，本文提出的方法可為嫦娥四號(hào)、嫦娥五號(hào)探月雷達(dá)數(shù)據(jù)的處理和解析提供技術(shù)方法參考。