張一凡,劉海飛,柳建新,郭 鵬,劉 昕

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院,長沙 410083;2.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,長沙 410083)

0 引言

地球物理觀測的數(shù)據(jù)都可視為空間、時(shí)間域離散的數(shù)字信號(hào)。在點(diǎn)距稀疏不均勻(如航空物探)、因客觀條件造成數(shù)據(jù)缺失(如地形因素)、儀器設(shè)備不穩(wěn)定造成部分壞記錄的情況下,往往需要通過有效插值方法盡量規(guī)律性補(bǔ)足觀測數(shù)據(jù)集。另外,為避免數(shù)據(jù)解釋得到錯(cuò)誤的認(rèn)識(shí),強(qiáng)噪音引起的異常數(shù)據(jù)需在數(shù)據(jù)處理中剔除,并盡量采用一種保光滑、保形狀的插值方法對(duì)真實(shí)數(shù)據(jù)進(jìn)行恢復(fù)。線性插值、樣條插值等常規(guī)插值方法在地球物理數(shù)據(jù)處理中應(yīng)用廣泛,但對(duì)于一些特殊場景,其插值光滑效果還需改進(jìn)。

MQ函數(shù)(Multi-quadric function)最早由Hardy[1]于1971年提出,Franke[2]發(fā)現(xiàn)在29類散亂數(shù)據(jù)插值中,MQ函數(shù)的精確性、有效性、穩(wěn)定性優(yōu)于其它方法,然而利用MQ函數(shù)需要求解線性方程組且系數(shù)矩陣通常是病態(tài)的。因此,對(duì)MQ函數(shù)進(jìn)行某種近似且避免求解線性方程組的方法逐漸引起重視,這種逼近函數(shù)被稱為擬MQ函數(shù)。Beatson 和 Powell[3]利用 MQ 函數(shù)提出了三種單變量的擬MQ插值算子QA、QB、QC,并得到了它們的誤差,同時(shí)證明擬插值算子QC具有很好的線性再生性和插值精度。Wu和Schaback[4]對(duì)QC進(jìn)行改進(jìn),提出了一種新的擬MQ插值算子QD,并討論了它們的逼近階和保形性。后來,該插值算子被廣泛應(yīng)用于求解偏微分方程[5-9]。此外,Ling[10]在Wu和Schaback工作的基礎(chǔ)上,提出了一種多水平單變量擬MQ插值算子。Ling[11]還利用空間維數(shù)分離方法,構(gòu)造了擬MQ插值算子QR。Jiang等人[12]在擬插值算子QD的基礎(chǔ)上構(gòu)造了一種高精度的擬插值算子QW,該算子具有保線性且具有很高的精度。高文武[13]討論了基于導(dǎo)數(shù)信息的擬MQ插值的構(gòu)造方法,并用該方法構(gòu)造了一個(gè)擬MQ插值算子,證明了該算子的收斂性及保形性。

擬MQ函數(shù)無需求解線性方程組,根據(jù)離散節(jié)點(diǎn)能快速直接給出插值結(jié)果,還具有良好的保形性、計(jì)算穩(wěn)定、計(jì)算量小等優(yōu)點(diǎn)。因此,在過去的幾十年里,擬MQ函數(shù)取得了一定的發(fā)展[14-16]。但整體而言,該函數(shù)在地球物理領(lǐng)域的應(yīng)用研究較少。

筆者利用Wu和Schaback[4]提出的擬MQ插值算子QD,分析整理了其在物探數(shù)據(jù)處理中的應(yīng)用,主要包括擬合地表高程、高密度電法數(shù)據(jù)插值與平滑、電測深曲線插值與求導(dǎo)等,取得了良好的處理效果,證實(shí)了該函數(shù)在物探數(shù)據(jù)處理中的有效性和實(shí)用性。

1 基本理論

1.1 MQ函數(shù)

f(xj),j=0,1,…,n

(1)

式中:φ(x-xi)為MQ函數(shù),對(duì)于一維插值問題,其形式為:

(2)

式中:s為形狀參數(shù),決定了曲線的銳度。若相鄰插值結(jié)點(diǎn)的最大間隔為h,s相對(duì)h保持相關(guān)且趨于零可獲得更高的準(zhǔn)確性[1]。將式(1)寫成矩陣形式:

Ac=z

(3)

其中:

A=

求解線性方程(3),即可求得解向量cn×1。將解向量cn×1及待插點(diǎn)xj代入式(1),便可計(jì)算出待插點(diǎn)xj的值f(xj)。

1.2 擬MQ函數(shù)

為避免求解線性方程組,Beatson 和 Powell[3]利用 MQ 函數(shù)提出了三種單變量擬插值算子QA、QB、QC。Wu和 Schaback[4]在此基礎(chǔ)上給出了逼近階和保形性較好的擬插值算子QD。

(4)

擬插值算子QB滿足:

f(xn)βn(x)

(5)

f(xn)αn(x)

(6)

為了進(jìn)一步提高插值精度,相較于擬插值算子QA和QB,還可以引入更多條件定義擬插值算子QC:

QCf(x)=QBf(x)+f′(x0)γ0(x)+

f′(xn)γn(x)

(7)

將其改寫為:

QCf(x)=f(x0)β0(x)+f′(x0)γ0(x)+

f(xn)βn(x)

(8)

Beatson 和 Powell[3]對(duì)三種擬插值算子進(jìn)行了大量測試,表明QC的插值精度最高。但應(yīng)用QC需要事先知道邊界點(diǎn)的導(dǎo)數(shù),實(shí)際應(yīng)用中很難滿足此條件。為克服這一缺點(diǎn),Wu和Schaback[4]對(duì)其進(jìn)行改進(jìn),定義了一種具有保單調(diào)性和保凸性且不需要知道邊界點(diǎn)的導(dǎo)數(shù)值的擬插值算子QD,它的逼近效果與QC相當(dāng)。擬插值算子QD的定義如下:

QDf(x)=f(x0)a0(x)+f(x1)a1(x)+

f(xn-1)an-1(x)+f(xn)an(x)

(9)

可將式(9)展開為:

(10)

由于大量實(shí)際工程數(shù)據(jù)均無法給出插值端點(diǎn)處的導(dǎo)數(shù)值,相比之下QD具有更好的應(yīng)用價(jià)值。

1.3 擬MQ函數(shù)的導(dǎo)數(shù)形式

常規(guī)插值方法往往只保證函數(shù)值插值精度,但很多地球物理方法需要對(duì)數(shù)據(jù)的數(shù)值微分進(jìn)行分析以突出弱異常,這對(duì)插值方法提出了更高的要求。對(duì)于式(10),其一階導(dǎo)數(shù)為:

QDf′(x)=

(11)

二階導(dǎo)數(shù)為:

(12)

通過式(11)、式(12),便能推算出離散點(diǎn)的導(dǎo)數(shù)值,從而在觀測曲線中發(fā)揮提取弱異常的優(yōu)勢。

2 擬MQ函數(shù)在物探數(shù)據(jù)處理中的應(yīng)用

2.1 在地表高程插值中的應(yīng)用

以山東某金礦地球物理勘查剖面的地表高程數(shù)據(jù)為例,該剖面已知的地表高程點(diǎn)有64個(gè),橫向分布范圍為[0～4 460],并且呈不等間距分布。以10 m為間隔從0到4 460均勻采樣作為待插結(jié)點(diǎn),采用三次樣條函數(shù)及擬MQ函數(shù)對(duì)地表高程點(diǎn)分別進(jìn)行插值。插值結(jié)果如圖1、圖2所示。

圖1 三次樣條函數(shù)的插值結(jié)果Fig.1 Interpolation results of cubic spline functions

圖2 擬MQ函數(shù)的插值結(jié)果Fig.2 Interpolation results of quasi-MQ functions

可以看出:在利用三次樣條函數(shù)進(jìn)行插值時(shí),當(dāng)函數(shù)值變化較大且插值結(jié)點(diǎn)較稀疏時(shí),會(huì)在插值結(jié)點(diǎn)附近引起插值曲線的震蕩(如3 160 m～3 260 m),而擬MQ函數(shù)則避免了這一缺陷。當(dāng)形狀參數(shù)s=0時(shí),插值曲線與原始高程曲線完全重合。當(dāng)s相對(duì)h(相鄰插值結(jié)點(diǎn)的最大間隔)逐漸增大時(shí),插值曲線在階躍點(diǎn)處的光滑程度也隨之增大,但總體形態(tài)與原曲線保持一致。因此,說明擬MQ函數(shù)具有良好的保形性和保單調(diào)性。對(duì)于形狀參數(shù)s的選擇問題,目前沒有一個(gè)量化的選擇依據(jù),但為了獲得更高的準(zhǔn)確性,s與h保持相關(guān)且趨于零時(shí)是最有利的[1]。若將擬MQ函數(shù)用于高程擬合,形狀函數(shù)不宜取得過大(s

2.2 在高密度電法數(shù)據(jù)處理中的應(yīng)用

在高密度電法勘探中,常常因電極接地電阻過大或淺地表不均勻體的影響,使得觀測數(shù)據(jù)含有一定程度的干擾噪聲(淺部因素也會(huì)影響深部數(shù)據(jù))。在布設(shè)的電極陣列中,同一根電極可能會(huì)作為供電電極或測量電極,如果某個(gè)電極接地電阻大,對(duì)于供電回路,將導(dǎo)致供電電流不穩(wěn)定對(duì)于測量回路,將產(chǎn)生讀數(shù)不穩(wěn)定,直接影響著數(shù)據(jù)的觀測精度,最終造成繪制的視電阻率擬斷面圖存在虛假或冗余異常。為壓制高密度電法數(shù)據(jù)的干擾噪聲,筆者嘗試采用擬MQ函數(shù)對(duì)高密度數(shù)據(jù)進(jìn)行降噪處理。

以邵懷高速公路某段高密度電法勘查為例,以檢驗(yàn)擬MQ函數(shù)的降噪效果。該勘查段測線長度595 m,120道電極,點(diǎn)距5 m,測量層數(shù)39層,圖3(a)為原始高密度數(shù)據(jù)繪制的擬斷面圖,可以看出在斷面圖中存在較多冗余異常。采用擬MQ函數(shù)對(duì)該數(shù)據(jù)斷面進(jìn)行逐層降噪處理,考慮高密度電法的觀測特點(diǎn),以觀測點(diǎn)距(5 m)作為形狀參數(shù),圖3(b)為處理后的高密度數(shù)據(jù)繪制的擬斷面圖,從斷面圖中可以看出,視電阻率等值線的平滑特征得到明顯改善,更有助于數(shù)據(jù)分析和解釋。

圖3 高密度視電阻率斷面處理前后對(duì)比圖Fig.3 Comparison before and after HD apparent resistivity section procession

2.3 在電測深曲線數(shù)據(jù)處理中的應(yīng)用

由于利用數(shù)值微分可以提取地球物理原始數(shù)據(jù)曲線中的弱異常信息的特點(diǎn),對(duì)新疆清河縣地下水的分布情況研究分析。以其中一條電測深曲線為例,利用擬MQ函數(shù)對(duì)其進(jìn)行計(jì)算一維加密插值(Fitting curve),進(jìn)而獲得更可靠的一階導(dǎo)數(shù)(FD)和二階導(dǎo)數(shù)(SD),具體如圖4 所示。

圖4 基于擬MQ函數(shù)的電測深數(shù)值微分曲線Fig.4 Numerical differential curve of electrical sounding based on quasi-MQ function

可以看出:在極距AB/2為90 m和125 m的位置存在微弱的低阻異常,通過對(duì)電測深曲線計(jì)算一維插值、一階導(dǎo)數(shù)和二階導(dǎo)數(shù),使這兩個(gè)位置的弱異常幅度得到增大。結(jié)合該區(qū)地質(zhì)資料及鉆探施工條件,在此處開展了鉆孔工作。發(fā)現(xiàn)在深度20 m以內(nèi)為第四系,20 m～75 m以內(nèi)為凝灰質(zhì)砂巖,75 m～130 m以內(nèi)為花崗巖,在深度130 m處打到構(gòu)造破碎帶,含水量較大,由于水壓過大導(dǎo)致沖擊鉆探方法難于施工,并于135 m深度處終孔,未打穿斷層破碎帶,終孔孔徑110 mm。采用32 m3/h水泵進(jìn)行抽水試驗(yàn),靜水位2.1 m,動(dòng)水位28 m,水位降深25.9 m,經(jīng)計(jì)算該鉆孔出水量為46.3 m3/h(1 111.2 m3/d)。因此,利用擬MQ函數(shù)可以較好提取物探原始數(shù)據(jù)曲線中的弱異常信息。

3 結(jié)論

通過擬MQ函數(shù)在物探數(shù)據(jù)處理中的應(yīng)用研究,得到以下結(jié)論:

1)利用擬MQ函數(shù)擬合地表高程,整體曲線較光滑,能很好地逼近地表形態(tài),并隨著形狀參數(shù)的增大光滑程度加大,仍能保留原始數(shù)據(jù)的整體特征,驗(yàn)證了擬MQ函數(shù)具有良好的保形性及保單調(diào)性。

2)利用擬MQ函數(shù)對(duì)高密度數(shù)據(jù)進(jìn)行插值平滑處理,能達(dá)到很好的降噪效果,從而有利于后續(xù)數(shù)據(jù)的分析與解釋。

3)利用擬MQ函數(shù)對(duì)電測深曲線進(jìn)行插值、一階求導(dǎo)、二階求導(dǎo)處理,能明顯放大地下弱異常的幅度,有利于提取物探原始數(shù)據(jù)曲線中的弱異常信息。

隨著理論研究和工程實(shí)踐的共同進(jìn)步,擬MQ函數(shù)有望在地球物理數(shù)據(jù)處理中發(fā)揮更大作用,顯示出更多優(yōu)勢,成為一種常規(guī)數(shù)學(xué)工具。