李 進，王小明，譚 磊

(1.南京市水利規(guī)劃設(shè)計院股份有限公司，江蘇 南京 210022；2.蘇交科集團股份有限公司，江蘇 南京 210000；3.浙江省水利河口研究院(浙江省海洋規(guī)劃設(shè)計研究院)，浙江 杭州 310020；4.浙江省水利防災(zāi)減災(zāi)重點實驗室，浙江 杭州 310020)

1 引 言

地球物理勘探是通過采集到的多物理場數(shù)據(jù)體來認識地下空間的地層構(gòu)造、礦產(chǎn)賦存、災(zāi)害分布以及工程地質(zhì)環(huán)境等[1-4]，具有高效、無損、透視、大數(shù)據(jù)的特點，在宏觀調(diào)查和精細勘探中都有著突出的優(yōu)勢，與鉆探手段相比，二維、三維的探測成果能連續(xù)的反映出勘探區(qū)的地質(zhì)情況，極大地提高了勘探密度[5,6]。但是，受限于野外工作環(huán)境、地形條件以及時間成本等因素的制約，在權(quán)衡精度與效率利弊的情況下，提交的物探圖像往往是適度特征點數(shù)據(jù)經(jīng)網(wǎng)格插值處理后的成果。目前，不規(guī)則離散數(shù)據(jù)網(wǎng)格化運算的插值技術(shù)較多，如美國Golden Software公司的商業(yè)Surfer系列軟件中就有12種不同類型的插值方式[7]，豐富的插值功能為數(shù)據(jù)顯示提供了便利，但具體到每種插值算法在數(shù)據(jù)可視化表達方面也存在一定的差異性。近年來，一些物探學者也針對不用插值算法的應(yīng)用及特點方面作出許多有益的工作。吳衛(wèi)國[8]利用Surfer軟件的克里金法對水系沉積物數(shù)據(jù)進行數(shù)據(jù)網(wǎng)格化，并結(jié)合白化功能減少離散數(shù)據(jù)的外推；劉曉霞等[9]推導了Kriging球面應(yīng)變計算公式，交叉檢驗及殘差分析結(jié)果表明，速度場濾波和網(wǎng)格插值是可效的；姚文等[10]以某磁鐵礦區(qū)高精度磁測數(shù)據(jù)為例，在比較了5種網(wǎng)格化平面圖的基礎(chǔ)上，認識到在保證網(wǎng)格化的精度和等值線圓滑的前提下，克里金插值法和最小曲率法網(wǎng)格化效果最好；王咸彬等[11]在總結(jié)了常見散亂數(shù)據(jù)插值方法的基礎(chǔ)上，采用克里金插值方法把實測測井數(shù)據(jù)插值成反演的三維速度場初始模型；郭崇光[12]采用Surfer軟件中的克里格方法對電測深數(shù)據(jù)進行網(wǎng)格化，并指出具體插值中的注意點。

高可靠度的原始數(shù)據(jù)是地球物理探測數(shù)據(jù)網(wǎng)格化成像的基礎(chǔ)，然而在觀測過程中顯然不能排除外來噪聲的干擾[13]，因此有必要采用一定的手段對網(wǎng)格數(shù)據(jù)體進行濾波降噪處理。李文杰[14]運用二維自動調(diào)平處理了頻率域航空電磁數(shù)據(jù)的零漂誤差，提高視電阻率和平面剖面圖質(zhì)量；余金煌等[15]采用小子域濾波技術(shù)對高密度電法探測的拋石層圖像進行增強處理，提高了圖像精度和可讀性；黃基文[16]根據(jù)基階模態(tài)與高階模態(tài)面波之間視速度存在差異，利用二維數(shù)字濾波法進行分離；王瑞雪等[17]把均值濾波、中值濾波以及小波變換對測井電導率數(shù)據(jù)進行濾波處理。近年來，高密度電法在水庫大壩滲漏隱患探測中得到廣泛應(yīng)用[18]，但較少研究不同插值方法對電阻率圖像上隱患的優(yōu)化識別，并且涉及到具體濾波方法的差異性以及相應(yīng)參數(shù)設(shè)置等內(nèi)容還需進一步探討。為此，本文以水庫大壩蟻穴滲漏高密度電法探測數(shù)據(jù)為試驗對象，比較了多種網(wǎng)格插值方法圖像的差異性，并對自然鄰點法插值圖像中的濾波方法及插值結(jié)果進行分析，從而得到適合于土石壩電阻率可視化表達的成像技術(shù)，以此提高隱患的識別能力及效果。

2 試驗工程概況

浙江某水庫始建于1957年10月，并于1958年7月建成蓄水，壩址以上集雨面積0.301 km2，總庫容為16.395萬m3，正常庫容為12.149萬m3，灌溉面積720畝，是一座以灌溉為主的小(2)型水庫。大壩壩型為類均質(zhì)壩，壩頂高程為55.70 m，壩頂長為103.5 m，最大壩高7.58 m，壩頂寬3.0 m，壩頂鋪廣場磚。

該水庫已服役運行60多年，大壩防滲結(jié)構(gòu)存在不同程度的老化，盡管多次對大壩進行加固處理，但在極端天氣作用下，于2019年在大壩左壩段出現(xiàn)滲漏、塌陷的突發(fā)險情，為盡早排查出大壩滲漏的原因以期保障水庫的安全運行，采用高密度電法對水庫大壩滲漏原因進行探測，并對異常區(qū)進行防滲處理，在灌漿過程中大壩迎水坡、背水坡均出現(xiàn)局部的冒漿現(xiàn)象，與探測得到的滲漏通道走向基本吻合。以此探測數(shù)據(jù)為基礎(chǔ)，從不同插值方法、濾波分析等角度，提升對探測成果的可識別性。

圖2是地下記錄視電阻率值的散點分布圖，記錄點數(shù)據(jù)共有946個，最小電阻率值為11.98 Ω·m，最大電阻率值339.15 Ω·m，均值95.93 Ω·m，中值86.39 Ω·m。

圖1 試驗現(xiàn)場示意Fig.1 Schematic diagram of the test site

圖2 原始數(shù)據(jù)記錄點分布Fig.2 Distribution of original data recording points

3 電阻率數(shù)據(jù)插值及濾波處理方法

3.1 常用數(shù)據(jù)插值方法的簡介

數(shù)據(jù)的網(wǎng)格化處理工作是空間離散測點原始數(shù)據(jù)可視化成像的基礎(chǔ)，在高密度電法勘探過程中，水平方向上各個電極測點的最小間距是保持固定的，垂直方向上相鄰深度點的距離又與水平長度有關(guān)，并且隨著深度不斷的加大，水平層間的測點呈等差數(shù)列不斷減少，總體表現(xiàn)為矩陣排列的形式。因此，為把縱、橫向上數(shù)據(jù)點之間的變化趨勢反映出來，對原始數(shù)據(jù)采用網(wǎng)格化插值是重要的環(huán)節(jié)。在地球物理數(shù)據(jù)網(wǎng)格插值方面應(yīng)用較多的方法[19]，通常有反距離加權(quán)插值法(Inverse Distance to a Power)、克里金插值法(Kriging)、最小曲率插值法(Minimum Curvature)、改進謝別德插值法(Modified Shepard’s Method)、自然鄰點插值法(Natural Neighbor)、最近鄰點插值法(Nearest Neighbor)、多項式回歸插值法(Polynomial Regression)、徑向基函數(shù)插值法(Radial Basis Function)、線性插值的三角剖分插值法(Triangulation with Linear Interpolation)、移動平均插值法(Moving Average)、數(shù)據(jù)度量插值法(Data Metrics)、局部多項式插值法(Local Polynomial)12種方法。同一組數(shù)據(jù)采用不同的插值方法，圖像中有用信息的顯示程度也存在差異，因此有必要通過比較選擇最優(yōu)的插值方法[20]，從而指導工程應(yīng)用中的數(shù)據(jù)處理。

3.2 常用網(wǎng)格濾波方法

數(shù)據(jù)網(wǎng)格化處理有效預(yù)測了測點之間的數(shù)據(jù)體，但不能消除原始信號中的噪聲信息。高密度電法數(shù)據(jù)噪聲來源復雜，包括天然電場、游散電流、接地不良、地形起伏以及儀器雜音等[21]，外來噪聲往往會降低數(shù)據(jù)信號的信噪比，從而降低圖像信息的真實度以及圖像自動識別的能力。因此，開展網(wǎng)格化數(shù)據(jù)體進行濾波工作對水庫大壩隱患的識別十分重要，但如何選擇適合于電阻率成像的濾波方法成為圖像可視化的關(guān)鍵。常用的濾波器分為線性濾波器和非線性濾波器，并且線性濾波器中有滑動平均濾波、距離加權(quán)平均濾波、反距離濾波、高斯低通法等。

4 結(jié)果分析

4.1 不同插值方法的網(wǎng)格化結(jié)果及分析

圖3是12種插值方法網(wǎng)格化后的電阻率圖像，所有圖像采用相同的色階歸一化。從圖像上可以直觀看出，不同的插值方法得到的電阻率值圖像存在較大的差異，主要表現(xiàn)在經(jīng)插值后的電阻率范圍、分布以及圖像的平滑度等。其中，多項式回歸插值法、移動平均插值法、數(shù)據(jù)度量插值法以及局部多項式插值法平滑度較大，不能有效揭示出大壩的電阻率值分布，不適合運用在水庫滲漏探測的高密度電法數(shù)據(jù)處理中；反距離插值圖像上出現(xiàn)較多的“牛角眼特征”，不利于對隱患的識別；改進的謝別德插值對圖像的邊界處理較差，局部出現(xiàn)負值，不符合電阻率的物理意義；最近鄰點法插值在圖像上形成連續(xù)的網(wǎng)格，網(wǎng)格之間的連續(xù)性較差；最小曲率插值的視電阻率范圍為8.29～518.94 Ω·m(表1)，比實際測量最大值339.15 Ω·m高約34.6%；線性插值的三角剖分插值中的電阻率值也出現(xiàn)負值?？死锝鸱?、自然鄰點法以及徑向基函數(shù)插值法圖像電阻率等值線平滑度較好，插值后的電阻率值與原始值范圍基本吻合，在測線上18～24 m段、核心深度3 m的高阻閉合區(qū)域與隱患位置基本一致，總體上能反映出大壩地下結(jié)構(gòu)的分布特征，也表明高密度電法在探測土石壩隱患中結(jié)果可靠度高。但是，自然鄰點法插值外延性較差，克里金法、徑向基函數(shù)插值能向無值區(qū)外推出電阻率的變化趨勢，但計算時間較長。因此，水庫滲漏高密度電法原始數(shù)據(jù)的插值優(yōu)選采用自然鄰點法，其次是克里金插值法、徑向基函數(shù)插值法，在工程應(yīng)用中采用合理的插值算法能更加真實地反映出大壩的電阻率變化特征。

圖3 不同插值方法的電阻率Fig.3 Resistivity plots of different interpolation methods

表1 不同插值方法的電阻率范圍

4.2 不同濾波方法的等值線結(jié)果及分析

地下空間的地電場信號頻率較低，因此應(yīng)采用低通數(shù)字濾波器對電阻率圖像進行噪聲處理。根據(jù)多種插值算法網(wǎng)格化成果的比較，采用自然鄰點插值的電阻率值進行濾波比較分析。圖4(a)是不進行濾波處理的等值線分布圖，在圖像上存在局部極值的小圈閉，并且線條平滑度較差；經(jīng)濾波處理后的圖像，圖中的噪聲信號明顯得到降低，滑動平均、距離加權(quán)法的濾波效果相對更好。

圖4 不同濾波方法的等值線Fig.4 Contours of different filtering methods

從濾波后的電阻率值范圍(表2)上來看，濾波之前的電阻率范圍為23.01～336.22 Ω·m，經(jīng)濾波之后高、低極值得到削減，其中滑動平均濾波電阻率范圍相對更收斂，有效降低了噪聲對數(shù)據(jù)的影響。

表2 不同濾波方法的電阻率范圍

4.3 不同過濾器尺寸的濾波結(jié)果及分析

數(shù)據(jù)濾波的效果不僅與濾波器的選擇有關(guān)，而且過濾器尺寸的大小也會對濾波效果產(chǎn)生影響。圖5是在垂向-3 m的橫向上以及水平24 m縱向上的電阻率變化圖。由圖5可知，過濾器尺寸越小，濾波效果越差，在曲線上的極值點處變化幅度越大，但隨著尺寸的增大，曲線的平滑度增加，從而導致相鄰異常區(qū)的邊界變得模糊。因此，在對水庫滲漏的高密度電法數(shù)據(jù)采用滑動濾波時，過濾器的尺寸7×7。

圖5 不同過濾器尺寸的電阻率變化Fig.5 Resistivity variations of different filter sizes

4.4 不同迭代次數(shù)的濾波結(jié)果及分析

濾波的迭代次數(shù)大，則可以提高等值線的平滑程度，但過于采用多次濾波計算也可能降低對異常區(qū)的識別。圖6是在垂向-3 m的橫向上以及水平24 m縱向上的電阻率變化圖。從圖像上可以看出，當?shù)螖?shù)為1時，電阻率曲線在極值點變化強烈，同時在曲線上有所偏離；當?shù)螖?shù)為7時，電阻率曲線明顯過于圓滑，在拐點位置電阻率的幅度過小，對異常區(qū)邊界的確定帶來不利影響。綜上所述，在對水庫滲漏的高密度電法數(shù)據(jù)采用滑動濾波時，迭代次數(shù)采用4次。

圖6 不同迭代次數(shù)的電阻率變化Fig.6 Resistivity variation of different number of iterations

5 結(jié) 論

1)高密度電法數(shù)據(jù)經(jīng)網(wǎng)格化插值處理后能夠把測量點可視化顯示出來，有利于評價地質(zhì)體的電阻率分布，建議在網(wǎng)格化插值中優(yōu)先選用自然鄰點法，其次是克里金插值法、徑向基函數(shù)插值法；

2)不同濾波方法、過濾器尺寸以及迭代次數(shù)對網(wǎng)格化數(shù)據(jù)的去噪效果存在差異，水庫大壩隱患的電阻率濾波優(yōu)先選用滑動平均法，過濾器尺寸的行列均為7，迭代次數(shù)為4次；

3)水庫大壩隱患類型復雜，電阻率等值線形態(tài)及大小也存在差異，在具體使用時，應(yīng)對不同電阻率分布特征的可視化方法進行綜合比較，從而提高對不良隱患的識別能力。