基于压缩并行非局部均值滤波的地质雷达噪声压制算法

崔亚彤，王胜侯，李金伟，吴宇豪，梁思维，马振宁

(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]。基于数据积分算法的高斯欧几里得范数公式如下[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的计算。基于CPU并行策略的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)优于本文列出的其他去噪方法。最后,通过利用并行算法,进一步提高计算效率。因此,本文方法在有效提高计算效率的同时,又提高了噪声压制效果。最后,通过对模型数据和实际数据的处理,验证了该方法的有效性、实用性。