楊望笑 竇銀科* 稂時楠 趙 博 王煜塵 左廣宇

①(太原理工大學電氣與動力工程學院 太原 030024)

②(北京工業(yè)大學信息與通信工程學院 北京 100124)

③(中國科學院空天信息創(chuàng)新研究院 北京 100190)

1 引言

調頻連續(xù)波(Frequency Modulated Continuous Wave, FMCW)雷達技術目前已在測距、氣象監(jiān)測、極地冰雪探測等領域得到廣泛應用。南極冰蓋密度對于估算冰蓋表面物質平衡是至關重要的[1,2]，文獻[3]闡述了冰蓋密度與介電常數之間的經驗公式，因此，冰蓋密度可經雷達反演冰蓋垂直剖面上的介電常數轉換得到。

目前通過雷達對冰蓋密度/介電常數進行反演的主要方法是基于共中心點的雷達反射信號。發(fā)射和接收天線之間的距離需在0～80 m范圍內不斷變化，提取某一反射層在不同天線間距下的傳播時間差，再通過Dix反演、相似性分析、干涉測量等方法估算冰蓋表面至反射層之間的等效介電常數，經文獻[3]中的經驗公式轉換為平均密度[4,5]。由于這種方法數據采集過程較為繁瑣，在進行連續(xù)的介電常數剖面反演時使用較少。

當采用單條雷達反射信號進行介電常數反演時，車輛或飛機在采集數據過程中可同時搭載發(fā)射和接收天線，高效獲取冰蓋內部結構和介電常數剖面，這對大范圍了解冰蓋物質平衡是至關重要的。但是由于少量液態(tài)水會導致雷達信號強烈的反射和衰減，因此基于單條反射信號的反演法無法在濕雪帶進行介電常數剖面反演，更適用于高海拔、低積聚率地區(qū)，且在反演前需分析研究區(qū)域中冰的介電特性以確定該類方法的可行性[6]。目前使用單條雷達反射信號反演介電常數的主要方法是高斯牛頓法、退火參數算法或粒子群算法尋找介電常數的最優(yōu)解[7,8]，但是這些方法的計算過程相對復雜，更適用于層位尺寸較大、數量少的應用條件。實際上冰蓋是由無數個薄層組成，每個層之間都會存在微小的介電常數差異，因此并不適用于傳統(tǒng)尋找最優(yōu)解的方法。剝層反演法可提取反射信號峰值的幅度與時延，當已知第1個層位的介電常數時，可根據幅度估算下一層的介電常數，進而遞推求出整個介電常數剖面。剝層反演法在應用過程中還存在方法誤差，導致反演結果在局部深度范圍內存在誤差[6]。因此有必要對基于單條反射信號的介電剖面反演算法進行改進。

為了彌補上述算法的不足，本文提出一種基于剝層反演法和密實化公式相結合的優(yōu)化算法。采用快速傅里葉變換(Fast Fourier Transform, FFT)對差頻信號進行頻譜估計后，首先使用剝層反演法對回波信號進行初步的密度反演，再根據冰蓋密實化機理進行校準，基于校準結果仿真差頻信號，通過與原始回波信號進行對比尋找校準方程的最優(yōu)解。對理論冰蓋模型和實測雷達數據的實驗結果表明，此算法可有效提高反演精度。

2 FMCW冰雷達工作模型

冰蓋是空氣與冰的混合物，可將其內部結構假設為L個平行且有介電特性差異的冰層。FMCW雷達差頻信號的時域表達式為

以6.25 MS/s的采樣率對差頻時域信號采樣，獲取N=25000個采樣點。為降低頻譜估計誤差，在N個采樣點之后填充40×N個0，再進行FFT。回波信號由L–1個峰值組成，如圖1所示，其中每個峰值表示層位交界面產生的反射。

圖1 FMCW雷達工作示意圖

回波信號中第k個峰值的位置表示在深度dk處交界面反射信號的時延τk，該峰值幅度AMP的表達式為

3 反演算法

本文算法是基于剝層反演法，并在此基礎上進行校準的一種優(yōu)化算法，下面分別對算法的應用方法和優(yōu)化過程進行介紹。

3.1 剝層反演法在冰蓋上的應用

由式(3)—式(8)可知，每個峰值的幅度由層位間介電常數差異決定。反演時可以提取回波信號中的峰值，當已知某一個層的介電常數εk即可根據差頻信號幅度和式(3)、式(4)、式(8)估算出下一層的介電常數εk+1，遞推求出整個垂直剖面上的介電常數[6]。但是該方法更適用于層位厚度明顯大于雷達分辨率的情況[10]。因此，對傳統(tǒng)的剝層反演法進行改進十分必要。

1927年的農村調查，也折射了讀書人在村民眼中的尷尬處境：“讀書成本太大，出來非但沒有官做，即教員位置亦粥少僧多，而況學些空架子，文不象秀才，武不象丁，手不能提籃，肩不能挑擔，不事生產，要吃要用?！盵21]

文獻[11]分別公布了冰芯B31的密度和介電常數測量數據。圖2(a)黑線揭示了冰芯5 mm深度分辨率的相對介電常數數據，對30 m附近的數據放大后表明相鄰的層位也會表現(xiàn)出介電常數差異，這意味著差頻信號中峰值所表示的信號反射并不是由一個冰層的介電特性突變造成的，而是無數個冰層之間微小的介電常數變化引起的許多微弱反射相互干涉形成的[12]。除此之外，當電磁波向由低向高介電常數，或由高向低介電常數的冰層傳播時，都會產生電磁波的反射，因此，直接使用差頻信號峰值進行剝層反演會產生明顯的反演誤差。

綜上，在進行剝層反演之前，需要通過對整個回波信號的峰值進行光滑樣條回歸并歸一化至[0,1]，可以忽視局部深度范圍內冰層介電常數波動產生的反射，從而提取出整個垂直剖面上的幅值變化趨勢。根據式(1)可知，第1層的差頻信號幅值僅取決于反射系數，所以可以使用冰蓋實測介電常數估算冰蓋表面反射系數，并將其用于幅值變化趨勢的反歸一化變換。通過上述步驟可將回波信號轉換為反射系數曲線，離散化后再使用剝層反演法，遞推求出整個介電常數剖面，進而轉化為密度剖面。

文獻[3]公布的密度ρ與相對介電常數ε之間的經驗公式被廣泛應用于冰蓋密度反演[4,5]，表達式為

通過式(10)將圖2(a)黑線所表示的相對介電常數轉化為密度，再分別使用寬度為4 m的窗口對實測密度(圖2(b)紅線)和由介電常數轉化的密度(圖2(b)藍線)進行滑動平均，兩者結果吻合，證明了式(10)的有效性。

圖2 冰芯B31密度和相對介電常數

實際上差頻信號峰值反映的是冰蓋內部介電常數標準差而不是密度增長率[13]，后者是反應密度隨深度增長過程的重要指標，而剝層反演法將兩者視為等效替代。通過對目前已公開的冰芯數據進行統(tǒng)計分析，兩者呈現(xiàn)相似的變化特征。以冰芯B31為例，分別使用寬度為0.7 m的窗口計算密度-深度變化率和介電常數標準差，并歸一化至[0, 1]后使用光滑樣條回歸提取變化趨勢，結果如圖3所示。可以發(fā)現(xiàn)兩者呈現(xiàn)出相似變化趨勢，即在0～30 m內大幅度降低，在30～70 m內升高后，緊接著又逐漸降低。但是在局部區(qū)域內又存在差異，介電常數標準差在50 m以后明顯高于密度變化率，這會導致50 m以后的剝層反演結果被高估。這屬于方法誤差，因此，對剝層反演得到的初步反演結果進行校準十分必要。

圖3 冰芯B31的介電常數標準差和密度變化率歸一化

3.2 密實化公式

初步反演結果中每個冰層的深度和密度均需要校準，其中包含了大量的未知參數，直接使用模擬退火算法、高斯牛頓法等全局優(yōu)化算法，計算過程過于復雜[14]。本文提出利用冰蓋密實化公式來表示密度變化過程，可有效降低未知參數的數量以提高計算速度，同時還可基于冰蓋密實化原理約束剝層反演得到的初步反演密度，從而達到數據校準的目的。

目前用于表示密度-深度關系的理論主要是基于密實化原理的半經驗公式，并且需要調節(jié)參數將理論密度剖面與實際冰芯數據相匹配。本文基于文獻[15]的密實化原理，對文獻[16]提出的密度-深度公式進行了優(yōu)化，表達式為其中，Ec,Eg和C是常數，g是重力加速度，ρi是純冰密度 (917 kg/m3)，Ts是平均冰蓋表面溫度，R是氣體常數，dz是冰層的厚度，dρ是相較于上一冰層的密度增長量。文獻[16]根據對南極Jurassic和Dyer 高原統(tǒng)計的密度特征對常數C進行擬合，結果表明在密度-深度變化率的第1個極小值點之前為0.03，而之后為0.07。

根據剝層反演法原理可知，方法誤差會隨著深度不斷的增加，這意味著淺層的反演誤差相對較小，因此可以用剝層反演法得到的初步反演結果確定極小值之前的dρ1，再使用優(yōu)化算法尋找出常數A的最優(yōu)解來確定dρ2，最終對dρ進行積分求出校準后的密度剖面。該方法在未依托式(11)中環(huán)境參數及經驗參數的情況下，利用密實化原理完成了數據校準，拓展了本文反演算法的適用范圍。

3.3 優(yōu)化算法

剝層反演法可以得到初步反演密度和密度-深度變化率極小值點所在的深度。通過指數函數對極小值之前的dρ1進行擬合，由于優(yōu)化過程中的計算量較小，本文通過窮舉法列舉常數A在[1, 15]范圍內可能產生的密度剖面，并在其中尋找校準密度的最優(yōu)解。

由3.1節(jié)可知，回波信號被視為密度增長趨勢，通過對回波信號進行積分可以近似的表達密度增長量。在圖3中100 m處密度增長率和介電常數標準差的積分結果相近，這表明雖然初步反演密度與冰芯實測密度之間在局部深度范圍內有所差異，但是兩者在100 m處的密度相近，因此可以將回波信號在0～100 m內的積分差異作為校準密度最優(yōu)解的判定標準。本文將對每個校準密度的回波信號進行仿真，并將其與雷達原始回波信號的積分相對比，當差異最小時則說明校準密度不僅反應了密度增長趨勢，還校正了局部深度范圍內的密度誤差。

為減少計算時間，本文將每個校準密度固定間隔離散化后基于FMCW雷達原理仿真差頻信號。首先根據式(3)、式(4)、式(8)確定每一層的反射系數、傳輸系數和衰減量，再使用式(9)計算每個層的時延，將上述參數代入式(1)求出每個冰層的反射信號，對其累加后得到差頻時域信號，使用FFT進行頻譜估計后提取峰值，對模擬的回波信號峰值進行光滑樣條回歸，歸一化至[0, 1]后進行積分，分別將每個校準密度的積分結果與雷達原始回波信號的積分結果相減后取絕對值，取差異最小的密度曲線作為校準密度的最優(yōu)解。綜上所述，冰蓋垂直剖面上密度和介電常數曲線的反演和校準流程如圖4所示。

圖4 冰蓋密度反演和校準流程圖

4 算法應用

為驗證本文所提出算法的有效性，這里給出一個反演實例，將優(yōu)化算法應用于冰蓋模型和雷達實測數據中。

4.1 理論模型實驗

本文根據圖2(a)冰芯B31在0～100 m內的實測介電常數-深度數據來建立1維水平層狀冰蓋模型[11]，整個垂直剖面上均使用平均電導率23.16 μS/m，且忽略冰的各向異性[17]?；跁r域有限差分方法(Finite Difference Time Domain, FDTD)對差頻時域信號進行仿真，頻譜估計后使用反演算法估算密度剖面，并與圖2(b)中的冰芯實測密度對比。

根據在南極實際使用的FMCW雷達系統(tǒng)參數設置激勵源，電磁波發(fā)射頻率為0.5～2 GHz，掃頻周期為4 ms；吸收邊界條件為完全匹配層(Perfectly Matched Layer, PML)，吸收層寬度為40；圖2(a)介電常數數據的深度分辨率為5 mm，本文設置空間步長dx=2.5 mm；根據Courant穩(wěn)定性條件計算得到時間步長dt=8.65×10–12s。仿真結束后，將接收信號與發(fā)射信號混頻，求出差頻時域信號，通過FFT獲得的回波信號如圖5(a)黑線所示。

使用1.5×10-3的表面反射系數和1.61的表面相對介電常數(320 kg/m3的表面密度)，得到初步密度剖面如圖5(b)中藍線所示。密度-深度變化率極小值點位于深度20 m。通過與窗口寬度為4 m的滑動平均密度(圖5(b)綠線)進行對比，可以明顯觀察到50～90 m深度范圍內的反演密度被高估，這是由于該深度范圍內介電常數標準差大于密度變化率(圖3)，導致均方根誤差(Root Mean Square Error,RMSE)為2.49%。

使用指數函數對0～20 m內的密度進行擬合，以0.1的間隔分別計算常數A在[1, 15]范圍內的密度剖面，并基于式(1)對密度剖面的差頻時域信號進行仿真，頻譜估計后進行光滑樣條回歸并歸一化后積分，如圖5(c)展示了每個回波信號與原始回波信號的積分差異，其中最小值為A=3.7，所代表的密度剖面如圖5(b)紅線所示。尋找到的最優(yōu)密度曲線的回波信號如圖5(d)紅線所示，300～800 ns處的回波信號更低，且50～90 m深度范圍內的反演密度被校正，RMSE降低至1.39%，與冰芯實測密度的相關性為0.999。

圖5 冰芯B31的密度反演結果

4.2 雷達數據實驗

中國第32次南極考察隊使用FMCW雷達探測了南極冰蓋，本文選擇了冰芯LGB69附近的雷達數據來反演平均密度剖面，并與冰芯密度進行對比，結果如圖6。

對雷達數據進行反演之前，至少需要使用40道差頻信號進行相干平均以減少非相干噪聲。使用表面相對介電常數1.82，即415 kg/m3的表面密度，得到初步密度剖面如圖6中黑線所示，其中30～60 m內的密度被低估，而70 m之后的密度被明顯的高估。當A=1.1時的校準密度為最優(yōu)解，如圖6中紅線所示，與冰芯密度(圖6綠線)的RMSE從3.86%降低至2.13%。上述驗證結果表明，本文的反演方法可以有效提升反演精度，校正了方法誤差。盡管密度RMSE僅降低了約1%，但是這對于研究冰蓋的密實化率有著至關重要的作用。

圖6 冰芯LGB69附近的密度反演結果

表面密度或表面介電常數決定了反演的平均密度，因此本文算法依賴表面密度數據。由于在反演算法中對回波信號進行了歸一化處理，當反演不同帶寬的FMCW雷達數據時，差頻信號之間的絕對電平差異被消除，因此反演結果仍能反映冰蓋的平均密度變化趨勢?？紤]到噪聲的存在，對實際雷達數據進行反演前相干平均是必不可少的，對于不同信噪比的雷達系統(tǒng)，可以在相干平均時適當調節(jié)雷達數據的道數以保證反演結果的準確性。

5 結束語

本文基于FMCW雷達差頻信號在層狀介質中的傳播特性和冰蓋的密實化原理，提出了一種將剝層反演法與密實化公式相結合的冰蓋密度反演方法。通過剝層反演法初步得到反演結果，根據冰蓋密實化規(guī)律對初步反演結果進行校準。理論模型和實際測量結果表明，該方法可有效提高反演精度，并校正了剝層反演法中存在的方法誤差，即介電常數標準差與密度變化率之間存在的差異。國內外學者已經獲取了大量的探冰雷達數據，基于本文的反演算法建立數據分析系統(tǒng)可以對雷達數據進行更充分的挖掘與利用，進而獲取整個南極的密度剖面，極大提升冰蓋物質平衡估算精度。