利用shearlet變換對探月雷達數(shù)據(jù)進行月壤結(jié)構(gòu)的重建

張領(lǐng), 劉財, 胡斌* , 張澤強, 徐晨陽,曾昭發(fā), 王德利, 徐懿

1 吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026 2 澳門科技大學(xué) 月球與行星科學(xué)國家重點實驗室， 澳門 999078 3 中水東北勘測設(shè)計研究院有限責(zé)任公司， 長春 130061

0 引言

作為中國月球和深空探測的重要組成部分，嫦娥三號于2013年12月著陸月球雨海地區(qū)，對著陸區(qū)域進行了深入的探測(Xiao et al., 2015).嫦娥三號所攜帶的探月雷達通過發(fā)射高頻電磁波脈沖對地下結(jié)構(gòu)進行探測.為了實現(xiàn)探測目標，探月雷達設(shè)計了兩個通道的天線：第一個通道(CH-1)用于探測淺月殼結(jié)構(gòu)(中心頻率為60 MHz，帶寬為40～80 MHz)，第二個通道(CH-2)用于探測月壤的厚度與結(jié)構(gòu)(中心頻率是500 MHz，帶寬為250～750 MHz)(Fang et al., 2014).

國家天文臺最初發(fā)表了有關(guān)探月雷達的原始與處理結(jié)果(Su et al., 2014).Xiao等(2015)根據(jù)探月雷達的CH-1和CH-2數(shù)據(jù)對嫦娥三號著陸區(qū)域進行了地下地質(zhì)結(jié)構(gòu)解釋，獲得了深度360 m以內(nèi)的9層地下結(jié)構(gòu).Zhang等(2015)基于探月雷達數(shù)據(jù)對玄武巖漿的幾次覆蓋進行了解釋.Yuan等(2017)利用雷達數(shù)據(jù)對月球的平靜期進行了研究.值得注意的是，雷達結(jié)果解釋的前提是數(shù)據(jù)的正確性，然而，Li等(2018)仔細分析了CH-1數(shù)據(jù)并與在地面上采集的實驗數(shù)據(jù)進行對比，認為CH-1數(shù)據(jù)中深部的反射同相軸可能并非是來自地下的信息而是由儀器本身所引起的.

根據(jù)CH-2的雷達數(shù)據(jù)，嫦娥三號著陸區(qū)基巖之上的月壤層主要分為三層：最上層是厚度小于1 m的風(fēng)化層，中間是厚度小于5 m的濺射物層，最下面是厚度約為5 m的古月壤層(Fa et al., 2015).Zhang等(2019)和Lai等(2016)推論出著陸點上包括月壤與古月壤的多層結(jié)構(gòu)及其接受風(fēng)化與撞擊坑濺射覆蓋的演化過程.Zhang等(2018a)和Dong等(2017)根據(jù)CH-2數(shù)據(jù)計算了月壤的電性參數(shù)與鐵鈦含量.Hu等(2019)基于CH-2數(shù)據(jù)利用相似性約束來定位了月壤中的巖石塊.Ding等(2020)利用CH-2雷達數(shù)據(jù)與地形約束了著陸區(qū)的介電特性，表明紫薇撞擊坑濺射物形成的小山丘是由愛拉托遜紀月海玄武巖構(gòu)成.

探月雷達第二通道擁有兩個接收器天線，因此，可以得到兩組CH-2數(shù)據(jù)(CH-2A和CH-2B).因為CH-2A數(shù)據(jù)被一種橫向噪聲所干擾，所以前人的文章都只利用了CH-2B數(shù)據(jù).CH-2A數(shù)據(jù)中這種橫向噪聲表現(xiàn)為中深部的大量近水平多次波.這種噪聲可能是由于接收天線本身以及月球上復(fù)雜的采集環(huán)境所引起的.

實際上，不僅CH-2A數(shù)據(jù)被這種水平橫向噪聲所干擾，CH-2B數(shù)據(jù)中也存在少量這類噪聲.這類水平橫向噪聲嚴重影響了CH-2數(shù)據(jù)的深度利用與解釋.為了解決水平噪聲問題，我們引入了多尺度幾何分析方法中的shearlet變換對噪聲進行分析和抑制.與常規(guī)傅里葉分析和小波分析相比，多尺度幾何分析(例如contourlet變換(Do and Vetterli, 2003)，curvelet變換(Candès et al., 2006)，shearlet變換(Easley et al., 2008)等)可以更好地處理高維信號.它從局部性、方向性和多尺度分析信號，以找到具有表面奇異性或線奇異性的高維信號的最佳或最稀疏表示.近年來，許多多尺度幾何分析方法已應(yīng)用于噪聲抑制中(Zhang et al., 2010).這些多尺度和多方向分析方法利用噪聲和信號在方向和尺度上的分布差異來抑制噪聲(Liu et al., 2014).在這些方法中，shearlet變換具有更好的方向靈敏度和稀疏表示(H?user and Steidl, 2014).

本文將shearlet變換應(yīng)用于嫦娥三號探月雷達CH-2A和CH-2B數(shù)據(jù)的噪聲分析和抑制中，并融合兩者的結(jié)果對嫦娥三號著陸區(qū)的月壤結(jié)構(gòu)進行重建.首先，設(shè)計了一套針對第二通道雷達數(shù)據(jù)的預(yù)處理流程.隨后，利用月壤綜合模型驗證shearlet變換抑制水平橫向噪聲的效果.其次，利用shearlet變換對CH-2A和CH-2B數(shù)據(jù)進行處理，并將兩套處理后的圖像進行融合.最后，根據(jù)探月雷達處理結(jié)果對嫦娥三號著陸點的月壤結(jié)構(gòu)進行重建，重建結(jié)構(gòu)中的濺射物厚度與撞擊坑濺射物經(jīng)驗公式結(jié)果得到了相互印證.

1 數(shù)據(jù)與方法原理

1.1 探月雷達數(shù)據(jù)預(yù)處理流程

玉兔號月球車攜帶著探月雷達對地下地質(zhì)結(jié)構(gòu)進行探測，總共在月表上行走了114.8 m的距離(如圖1所示).圖1中的底圖由嫦娥三號著陸器搭載的降落相機所拍攝.在對探月雷達數(shù)據(jù)進行進一步的分析與解釋之前，需要對雷達數(shù)據(jù)進行預(yù)處理.

圖1 玉兔車行進路徑圖Fig.1 The Yutu′s path on the Moon

根據(jù)數(shù)據(jù)的質(zhì)量、采集環(huán)境以及采集參數(shù)設(shè)置，本文設(shè)計了針對嫦娥三號探月雷達數(shù)據(jù)的預(yù)處理流程.經(jīng)過預(yù)處理后，得到了兩套高分辨率的第二通道數(shù)據(jù)(圖2).

圖2 探月雷達第二通道數(shù)據(jù)預(yù)處理流程及結(jié)果Fig.2 The pipeline chart of CH-2 LPR data processing. N103—N209 denote the positions

N103—N209表示航點對應(yīng)位置.盡管針對CH-2數(shù)據(jù)進行了一系列的常規(guī)方法處理(例如，帶通濾波與去背景等)，CH-2A與CH-2B仍然被一種橫向噪聲所干擾(圖2中黃框部分).這種明顯的近水平的橫向多次波嚴重影響了中深部數(shù)據(jù)的解釋.這種橫向噪聲很可能來源于接受天線儀器本身和周圍采集環(huán)境的共同干擾.為了解決這些橫向噪聲的干擾，本文引入了多尺度幾何分析方法中的shearlet變換來分析與抑制噪聲.

1.2 Shearlet變換

Shearlet變換作為一種有效的多尺度與多方向的方法，可以基于多分辨率分析準確地生成shearlet系數(shù)(Yang et al., 2014).與使用仿射系統(tǒng)理論的方法相比，shearlet變換具有更好的稀疏性和更簡單的數(shù)學(xué)結(jié)構(gòu)(Guo et al., 2004, 2006; Heil, 2006).并且，即使在各向異性條件下，shearlet系統(tǒng)也能夠區(qū)分奇異分量并表征多維信號(H?user et al., 2014).shearlet變換的連續(xù)仿射系統(tǒng)可以描述如下：

{TtDMψ,M∈G,t∈Rn}，

(1)

(2)

膨脹矩陣Mas是兩個操作的結(jié)合Mas=SsAa，其中S代表剪切操作，A表示尺度操作：

(3)

其中s和a分別代表剪切因子和尺度因子.

此外，對于ξ=(ξ1,ξ2)∈R2,ξ1≠0，ψ定義為

(4)

(5)

t∈R2}.

(6)

考慮到公式(4)滿足連續(xù)shearlet變換的要求，定義f的連續(xù)shearlet變換如下：

SH(f)(a,s,m)∶=〈f·φa,s,m〉

(7)

圖3表示出了頻域中的shearlet基函數(shù)的劃分.當在尺度a，角度s和位置t上基函數(shù)有所不同時，Aa將會把數(shù)據(jù)分解為具有不同頻帶的若干子集.大尺度shearlet結(jié)果描述數(shù)據(jù)細節(jié)信息，小尺度shearlet結(jié)果表示趨勢.同理，Ss按方向劃分數(shù)據(jù)，s用來控制傾角范圍.

圖3 (a) 不同參數(shù)的水平shearlet基函數(shù)分頻示意圖； (b) 不同參數(shù)的垂直shearlet基函數(shù)分頻示意圖Fig.3 The frequency division of the horizontal shearlets (a) and vertical shearlets (b) for different parameters

2 模型驗證

為了驗證shearlet變換去除橫向噪聲的有效性，通過時間域有限差分方法(FDTD)模擬探月雷達合成數(shù)據(jù)(Zhang et al., 2018b).月壤模型和雷達模擬結(jié)果如圖4所示，雷達模擬的有關(guān)參數(shù)如表1所示.將背景噪聲和低切濾波的高斯隨機噪聲添加到探月雷達數(shù)據(jù)中(圖5)，形成含噪的合成數(shù)據(jù).本文的目標是衰減背景噪聲(低頻，水平方向)和隨機噪聲(中高頻)，因此設(shè)計了包括直流分量消除、去背景和帶通濾波在內(nèi)的常規(guī)處理方法組合，與shearlet變換去噪效果進行對比.

圖4 (a)月壤模型； (b) 雷達模擬結(jié)果Fig.4 (a) The lunar regolith model； (b) The simulation result.

圖5 (a) 噪聲； (b) 加噪數(shù)據(jù)Fig.5 (a) Noise； (b) Noisy data

表1 正演模擬參數(shù)Table 1 Parameter of simulation

圖6展示了常規(guī)處理手段方法和shearlet變換方法的對比結(jié)果圖.整體上，兩者對大部分噪聲都有較好的抑制作用.但是在常規(guī)處理結(jié)果中仍存在一些殘留噪聲，這是因為有用信號和噪聲信號之間存在頻帶的重疊.而shearlet變換能將數(shù)據(jù)映射到幾個子域中，這些子域可以稀疏地表示數(shù)據(jù)，使得噪聲和信號更好的分離，從而得到更好的處理效果.圖7通過頻譜可以看出，shearlet變換能夠更加真實地還原原始信號.

圖6 (a) 常規(guī)處理結(jié)果； (b) shearlet變換處理結(jié)果Fig.6 Denoised data, (a) the conventional method and (b) the shearlet method

圖7 頻譜分析結(jié)果(a) 原始數(shù)據(jù)； (b) 含噪數(shù)據(jù)； (c) 常規(guī)處理方法； (d) shearlet變換.Fig.7 Frequency spectrum of (a) original data, (b) noise, (c) conventional denoised data and (d) shearlet denoised data

本文所提出的shearlet變換方法不僅是一種抑制噪聲信號的方法，還是信號分析的有利工具.準確的數(shù)據(jù)分析，有針對性的策略和稀疏的表示方式，使得shearlet變換帶來更理想的去噪結(jié)果.

3 探月雷達數(shù)據(jù)分析、處理及融合

3.1 利用shearlet變換進行多尺度幾何分析

為更好地說明shearlet變換分析的過程，本節(jié)以CH-2A數(shù)據(jù)中N203到N204航點的數(shù)據(jù)為例(圖10a).嫦娥三號探月雷達數(shù)據(jù)，尤其是CH-2A數(shù)據(jù)，存在明顯的橫向水平噪聲(80～150 ns)，這嚴重影響深部信息的獲取.為了更準確地研究月壤深部結(jié)構(gòu)信息，使用shearlet變換(a=4,s=4)來分析探月雷達數(shù)據(jù)，并研究噪聲和信號在頻率和方向上的分布情況.Shearlet變換按頻率和方向?qū)⑿盘栍成涞讲煌淖蛹?將相同比例或角度的子集相加，可獲得具有不同頻帶和不同方向帶的子信號(見圖8和圖9).

圖8 CH-2A中N203—N204段數(shù)據(jù)不同尺度下的shearlet分量Fig.8 The shearlet coefficients of the selected portion (N203 to N204 in CH-2A data) in different scales

圖9 CH-2A中N203—N204段數(shù)據(jù)不同角度下的shearlet分量Fig.9 The shearlet coefficients of the selected portion (N203 to N204 in CH-2A data) in different angles

圖8顯示了雷達數(shù)據(jù)三個主要尺度(頻率)分量.可以看到信號的頻帶低于噪聲，信號集中在中低頻段，而噪聲集中在中高頻段，因此可以在低尺度上觀察到月壤的結(jié)構(gòu).圖9顯示了雷達數(shù)據(jù)三個主要角度分量.水平噪聲主要集中在水平角度分量中.總之，對于探月雷達的第二通道數(shù)據(jù)，抑制噪聲和提高降結(jié)構(gòu)識別能力的關(guān)鍵是消除數(shù)據(jù)信號中的水平分量.

3.2 利用改進的水平滑動濾波在shearlet域中去除噪聲

從圖9b中可以看到，在水平的分量中仍然有許多有用的信號(0～60 ns).本文認為水平shearlet分量中的水平橫向噪聲是背景場.因此，提出了一種基于shearlet變換的水平噪聲抑制方法.該方法實際上是在傳統(tǒng)的背景消除方法，即水平滑動濾波器的基礎(chǔ)上，代替直接去除水平方向上的小波系數(shù)，使用水平滑動濾波來盡可能多地保留有用信號(曾昭發(fā)等，2010).水平滑動濾波的目的是去除水平方向上的慢速振蕩分量，以突出顯示快速振蕩的分量.該方法的原理是確定窗口的水平寬度，對窗口中的所有道集求平均，并從窗口中間減去平均值，如公式(8)所示：

(8)

其中，x′(m,n)是處理后的結(jié)果，xc是窗口中間的數(shù)據(jù)，x(m,n)是原始數(shù)據(jù)，m和n是坐標，N是道集數(shù)，用于控制處理數(shù)據(jù)的平滑度.

然而，時空域中的水平滑動濾波仍然會破壞淺部的有用信號．因在shearlet域中，子集的大小與原始數(shù)據(jù)的大小相同，為了抑制shearlet域中水平子集分量的背景，shearlet域中的水平滑動濾波可以表示為

(9)

其中，c′s(m,n)是處理的shearlet系數(shù)，csc是窗口中間的shearlet系數(shù)，cs(m,n)是原始的shearlet系數(shù).

圖10顯示了處理后的結(jié)果.顯然，大多數(shù)的水平噪聲都得到了有效抑制，有效信號(90～130 ns)得以恢復(fù).根據(jù)處理結(jié)果，可以勾畫出一個新的分界面(圖10c中箭頭所示).所提出的方法會稍微破壞淺部的有用信號(0～60 ns).去噪后的數(shù)據(jù)為研究月壤的結(jié)構(gòu)提供了良好的圖像.

圖10 N203—N204段數(shù)據(jù)的噪聲抑制(a) 原始數(shù)據(jù)； (b) shearlet變換分離出的噪聲； (c) 去噪結(jié)果.Fig.10 Demonstration of horizontal noise attenuation for the selected portion (N203 to N204 in CH-2A data)

3.3 CH-2A與CH-2B處理結(jié)果融合

將shearlet變換去噪方法應(yīng)用至第二通道的全部數(shù)據(jù)，圖11a和圖11b分別顯示的是CH-2A與CH-2B數(shù)據(jù)shearlet變換去噪后的結(jié)果.經(jīng)過處理后，CH-2A大部分的橫向近水平噪聲都得到了有效抑制(紅框).特別是在10～30 m區(qū)域，有用信號(受噪聲干擾)在一定程度上得以恢復(fù).雖然CH-2A數(shù)據(jù)的處理結(jié)果中存在一些噪聲殘留，這是因為為了保留更多有用信號，只刪除了最水平的shearlet系數(shù).與CH-2A數(shù)據(jù)不同，CH-2B數(shù)據(jù)中的水平噪聲雖然幅度較小，但是它仍然會影響深部有用信號.CH-2B數(shù)據(jù)經(jīng)過處理后，深部數(shù)據(jù)信號更清晰(紅框)，這對于確定月壤的下邊界具有一定意義.

圖11 去噪后的CH-2數(shù)據(jù)(a) CH-2A； (b) CH-2B.Fig.11 Horizontal noise attenuation for LPR CH-2 data

盡管CH-2A數(shù)據(jù)受到噪聲困擾，但經(jīng)過處理后，嘈雜的CH-2A數(shù)據(jù)信噪比有了顯著的提高，其中淺部的信號(30～50 ns)比CH-2B數(shù)據(jù)清晰.這些淺部信號對于后面將要解釋到的撞擊坑濺射物的厚度至關(guān)重要，因此為了全面利用CH-2A和CH-2B兩個數(shù)據(jù)集的整體優(yōu)勢，本文借鑒了圖像融合的思想.具體步驟如下(圖12)：

(1)使用shearlet變換將探月雷達CH-2A和CH-2B數(shù)據(jù)映射到不同的稀疏域.

(2)將稀疏域系數(shù)分為低頻尺度和高頻尺度.

(3)對于低頻尺度的系數(shù)，選擇兩個數(shù)據(jù)的平均值作為融合系數(shù).

(4)對于高頻尺度的系數(shù)，選擇兩個數(shù)據(jù)中的最大值作為融合系數(shù).

(5)將融合系數(shù)映射回時空域，獲得最終融合結(jié)果.

圖12 CH-2A和CH-2B數(shù)據(jù)圖像融合流程Fig.12 Image of LPR CH-2A and B data fusion flow chart

經(jīng)過圖像融合以后，獲得一個深部信息與淺部信息更加均衡的結(jié)果(圖13).

圖13 CH-2A和CH-2B融合結(jié)果，時深轉(zhuǎn)換所用的相對介電常數(shù)為3Fig.13 The stacked result of LPR CH-2A and B data. The relative dielectric constant is set as 3

4 結(jié)果解釋

根據(jù)月壤的形成機制、嫦娥三號著陸區(qū)域特點以及前人對探月雷達的數(shù)據(jù)分析，對融合后的CH-2數(shù)據(jù)進行分層解釋(圖14)：

(1)基巖.基底巖石是在雨海紀晚期或埃拉托遜紀早期所形成的玄武巖覆蓋層(McKay et al., 1991).

(2)古月壤層.在埃拉托遜紀和哥白尼紀，各類隕石撞擊月球表面的玄武巖蓋層形成古月壤層(McKay et al., 1991；歐陽自遠，2005).圖中的藍點代表未徹底風(fēng)化的巖石碎塊或隕石碎塊(Hu et al., 2019；Fa, 2020).

(3)濺射層.在哥白尼紀的某個時間點，嫦娥三號著陸點附近的撞擊坑形成(圖14)，其所形成的濺射物覆蓋了其周圍的月壤層(Xiao et al., 2015).圖中的黃點代表濺射出的巖石碎塊(Hu et al., 2019；Fa, 2020).

圖14 嫦娥三號著陸區(qū)域月壤層結(jié)構(gòu)解釋重建圖Fig.14 The interpretation of regolith structure on the CE-3 landing site

另一方面，嫦娥三號降落在一個直徑約為450 m的撞擊坑的附近，距撞擊坑的邊緣約50 m處.通過眾多撞擊坑及其濺射物厚度的研究及公式(Mcgetchin et al., 1973)可以計算出著陸點的濺射層的厚度為4 m左右，與上述的分層結(jié)果(圖14)一致，也側(cè)面驗證了融合數(shù)據(jù)的有效性與分層結(jié)構(gòu)的可信度.

5 結(jié)論

本文利用shearlet變換對探月雷達數(shù)據(jù)進行了噪聲分析與去除，并對嫦娥三號著陸區(qū)的月壤結(jié)構(gòu)進行重建.首先，在shearlet域中，觀察橫向噪聲的分布，并分離出以噪聲為主的shearlet分量，得到重建后的雷達數(shù)據(jù).隨后，為充分利用嫦娥三號第二通道兩套數(shù)據(jù)的整體優(yōu)勢，對兩套經(jīng)過shearlet變換去噪后的雷達數(shù)據(jù)進行融合，得到綜合探月雷達圖像，有效地去除了水平噪聲并增強了來自淺月表的有效信號.最后，根據(jù)探月雷達處理結(jié)果對嫦娥三號著陸點的月壤結(jié)構(gòu)進行重建，重建結(jié)構(gòu)中的濺射物厚度與撞擊坑濺射物經(jīng)驗公式結(jié)果得到了相互印證.

在分析了CH-2A和CH-2B的融合數(shù)據(jù)之后，獲得了嫦娥三號著陸區(qū)域的月壤結(jié)果與演化模型.其月壤層可大致分為兩層：覆蓋在基質(zhì)玄武巖之上的古月壤層和由附近撞擊坑濺射物形成的濺射層.通過濺射物厚度經(jīng)驗公式，驗證了融合數(shù)據(jù)的有效性與分層結(jié)構(gòu)的可信度.

嫦娥四號探月雷達到目前為止也取得了階段性成果，包括利用第一通道對深部結(jié)構(gòu)進行解釋(Zhang et al., 2020; Lai et al., 2020)，以及利用第二通道對嫦娥四號著陸區(qū)域濺射層與月壤層的解釋(Li et al., 2020; Lai et al., 2019).其中嫦娥四號利用的第二通道也均是CH-2B數(shù)據(jù).本文方法為嫦娥四號第二通道的兩套數(shù)據(jù)(CH-2A和CH-2B)的綜合利用提供了依據(jù)，同時為后續(xù)的嫦娥五號以及天問一號的雷達數(shù)據(jù)處理提供了處理方法支撐.