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

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

1 引言

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

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

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

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

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

2 方法

2.1 非局部均值濾波方法

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

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

(1)

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

(2)

(3)

(4)

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

2.2 數(shù)據(jù)壓縮算法

為了提高NLM方法的計(jì)算效率,前人基于數(shù)據(jù)積分算法來加速NLM方法[19]?；跀?shù)據(jù)積分算法的高斯歐幾里得范數(shù)公式如下[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)為積分矩陣,用來計(jì)算差值矩陣s(t,x)的積分,表示為[19]:

(6)

(7)

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

(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 自適應(yīng)濾波參數(shù)估計(jì)算法

NLM方法中的濾波參數(shù)h值決定了NLM方法去噪效果的好壞。常規(guī)NLM方法中,濾波參數(shù)h為全局固定的常數(shù)值,不能保證得到全局最優(yōu)解,這給NLM方法造成了很大的局限性。實(shí)際地質(zhì)雷達(dá)數(shù)據(jù)的有效信號(hào)能量差異大,甚至隨機(jī)噪聲的分布也不是完全隨機(jī)的,因此濾波參數(shù)h的取值通常難以確定。通過最小方差估計(jì)算法可以為含噪數(shù)據(jù)每個(gè)數(shù)據(jù)點(diǎn)估計(jì)濾波參數(shù)h,該方法認(rèn)為參數(shù)h的估計(jì)是噪聲σ的標(biāo)準(zhǔn)差,滿足下式[23]

(11)

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

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

(12)

(13)

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

(14)

本文基于快速自適應(yīng)非局部均值濾波方法,利用數(shù)據(jù)的中心對(duì)稱性與積分算法構(gòu)成了數(shù)據(jù)壓縮算法,極大地提高了常規(guī)NLM方法的計(jì)算效率,通過計(jì)算相似度標(biāo)準(zhǔn)差,實(shí)現(xiàn)了自適應(yīng)濾波參數(shù)估計(jì)算法,有效地提高了去噪效果和計(jì)算效率。

2.4 并行計(jì)算

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

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

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

3 模型數(shù)據(jù)試驗(yàn)

本文結(jié)合常規(guī)NLM方法、中值濾波方法、維納濾波方法與本文快速自適應(yīng)NLM方法進(jìn)行對(duì)比,分別利用簡單模型和復(fù)雜模型,從計(jì)算效率和精度兩方面分別驗(yàn)證本文方法的有效性。

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

表1 模型1的物理幾何參數(shù)

圖1 模型試驗(yàn)1的合成模型數(shù)據(jù)噪聲壓制結(jié)果Fig.1 Noise attenuation results of synthetic model data from model data test 1

圖1(a)為200道無噪聲模型數(shù)據(jù),圖1(b)為含噪數(shù)據(jù),信噪比為1.5,圖1(c)～圖1(f)分別表示中值濾波方法、維納濾波方法、常規(guī)NLM方法與本文快速自適應(yīng)NLM方法的去噪結(jié)果,圖2為去噪數(shù)據(jù)與原始無噪聲數(shù)據(jù)的殘差剖面。上述方法去噪后的均方根誤差(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,利用本文方法去噪結(jié)果的RMS最低,信噪比最高。從圖1所示的去噪結(jié)果可以看出以上四種方法均可有效壓制隨機(jī)噪聲,但從圖2殘差剖面看出,利用中值濾波方法與維納濾波方法去噪,容易對(duì)有效信號(hào)的能量造成一定損失,而常規(guī)NLM方法和本文方法相對(duì)表現(xiàn)較好,在有效去除噪聲的同時(shí),減少了對(duì)有效信號(hào)的損傷。

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

表2 不同去噪方法的mSNR和mRMS對(duì)比

圖2 模型試驗(yàn)1的殘差剖面Fig.2 Residual profile in model test 1

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

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

圖3 模型試驗(yàn)2的合成模型數(shù)據(jù)噪聲壓制結(jié)果Fig.3 Noise attenuation results of synthetic model data from model data test 2

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

圖5 實(shí)際數(shù)據(jù)去噪結(jié)果Fig.5 Noise attenuation results of actual data

4 實(shí)際數(shù)據(jù)試驗(yàn)

為了驗(yàn)證本文方法的實(shí)用性、有效性,將本文方法應(yīng)用于實(shí)際大型地質(zhì)雷達(dá)數(shù)據(jù)的隨機(jī)噪聲壓制處理中。

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

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

5 結(jié)論

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