蘇鵬，楊進(jìn)

(中國(guó)地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083)

0 引言

電阻率法廣泛應(yīng)用于水文、工程、環(huán)境地質(zhì)問(wèn)題的探測(cè)工程中。通常，電阻率法用于對(duì)靜態(tài)的地下電性結(jié)構(gòu)進(jìn)行成像，然而在環(huán)境監(jiān)測(cè)、地質(zhì)災(zāi)害監(jiān)測(cè)及地下水污染監(jiān)測(cè)中，由于地下介質(zhì)處于動(dòng)態(tài)變化的狀態(tài)，如沿海地區(qū)地下水的海水入侵、山區(qū)滑坡巖移等，使用動(dòng)態(tài)監(jiān)測(cè)的方法更為有效。針對(duì)上述這些動(dòng)態(tài)問(wèn)題，需要通過(guò)不同時(shí)間的多次觀測(cè)數(shù)據(jù)，利用時(shí)移電阻率法(time-lapse,ERT)來(lái)解決。目前，時(shí)移電阻率法正處于快速發(fā)展階段，在水文、環(huán)境領(lǐng)域內(nèi)的應(yīng)用范圍逐漸擴(kuò)大,例如：利用時(shí)移監(jiān)測(cè)方法研究滲流運(yùn)移情況[1-3]，研究水文地質(zhì)參數(shù)變化規(guī)律[4-5]、溶質(zhì)運(yùn)移[6-7]及研究地下熱傳導(dǎo)情況[8-10]；在國(guó)內(nèi)，李飛等[11]利用時(shí)移高密度法研究地下煤層開(kāi)采所引起的上覆巖層變形破壞，孫大利等[12]使用高密度時(shí)移監(jiān)測(cè)分析堤壩隱患的時(shí)移特征。

時(shí)移電阻率反演算法基于常規(guī)2D/3D電阻率反演算法，但更加關(guān)注不同時(shí)刻觀測(cè)數(shù)據(jù)差異以及相應(yīng)地下電性結(jié)構(gòu)模型參數(shù)的差異情況，旨在消除多個(gè)數(shù)據(jù)中的隨機(jī)誤差，凸顯地下介質(zhì)隨時(shí)間變化而發(fā)生的真實(shí)變化。Daily等[13]反演初始數(shù)據(jù)集和后續(xù)數(shù)據(jù)集的比值以凸顯地下電性變化；LaBrecque等[14]試圖最小化初始數(shù)據(jù)集和后續(xù)數(shù)據(jù)集的數(shù)據(jù)差值及其與響應(yīng)模型差值；Kim等[15]提出將數(shù)據(jù)集和模型參數(shù)放入時(shí)空域(4D)中進(jìn)行離散，以進(jìn)行全4D反演；Hayley等[16]將多個(gè)不同時(shí)間點(diǎn)的數(shù)據(jù)集和模型參數(shù)同時(shí)反演迭代；Karaoulis等[17]對(duì)全4D時(shí)移電阻率反演方法進(jìn)行改進(jìn)，將時(shí)間正則項(xiàng)設(shè)置為可變的，并對(duì)超參數(shù)調(diào)優(yōu)方式加以改善；Loke等[18]利用光滑約束的最小二乘法結(jié)合L-curve的參數(shù)調(diào)優(yōu)方法來(lái)提高4D反演速度。

本文對(duì)高密度電法監(jiān)測(cè)數(shù)據(jù)時(shí)移反演做了初步研究，在對(duì)比分析上述不同時(shí)移算法的基礎(chǔ)上，采用Karaoulis的思路推導(dǎo)并實(shí)現(xiàn)了時(shí)移反演算法；通過(guò)模擬數(shù)據(jù)對(duì)比了常規(guī)電阻率反演和時(shí)移反演的效果，展示了時(shí)移反演算法在對(duì)動(dòng)態(tài)地下目標(biāo)體的監(jiān)測(cè)追蹤中的優(yōu)越性。

1 時(shí)移電阻率反演

1.1 電阻率法反演理論

對(duì)于常規(guī)的2D/3D電阻率法而言，反演是通過(guò)所采集數(shù)據(jù)獲得可靠地下電性結(jié)構(gòu)的過(guò)程。在電阻率反演中，地下電性結(jié)構(gòu)為離散的模型參數(shù)電阻率，其反演問(wèn)題是一個(gè)病態(tài)問(wèn)題，故需要采取正則化的方法來(lái)避免反演的病態(tài)問(wèn)題：

Φ(m)=Φd+βΦm，

(1)

式中：Φd為數(shù)據(jù)擬合差；Φm為模型正則項(xiàng)；β為正則化參數(shù)，用以平衡兩項(xiàng)的貢獻(xiàn)。上式具體展開(kāi)形式如下[19]：

(2)

式中：Wd是一個(gè)數(shù)據(jù)加權(quán)矩陣，它是一個(gè)對(duì)角矩陣；d(m)為正演所得預(yù)測(cè)數(shù)據(jù)；dobs為觀測(cè)數(shù)據(jù)；Wm為模型加權(quán)矩陣；m為迭代模型；mref為參考模型。

(3)

1.2 時(shí)移電阻率反演

相較于常規(guī)電阻率方法，時(shí)移電阻率法反演的目標(biāo)函數(shù)中添加了一個(gè)時(shí)間約束項(xiàng)，即目標(biāo)函數(shù)為：

(4)

式(4)可展開(kāi)為如下形式：

(5)

(6)

對(duì)于時(shí)間加權(quán)項(xiàng)Wt，形式如下(以3個(gè)模型為例)：

對(duì)比式(3)與式(6)，兩者不同之處在于時(shí)移反演算法中需要構(gòu)建時(shí)間加權(quán)矩陣，改加權(quán)矩陣的作用是體現(xiàn)出不同時(shí)間的模型差異對(duì)目標(biāo)函數(shù)的貢獻(xiàn)。筆者編寫(xiě)了時(shí)移電阻率反演程序，實(shí)現(xiàn)2D高密度電法時(shí)移反演。由于需要對(duì)模型約束和時(shí)間約束兩個(gè)正則化參數(shù)調(diào)優(yōu)，對(duì)模型約束項(xiàng)使用ACB方法[20]進(jìn)行自動(dòng)優(yōu)化。實(shí)際上，模型參數(shù)在時(shí)間域上的變化要小于空間域的分布，而且各不同時(shí)間點(diǎn)的基本地下結(jié)構(gòu)是保持不變的，所以時(shí)間正則項(xiàng)對(duì)目標(biāo)函數(shù)的貢獻(xiàn)要小于空間正則項(xiàng)，本文假設(shè)模型參數(shù)在時(shí)間域上的變化程度是中等的，時(shí)間正則化參數(shù)α=0.01[17]。

2 模擬數(shù)據(jù)反演分析

為了研究時(shí)移電阻率法反演的效果，首先建立一組(多個(gè))正演模型(圖1)。其中，T0模型(圖1a)由兩部分組成，位于0～19 m，埋深4.44～7.88 m區(qū)域的電阻率設(shè)置為200 Ω·m，其余區(qū)域設(shè)定為100 Ω·m；后續(xù)的T1、T2、T3在T0模型基礎(chǔ)之上，設(shè)置一個(gè)隨著時(shí)間變化的目標(biāo)低阻體，電阻率為10 Ω·m，隨著時(shí)間的推移，目標(biāo)體逐漸擴(kuò)散，在橫向和縱向均有變化，以便于檢驗(yàn)反演算法的橫向(縱向)分辨能力。模擬常用的溫納陣列類型的數(shù)據(jù)，測(cè)線共有51個(gè)電極，電極距為1 m，在模擬數(shù)據(jù)中加入了5%的隨機(jī)噪聲。

圖1 不同時(shí)間點(diǎn)的模型Fig.1 The models of different time steps

在此基礎(chǔ)上，對(duì)以上各個(gè)模擬數(shù)據(jù)分別進(jìn)行單獨(dú)反演和時(shí)移反演，經(jīng)過(guò)對(duì)所編寫(xiě)的反演程序反復(fù)實(shí)驗(yàn)發(fā)現(xiàn)，通常情況下反演迭代過(guò)程擬合差下降速度較快，一般僅需迭代6次左右即可收斂，故將最大迭代次數(shù)設(shè)置為6，兩種不同反演方式的結(jié)果如圖2所示。

圖2 模擬數(shù)據(jù)單獨(dú)和時(shí)移反演結(jié)果Fig.2 The single inversion and time-lapse inversion results on simulation data

從圖2可以看到，單獨(dú)反演與時(shí)移反演兩種方法對(duì)不同數(shù)據(jù)均有較好的約束。其中，不同數(shù)據(jù)單獨(dú)反演的擬合差在3%～4%之間，而時(shí)移反演的總擬合差則小于1%。由于噪聲和誤差的影響，反演結(jié)果中存在與模型無(wú)法對(duì)應(yīng)的假異常。對(duì)比兩種反演方法的結(jié)果，時(shí)移反演的假異常明顯比單獨(dú)反演的要少，而且對(duì)目標(biāo)異常體的約束效果較好。為了突出地下目標(biāo)體的變化情況，兩種方法均以第一個(gè)模型的反演結(jié)果ρT0為基準(zhǔn)，通過(guò)ρTi-ρT0的方式來(lái)顯示各個(gè)模型的變動(dòng)情況，結(jié)果如圖3所示。

圖3 不同模型反演結(jié)果的差異Fig.3 The difference images of inversion results

圖3中，代表了T1、T2、T3相對(duì)與T0時(shí)的地下目標(biāo)變化在反演結(jié)果中的相應(yīng)體現(xiàn)，由于所設(shè)定的變化目標(biāo)為一低阻體，故隨著該目標(biāo)的體積不斷增大，圖3中數(shù)值-0.3～-0.5所占的面積不斷增大，說(shuō)明兩種方法均可以追蹤地下電阻率動(dòng)態(tài)變化的目標(biāo)體。但是，對(duì)比時(shí)移反演(圖3b)和單獨(dú)反演(圖3a)的結(jié)果可以看出，時(shí)移反演對(duì)其變化的目標(biāo)具有更好的效果，而單獨(dú)反演的效果較差，其對(duì)變化的目標(biāo)體的約束不夠，并且依然存在不同程度的假異常。

3 結(jié)論

1) 通過(guò)對(duì)反演算法的理論推導(dǎo)可發(fā)現(xiàn)，時(shí)移反演算法有別于常規(guī)算法之處在于該算法在目標(biāo)函數(shù)中添加了一個(gè)時(shí)間約束項(xiàng)，所以時(shí)移反演算法的效果主要取決于時(shí)間約束項(xiàng)構(gòu)建得是否合理。

2) 監(jiān)測(cè)的多個(gè)數(shù)據(jù)集中存在不同程度的隨機(jī)噪聲，是導(dǎo)致反演中出現(xiàn)假異常的原因之一。由于時(shí)移反演算法中時(shí)移約束項(xiàng)的存在，將不同時(shí)間點(diǎn)的數(shù)據(jù)集統(tǒng)一放入一個(gè)目標(biāo)函數(shù)中進(jìn)行綜合考量，從而使得時(shí)間序列上相鄰的響應(yīng)模型構(gòu)成一定的聯(lián)系，在時(shí)空域(4D)尺度上進(jìn)行最優(yōu)化，而不是孤立地對(duì)各個(gè)數(shù)據(jù)集和模型進(jìn)行單獨(dú)反演。

3) 模擬數(shù)據(jù)的實(shí)驗(yàn)表明，對(duì)于隨時(shí)間變化的電阻率異常體，電阻率法時(shí)移反演克服了多個(gè)數(shù)據(jù)集存在不同隨機(jī)噪聲的影響，能夠較好地約束地下動(dòng)態(tài)目標(biāo)，追蹤其隨時(shí)間的變化情況。