李陽平 黃 玲 王 珂 趙海峰*
①(中國科學院大學 北京 100086)
②(中國科學院空間應用工程與技術中心中國科學院太空應用重點實驗室 北京 100094)
③(中國科學院空天信息創新研究院 北京 100094)
月球是與地球最近的天體,但時至今日,人類對月球表面月壤結構的勘探能力和成果仍十分有限。在美國Apollo任務中[1],宇航員僅在近地表提取了月壤與月巖樣本,月球更深部的原狀結構信息和深部樣品無法獲取。我國嫦娥探月工程項目在月面上采用了探地雷達系統進行探測[2–4],科學研究結果表明探地雷達系統不僅能夠對月球地下介質分布開展有效的探測,也能為月面深部鉆采的場地優選和路徑規劃提供依據,增強月面鉆探的能力,降低鉆探取樣的風險。由于探地雷達的高頻電磁波在月壤及月巖砂礫介質中具有很強的散射特征和各向異性特征,探測回波復雜,并且數據有限,采用傳統的數據處理與目標成像技術難以對月壤中的月巖體形態進行精確成像,無法高效地適用于未來月球探測任務中。本文針對月面淺層(5~10 m深)鉆進采樣的科學需求,在雷達探測數據處理引入機器學習算法,增強探月雷達探測數據實時地層成像的效率,最終為月球探測任務的月巖鉆進避障和鉆進策略優化提供一種可行的解決方案。
由于近年來計算機計算性能的提高,機器學習算法的爆發性發展,人工智能神經網絡算法相對于傳統的模型反演算法展現出了巨大的發展潛力和優勢。研究成果表明人工智能神經網絡在探地雷達信號處理中具有較高效準確的優勢[5–7]。深度神經網絡由輸入層、隱藏層和輸出層組成,輸入數據的特征變量作為多變量線性組合被傳遞給隱藏層的節點,隱藏層的每個節點設置非線性激活函數,通過每層的線性組合傳遞至輸出層,從而實現神經網絡對復雜非線性關系的建模。
隨著深度神經模型的不斷優化和擴展,已有學者將深度學習應用于探地雷達回波信號分類。Gao等人[8]利用基于區域的快速卷積網絡方法(faster R-ConvNet)進行路面缺陷檢測,探測精度可達89.13%。Li等人[9]利用快速目標檢測(You Only Look Once, YOLO)模型檢測瀝青路面的隱藏裂隙,并在小數據集上也取得良好效果。2020年,Khudoyaro等人[10]將深度3維卷積神經網絡應用于3維探地雷達數據,識別地下物體,將城市道路分類為管道、孔洞、檢修孔和土壤,其準確率達到了97%。然而,當前的研究方法需要花費大量人力進行特征標定,而對于探地雷達圖像識別而言,目標的精確幾何特征仍無法定量描述,因此只能將其輸出特征分類為不同的地質結構,分類結果單一且粗糙。
為了獲得更準確的掩埋物體特征,Liu等人[11]提出基于全波形反演的檢測算法,能夠準確恢復地下掩埋的圓柱體位置和直徑,但是麥克斯韋方程組的計算需要消耗大量資源,無法實時獲得結果。于是,Giannakis等人[12]將深度學習與探地雷達(Ground-Penetrating Radar, GPR)信號處理結合,實時預測掩埋物體的深度,但目前僅支持對規則圓形進行處理,無法針對復雜表面預測,地質環境也局限于混凝土與鋼筋結構,難以處理月球環境。
因此,本文結合深度神經網絡在非線性問題中的優勢,搭建了一種能分析預測不規則介質上表面幾何特征與位置的深度學習模型。首先利用邊緣檢測等數字圖像處理方法獲取月巖輪廓,構建月下地質結構環境,模擬仿真得到探地雷達回波信號數據集;其次,利用主成分分析法提取雷達回波數據的主要特征,提取后的特征能夠加強信號的表現力,作為神經網絡模型的輸入數據;最后,根據反向傳播算法,搭建了一種深度神經網絡模型,在訓練過程中自主調整神經元節點權重,實現對掩埋介質的上表面幾何特征與位置的回歸預測。
與傳統探地雷達數據處理方法相比,本文所提模型有以下優點:(1)本文針對月球環境,采用數字圖像處理方法獲取掩埋介質輪廓,并構建大量接近真實幾何特征的月球地質模型;(2)本文采用的深度神經網絡模型,其中隱藏層神經元節點的激活函數具有非線性特征,從而對非線性數據處理具有較好的預測結果;(3)本文對于雷達信號進行分類,利用訓練完成后的數字化成像模型,使用時直接輸入回波信號數據獲得成像結果,計算消耗小,耗時低,能夠滿足未來航天任務設計的實時性要求。
在探地雷達的人工智能神經網絡訓練中,訓練數據通常來自兩種情況:(1)相對有限的一組探地雷達真實數據;(2)基于簡單規則形狀模型而獲取的仿真模擬數據。這些數據有如下缺點:(1)真實數據的數據量不足、環境或系統影響因素多;(2)簡單模型獲得的模擬仿真數據不能反映實際情況。在當前更快速的全波形正向求解器和改進的分析方法背景下,如何獲得更逼真的月壤模型模擬數據是解決問題的關鍵。本文根據月面環境,基于電磁波在地下的傳播特性,結合嫦娥3號和4號搭載的月球探測器參數,以及Apollo系列任務返回的月壤與月巖樣本圖像,構建了一種接近真實的月面地質環境的建模和回波數據分類方法。
探月雷達的工作原理是發射天線將超寬帶的毫微秒電磁波脈沖向月面輻射,當遇到電性特征不同的物質(如熔巖層和濺射物)時,電磁波將會發生反射和散射現象。探月雷達的接收天線將接收到的信號進行放大、采樣,隨后對獲得的探測信號進行處理、反演成像,最終可獲得月球地下介質的結構分布。
為深入探究不同深度下復雜月球地質結構特征與探月雷達回波信號的關系,本文首先針對月球淺表層0.5 m深單個掩埋非規則形狀月巖塊體進行仿真模擬,取得預測結果后,再進行10 m深月球復雜地質環境模擬。其月面雷達探測系統的主要技術指標如表1所示。

表1 探月雷達的參數設計
針對月面下0.5 m的淺層月面雷達穿透深度,設計月面雷達的中心頻率為1.2 GHz。在月球環境條件下,其探測距離向分辨率約為2 cm,該月面雷達系統可用于探測月面表層掩埋月巖的輪廓和分布位置。針對10 m深的探測目標,設計月面雷達中心頻率為700 MHz。在月壤條件下,其探測距離向分辨率約為20 cm,該月面雷達系統可用于探測月壤厚度及月巖分層結構。
在正演仿真模型中的月壤中月巖介質物性參數和幾何參數設計,則采用基于Apollo 14任務返回的巖石樣本作為參考模型[13],從而可使得仿真模型和月巖表面形態更接近真實的月面條件與環境,能更有效地檢驗方法的有效性。其巖石樣本部分照片如圖1所示。

圖1 Apollo 14號返回樣品圖像
構建幾何參數的過程中,首先通過提取真實月巖照片的月巖輪廓來構建2維gprMax仿真模型[14],從而可精確地模擬月巖樣品表面較為真實的粗糙形態特征。在2維月巖圖像中,圖像存儲的大部分信息主要由邊緣組成,表現為局部特征不連續,可將其理解為圖像灰度變化劇烈的區域,因此利用梯度算子來進行邊緣檢測是效果較好的一種方式。利用Canny梯度算子[15]對圖像進行邊緣檢測后,獲得結果如圖2(a)所示。
可以看出,圖像邊緣的噪聲較多,中間存在許多的空穴,其與真實月巖表面形態仍存在較大的差別,因此采用膨脹腐蝕算法進行處理,然后進行孔洞填充,最后得到與真實月巖表面形態更為接近的月巖截面輪廓,如圖2(b)所示。

圖2 獲取月巖輪廓
在獲取月巖截面圖像后,利用gprMax軟件進行2維有限差分數值模擬。gprMax軟件是一個模擬電磁波傳播的開源軟件[14]。它基于Yee算法,通過有限時域差分方法求解麥克斯韋方程組。在模型的構建中,還能導入便攜式網絡圖形(Protabel Network Graphics, PNG)格式的圖片,通過圖片中不同顏色像素來定義其材料特征,由此可以構建出月巖和月壤二相介質的模型。
在使用gprMax進行正演模擬中,首先需要確定仿真模型區域,并劃分網格大小。本文所使用的環境設定為:CPU Intel Xeon Gold 5118, 256 GB內存,Tesla V100 16 GB顯卡的服務器,由于數據量和計算時間的限制,本文將淺層2維月球地質模型區域設計在0.7 m×0.5 m的空間中,網格大小為0.005 m,然后將巖石位置隨機分布在空間中的任意坐標。其中一個代表性模型如圖3所示。

圖3 模擬二相不規則月巖環境
在此模型中(圖3),月巖體底部距月壤表面的埋深為0.25 m。波形的中心頻率設置為1.2 GHz。而時間窗口設置為7.5 ns。雷達的發射天線與接收天線的間距為2 cm,測線上的測點間距為1 cm,即雷達自模型最左側開始行進,每隔1 cm進行電磁波的發射與接收,從而獲得雷達探測的B-Scan數據。
為進一步研究含有多月巖塊體地層的回波信號特征與數據處理方法,本文將2維月球地質模型區域設置為10 m×20 m。在模型中隨機掩埋不規則形態和大小的月巖介質,生成不同巖石含量的地質環境,其中兩個代表性模型如圖4所示。

圖4 多目標非均勻月巖地層
在此模擬環境中,模擬區域設置為20 m×10 m,網格大小為0.05 m,巖石大小分為3種:10 cm左右、30 cm左右和50 cm左右??刂拼笾行∈^比例,均勻地將它們散布在月壤地層中。而當月巖個數逐漸增多時,需要控制巖石面積百分含量的范圍。本文按照黑色像素點/白色像素點的比值,來計算巖石含量。
構建月球地質環境的過程中,需要考慮的月巖電性特征主要包括介電常數、電導率、磁導率、磁損耗。Apollo系列任務返回的月壤與月巖樣本均進行了介電常數、密度的測量。采用的月壤參數如表2,月壤參數基于450MHz雷達下測量[16]。

表2 月壤電性參數表
不規則形態月巖的電性參數設置如表3所示。

表3 月巖電性參數表
研究表明,介電常數的實數部分大約與頻率無關,在1 MHz下的探測結果與450 MHz的結果非常接近。而損耗正切角通常是頻率的復雜函數,不同頻率下的結果各不相同[16]。因此,本實驗中,只模擬了月巖與月壤的介電常數,而并未模擬損耗正切值。
本文基于模擬月巖的幾何與形態參數進行gprMax仿真模擬。在獲得大量雷達信號數據集的前提下,提取數據中的主要特征,作為神經網絡的輸入數據。最后采用反向傳播算法(Back Propagation, BP),構建深度神經網絡,對地質結構特征進行預測,繪制表面輪廓圖像。
神經網絡模型的輸入數據為探地雷達回波數據,與其對應的輸出數據為此收發射天線垂直界面巖石的上表面各像素位置。然而,0.5 m深雷達單道數據的回波量為339維,如果全部作為深度學習的輸入數據,則數據量過大,導致模型泛化性能差。因此需要提取雷達回波數據的主要特征。主成分分析(Principal Component Analysis, PCA)是機器學習研究中的一種常用的降維方法[17]。主成分分析法提取特征時構建一個n×m的變換矩陣,該矩陣能夠將樣本x(即339維的單道回波數據)映射到一個新的特征子空間中,并使得變換后的特征子空間維數m少于原始特征空間維數n,其主要思路為

在本文的研究中,解釋的特征值比率累計和在95%以內,即盡可能降低信息損失。通過分析對比多參數、多模型條件下的訓練結果,選擇了最合適的特征維數,將0.5 m地質模型特征的339維回波數據降為39維主特征。
在完成雷達觀測數據的降維壓縮之后,本文利用反向傳播算法的思想,搭建神經網絡模型。模型如圖5所示。

圖5 神經網絡模型示意圖
輸入層有N個訓練樣本,每個訓練樣本中包含回波數據經過主成分分析后的39維主特征,對于訓練樣本集X=[x1,x2,···,xk,···,xn]來說,預測輸出為Y=[y1,y2,···,yk,···,yn],目標輸出則為T=[t1,t2,···,tk,···,tn]。初始權重ω的選取則一般采用統計學上隨機抽取的思想,在算法中隨機產生大量的初始值,逐次進行迭代優化。在BP算法中,由于函數大多輸出值在[–1,1],因而初始權重選擇一般也在這個范圍內。采用逐步搜索法來解決初始權值選擇的隨機性。先將初始權值區域N等分,在N個區域內再重復上述步驟,直到誤差函數不再減小,認為找到了最優解。只要區域選擇得足夠小,這種方法能夠有效避免局部極小問題。算法核心就是將復雜的雷達回波信號識別轉化為回歸問題,從而預測不規則介質的上表面幾何特征與位置。
其模型的主要參數如表4所示。

表4 神經網絡模型結構參數
數據輸入輸入層后,需要衡量整個模型是否找到了滿足誤差最小的最優解,總誤差采用MSE函數,其表達式為

在本文中,經過模擬實驗,隱藏層選擇過少會導致模型擬合度降低,而隱藏層過多時模型泛化能力較差,因此隱藏層被確定為8層。其中,第1層與第2層的激活函數為ReLU[18],神經元共64個。第3層和第4層激活函數為Sigmoid[19],神經元共64個。接下來第5層~第8層的激活函數都為ReLU,而第5層與第6層的神經元個數為32個,第7與第8層的神經元個數為16個。其結構如表5所示。

表5 神經網絡結構
根據探地雷達電磁波的工作方式,在前期實驗中,為了研究神經網絡算法對月巖幾何形狀預測的效果,因此首先將不規則月巖形狀掩埋物體的介電常數設置為與月壤差異較大的數值,以便逐步深入研究方法的可行性和效果。利用gprMax對圖3所示的模型結構進行仿真,掩埋的月巖上表面深度為0.11,由此得到一道A-Scan,與多道A-Scan組成的截面圖像B-Scan。其結果如圖6所示。

圖6 不規則月巖雷達回波信號
數據樣本的豐富度影響模型是否能預測此條件下的任意月巖形態,本文首先隨機生成了8000個模擬月巖環境,由于數據量不足,因此訓練后的神經網絡模型預測精度非常低。增加數據量至26000個模擬月巖地質環境。圖7是不同數據量對訓練預測結果的影響。

圖7 不同數據量對訓練結果影響
獲取足夠的數據量后,利用人工神經網絡開始進行學習預測,模型的評價標準則由R-square確定系數[20]來確定。其表達式為

其中,y1, y2,···,yn表示測量值,f1, f2,···,fn表示預測值。R-square取值在0~1,越接近1時則表明回歸模型表現效果越好。由此確定系數,得到較大差異介電常數不規則形狀掩埋月巖體上表面的位置預測結果,R-square為0.93,如圖8所示。

圖8 較大差異介電常數不規則形狀掩埋月巖體上表面輪廓的預測對比圖
由此可見,此方法預測的上表面較為清晰,誤差較小,誤差范圍一般在10 cm以內。最后,從預測結果中隨機選擇3個地質環境,對其上表面輪廓進行繪制,其結果如圖9所示。

圖9 較大差異介電常數不規則形狀掩埋月巖體上表面輪廓預測結果的誤差直方圖
圖10(a)的上表面形態預測準確率為93.76%,圖10(b)的上表面預測準確率為96.04%,而圖10(c)的上表面預測準確率是95.13%左右。

圖10 較大差異介電常數不規則形狀掩埋月巖體上表面輪廓預測示例
在完成較大差異介電常數不規則形狀掩埋月巖體模型的仿真與成像研究,并取得良好效果之后,進一步將模型月巖介電常數介質材料變換為表3所示(Apollo號返回的月巖樣本材料數據資料)。模擬仿真如圖11所示。
由圖11可知,月巖樣本材料掩埋介質的上表面位置預測的R-square為0.78,與大介電常數掩埋介質的預測準確率有一定差別。

圖11 月巖材料不規則月巖上表面輪廓預測
在獲取單個巖石的上表面輪廓預測結果后,針對含有多月巖塊體地層的模型(對應圖4(a)與圖4(b)),其回波信號B-Scan圖像如圖12所示。
由于掩埋月巖在深度方向上互相交疊,出現反射、散射、多次回波等現象,因此回波信號十分雜亂。且回波數據特征空間維度較高,利用PCA提取主要特征后,降維的輸入數據為56維。巖石含量為38%的預測結果如圖12(b)所示,深度學習模型的預測準確率R-square值僅為0.48。

圖12 含有多月巖塊體地層模型雷達回波信號的B-Scan多道雷達回波數據圖像
對月巖上表面輪廓進行模型和預測結果對比分析,可見0.5 m深的大介電常數模型可以看出大部分區域擬合得較好,誤差大約在10 cm以內,但下表面輪廓由于預測精度較低,因此無法進行擬合。其原因是: (1)受上表面強反射影響,下表面的反射信號回波幅值弱;(2)部分下表面的回波信號與上表面的2次回波信號重疊,難以區分上下表面回波信號,這也是探地雷達數據處理的共性問題。采用本文所提方法,其優勢在于,在將模型訓練好之后,只需要約為0.0032 s的時間即可獲得預測結果,極大提高了數據處理效率,達到現場處理與解釋的能力。分析結果還表明當月巖的介電常數與月壤接近時,預測精度有所下降,但仍可預測其表面輪廓的形態趨勢。
如圖13所示,在多月巖塊體地層模型的仿真與模型形態預測中,雷達仿真數據復雜程度增大,回波數據包含多目標散射、多次回波信息。月巖表面的模型形態與預測結果形態之間的匹配關系下降較為嚴重,但其仍顯示出一定的線性相關性,說明預測方法仍能反映出月巖表面形態的主要特征。為進一步提高預測精度,需要在后續研究中首先從大量回波數據中提取反射信息,壓制目標間的多次反射與回波,之后再進行主成分分析,提取信號主要特征。另外,有限訓練數據集無法滿足復雜模型的回歸擬合優化需求,需要增大多月巖復雜模型的仿真訓練數據多樣性和數據量,以使得模型泛化能力強,增強預測的穩定性。

圖13 含有多月巖塊體地層模型上表面輪廓預測結果
本文基于深度學習方法的優勢,針對月壤下月巖形態探測實時處理的需求,將探地雷達回波信號的解釋和幾何重建問題轉化成深度學習中的回歸問題,并通過仿真分析來驗證該思路的可行性。
(1)本文通過提取Apollo 14號返回的月巖的2維輪廓,構建月巖形狀的2維模型,并根據Apollo返回樣品的電性參數來設置模型的電學特性,從而更真實地模擬了月巖結構的真實形態。
(2)利用gprMax軟件構建探月雷達觀測數據集,針對獲得的高維回波數據,采用主成分分析方法來降低模型輸入數據的復雜程度,防止過擬合現象。然后,構建神經網絡模型對月巖幾何形態進行直接預測。
(3)在基于2維模擬數據建立的神經網絡模型預測中,當月壤中掩埋的月巖具有較大介電常數材料時,上表面的預測準確率達到0.93,能較好地重構出上表面輪廓,與Giannakis等人[12]的混凝土鋼筋規則圓形介質預測實驗相比,本文工作已能夠對不規則介質的上表面進行有效的預測重構;利用已訓練完成的模型進行數據處理,預測時間約為0.003 s,能夠達到實時性預測的要求。當月巖樣本介電常數材料接近月壤時,上表面的預測準確率為0.78,表明月壤和月巖介電常數對比度較小時,模型預測精度下降。
(4)在含復雜多月巖塊體分布的地層模型上表面輪廓預測中,由于回波信號復雜,深度學習模型的預測準確率僅為0.48,其原因可能是復雜月巖分布導致雷達信號散射疊加,在后續工作中,可以研究不同掩埋介質深度、不同巖石含量下的模型擬合效果,并對比其他神經網絡方法,選擇更優化的數學模型,提高預測準確率。
在后續研究中,需要對更多地質介質的物性參數數據庫進行統計和分類,建立完整的基于數據驅動的神經網絡地層結構成像方法的數據集和理論體系。此方法對地球或者地外地質探測工作有著非常大的應用價值和科學意義。