基于壓縮并行非局部均值濾波的地質雷達噪聲壓制算法

崔亞彤，王勝侯，李金偉，吳宇豪，梁思維，馬振寧

(1.天津市勘察設計院集團有限公司,天津 300191;2.中國地質大學 資源學院,湖北 武漢 430074;3.中國地質大學 地球物理與信息技術學院,北京 100083)

1 引言

地質雷達(Ground Penetrating Radar,GPR)作為一種探測充分、操作簡單、采樣快速、分辨率高的地球物理技術,被廣泛地用于各種地質目標體的探測,包括城市道路空洞和路面層厚度估計[1,2]、考古調查[3]、冰川內部結構的調查[4]、污染物環境調查[5]等方面。

地質雷達的發射天線和接收天線通常與地面接觸,接收到的反射波用于識別目標地質體的埋深、大小、形狀以及賦存狀態。地質雷達數據中存在多種類型的噪聲。其中,隨機噪聲是由在任意時間和頻率上的隨機波動造成的。隨著探測深度的增加,信號被極大地減弱,當雷達數據受到隨機噪聲嚴重污染時,有效信號將被大量掩蓋。因此,噪聲壓制在數據處理中起著至關重要的作用,去噪效果將直接影響后續地質雷達資料的解釋工作。

為了有效壓制隨機噪聲,提高信噪比,國內外許多學者提出不同的方法,如曲波變換[6],拉東變換[7]、維納濾波[8]、中值濾波[9]、小波變換[10,11]、非局部均值濾波(Non-Local Mean Filtering,NLM)方法[12]等,以及基于人工智能的方法,如數據驅動緊密框架[13]和機器學習[14]。NLM方法并不基于線性假設,因此在處理彎曲同相軸時,該方法可以有效地保護有效信號,壓制隨機噪聲。常規NLM方法起源于圖像的隨機噪聲壓制處理[15],該方法與其他濾波方法相比,計算效率較低,尤其在處理大型地質雷達數據時,效率低下。

前人提出了并行分塊NLM法[16]、基于隨機投影算法的NLM法[17]、變窗口的快速NLM法[18]等快速算法,來提高計算效率。但由于前述方法為了增加計算效率,在計算過程中并沒有完全遍歷每個數據點,因此可能會損失部分有效信號,降低計算精度。Froment提出的基于數據積分算法的快速NLM法[19]可以在提高計算效率的同時,計算所有數據點,因此相較于前述方法,避免了犧牲精度的可能。常規NLM方法在處理實際數據時,所選擇的濾波參數值通常為一常數,為了進一步提高去噪效果,針對不同區域可選擇不同的濾波參數來提高去噪效果[20,21],但會明顯增加計算量。

本文給出了一種基于非局部均值濾波的自適應數據壓縮并行去噪算法,可快速且有效地壓制地質雷達隨機噪聲。本文利用一種中心對稱數據積分算法以及自適應調整濾波參數分布算法,在提高計算精度的同時,有效地降低了計算成本。最后,通過對簡單模型數據、復雜模型數據和美國地質調查局公布的實際數據進行去噪試驗,驗證了該方法的可行性、有效性。

2 方法

2.1 非局部均值濾波方法

噪聲數據Dnoise(t,x)可由下式表示:

Dnoise(t,x)=Dtrue(t,x)+n

(1)

其中,Dtrue(t,x)表示無噪聲數據,n并表示隨機噪聲。Ddenoise(t,x)表示去噪后的數據,通過對Dnoise(t,x)進行加權平均計算即可得到每個點處的Ddenoise(t,x)數據。因此,Ddenoise(t,x)可以表示為[22]:

(2)

(3)

(4)

假設Dnoise(t,x)總點數為N=NtNx,鄰域窗口之間的相似度計算量為(2ds+1)2,搜索窗內計算(2Ds+1)2次相似度,NLM濾波的計算復雜度為O(N(2Ds+1)2(2ds+1)2),計算量巨大。因此,常規NLM方法在處理大型實際地質雷達數據時,會明顯受到計算效率的制約。

2.2 數據壓縮算法

為了提高NLM方法的計算效率,前人基于數據積分算法來加速NLM方法[19]?；跀祿e分算法的高斯歐幾里得范數公式如下[19]:

[St(t+ds,x+ds)+St(t-ds-1,x-ds-1)-

St(t+ds,x-ds-1)-St(t-ds-1,x+ds)]

(5)

S(t1,x1)為積分矩陣,用來計算差值矩陣s(t,x)的積分,表示為[19]:

(6)

(7)

為了進一步降低計算復雜度和提高計算效率,本文利用差值矩陣s(t,x)的中心對稱性來計算高斯歐幾里得范數,即

(8)

s(t-r1,x-r2)[+r]=‖Dnoise(t-r1,

(9)

[St(t+ds-r1,x+ds-r2)[+r]+

St(t-ds-r1-1,x-ds-r2-1)[+r]-

St(t+ds-r1,x-ds-r2-1)[+r]-

St(t-ds-r1-1,x+ds-r2)[+r]]

(10)

2.3 自適應濾波參數估計算法

NLM方法中的濾波參數h值決定了NLM方法去噪效果的好壞。常規NLM方法中,濾波參數h為全局固定的常數值,不能保證得到全局最優解,這給NLM方法造成了很大的局限性。實際地質雷達數據的有效信號能量差異大,甚至隨機噪聲的分布也不是完全隨機的,因此濾波參數h的取值通常難以確定。通過最小方差估計算法可以為含噪數據每個數據點估計濾波參數h,該方法認為參數h的估計是噪聲σ的標準差,滿足下式[23]

(11)

其中,V0代表無噪聲數據,濾波參數矩陣h2的估計即可表示為[23]

h2(t,x)≈σ2=min(E‖V(t,x)-

(12)

(13)

其中,TDs(t,x)代表搜索窗口Ω={i,j∈Ω:|i-t|≤Ds,|j-x|≤Ds}范圍內的相似度標準差,即

(14)

本文基于快速自適應非局部均值濾波方法,利用數據的中心對稱性與積分算法構成了數據壓縮算法,極大地提高了常規NLM方法的計算效率,通過計算相似度標準差,實現了自適應濾波參數估計算法,有效地提高了去噪效果和計算效率。

2.4 并行計算

傳統NLM遍歷搜索框內數據的計算是相互重疊的,因此難以應用并行計算。在串行通信中,由于只有一條鏈路,由于使用單鏈路,并且在傳輸和接收數據時不太復雜,串行通信更加簡單和容易執行。由于數據傳輸方式單一,因此計算所需時間較長。在并行通信中,數據在多個鏈路上同時傳輸、計算。這種類型的通信由于使用了多條鏈路,可以同時進行數據計算,因此傳輸所需的時間比串行通信要短得多。

并行計算主要分為CPU(Central Processing Unit,CPU)并行與GPU(Graphics Processing Unit,CPU)并行,其中CPU并行具有共享內存并行編程(Openmulti-processing,OpenMP)、消息傳遞接口(Message Passing Interface,MPI)等。Matlab并行計算提供了數據并行性和任務并行性這兩個功能,可以在更高的級別上執行數據和任務并行算法,而不需要為特定的硬件和集群架構編寫程序[24]。

本文并行算法是基于OpenMP的CPU并行策略,采用Matlab Parfor并行循環進行計算。本文將不同搜索距離的數據進行分塊并記錄在索引矩陣中,以保證數據積分的計算為非重疊的;根據索引矩陣,將不同分塊的數據分配給10個CPU核心進行并行計算;最后所有的計算完成后,通過索引矩陣合并所有計算出的歐幾里得距離,進而同時計算每個數據點對應的加權平均,實現NLM的計算?；贑PU并行策略的NLM算法,通過把每個搜索窗口的處理放到不同的線程上執行,對各個搜索窗口同時進行處理,計算、時間復雜度大大降低。

3 模型數據試驗

本文結合常規NLM方法、中值濾波方法、維納濾波方法與本文快速自適應NLM方法進行對比,分別利用簡單模型和復雜模型,從計算效率和精度兩方面分別驗證本文方法的有效性。

gprMax是一個模擬電磁波傳播的開源軟件,該軟件利用時域有限差分(Finite Difference Time Domain,FDTD)方法求解麥克斯韋方程組,實現對模型的正演[25-28]。數據試驗1為道路空洞正演模擬,物性與幾何參數見表1,激發源采用雷克子波,中心頻率為1 000 MHz,時窗為30 ns,道間距為0.02 m,模擬200道雷達數據,數據大小為1 024×200。

表1 模型1的物理幾何參數

圖1 模型試驗1的合成模型數據噪聲壓制結果Fig.1 Noise attenuation results of synthetic model data from model data test 1

圖1(a)為200道無噪聲模型數據,圖1(b)為含噪數據,信噪比為1.5,圖1(c)～圖1(f)分別表示中值濾波方法、維納濾波方法、常規NLM方法與本文快速自適應NLM方法的去噪結果,圖2為去噪數據與原始無噪聲數據的殘差剖面。上述方法去噪后的均方根誤差(Root Mean Square,RMS)分別為0.28、0.24、0.23、0.21,信噪比(Signal to Noise Ratio,SNR)分別為24.86、30.49、31.80、32.05,利用本文方法去噪結果的RMS最低,信噪比最高。從圖1所示的去噪結果可以看出以上四種方法均可有效壓制隨機噪聲,但從圖2殘差剖面看出,利用中值濾波方法與維納濾波方法去噪,容易對有效信號的能量造成一定損失,而常規NLM方法和本文方法相對表現較好,在有效去除噪聲的同時,減少了對有效信號的損傷。

為了對比不同噪聲水平數據的去噪效果,模型1數據信噪比分別分布于0.4～0.6,1.4～1.6,2.4～2.6之間,利用中值濾波方法、維納濾波方法、常規NLM方法與本文快速自適應NLM方法進行去噪,每種方法各執行5次,計算得到平均均方根誤差(Mean Root Mean Square,mRMS)和平均信噪比(Mean Signal to Noise Ratio,mSNR),結果如表2所示。通過多次數據試驗發現,利用本文方法去噪結果的mSNR高于其他方法,mRMS低于其他方法。

表2 不同去噪方法的mSNR和mRMS對比

圖2 模型試驗1的殘差剖面Fig.2 Residual profile in model test 1

Philipp與Tronicke使用一組已公布的河流—冰川沉積物含水層模擬研究的水相數據,建立了一個三維沉積模型,采用赫茲偶極子源激發,中心頻率為100 MHz,時窗為200 ns,道間距為0.05 m,并將該模型進行正演,得到了一個三維沉積模型地質雷達響應數據[29]。模型數據試驗2為該三維地質雷達數據中的一個二維切片數據。

該切片數據大小為4 156×310,圖3(a)為310道無噪聲模型數據,圖3(b)為去初至后含噪數據,其信噪比為2.5,圖3(c)～圖3(f)分別展示了中值濾波方法、維納濾波方法、常規NLM方法和本文方法去噪后的結果,相較于其他方法,本文方法能夠較好地恢復深層信息。上述方法去噪后的RMS分別為0.306 7、0.318 2、0.232 1、0.204 6,信噪比分別為26.89、27.42、28.43、29.02,利用本文方法去噪結果的RMS最低,信噪比最高。模型數據試驗1和試驗2的結果表明,本文方法除在殘差剖面中泄漏少量能量外,在壓制隨機噪聲和保留有效信號方面具有較好的效果。

圖3 模型試驗2的合成模型數據噪聲壓制結果Fig.3 Noise attenuation results of synthetic model data from model data test 2

最后,本文利用一系列不同大小的數據對比不同方法的計算時間。本試驗所使用的計算機處理器為英特爾Corei9-10900K,內存為32 G,處理器數量10,圖4顯示了本文方法和常規NLM方法的計算時間對比。當模型數據量從64×64增加到1 024×1 024時,常規NLM方法的計算時間隨數據量的增加呈指數增加,而本文方法在計算效率上明顯高于常規NLM方法,計算效率高。

圖5 實際數據去噪結果Fig.5 Noise attenuation results of actual data

4 實際數據試驗

為了驗證本文方法的實用性、有效性,將本文方法應用于實際大型地質雷達數據的隨機噪聲壓制處理中。

2013年4月13日至20日,美國地質調查局圣彼得堡海岸和海洋科學中心在阿拉巴馬州多芬島進行了地球物理和沉積物采樣調查。該研究的目的是量化大陸沼澤和河口環境的無機物和有機物的加積率。為了實現這一研究目的,研究人員使用了各種現場和實驗室方法,包括使用地質雷達進行地下成像、沉積物取樣、巖性和微化石分析等,通過上述定量、定性數據,對該研究區域進行地質演化推算[30]。

本文采用的實際數據為其中一個地質雷達剖面,該數據大小為381×2 300,為2 300道數據,去噪結果如圖5所示。從圖5(a)可以看出,隨機噪聲污染了實際地質雷達資料,使得地層模糊,斷層構造不清晰,難以進行后續的解釋工作。圖5(b)～圖5(c)展示了常規NLM方法與本文方法的去噪結果,圖5(d)～圖5(e)分別為以上兩種方法的殘差剖面。常規NLM方法在殘差剖面中殘留部分有效信號,且運行時間約3.7 min,本文方法處理結果中的構造信息清晰,有效信號的信息得到了很好的保存,能夠有效壓制隨機噪聲,計算時間僅為9.7 s,計算效率也有明顯的提高。

5 結論

本文給出了一種基于非局部均值濾波的自適應數據壓縮并行去噪算法,用于地質雷達數據隨機噪聲壓制。本文利用數據積分算法與數據的中心對稱性算法,加速傳統NLM濾波方法,提升計算效率,計算量由O(N(2Ds+1)2(2ds+1)2)降低至O(N(2Ds+1)2/2)。其次,本文采用自適應濾波參數估計算法,提高了噪聲壓制效果,采用本文方法去噪結果的平均均方根誤差(mRMS)和平均信噪比(mSNR)優于本文列出的其他去噪方法。最后,通過利用并行算法,進一步提高計算效率。因此,本文方法在有效提高計算效率的同時,又提高了噪聲壓制效果。最后,通過對模型數據和實際數據的處理,驗證了該方法的有效性、實用性。