郭 鵬, 劉海飛, 劉聲凱, 柳建新, 張一凡, 李昶萱

(1.中南大學(xué) a.地球科學(xué)與信息物理學(xué)院; b.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)沙 410083;2.湖南省水文地質(zhì)環(huán)境地質(zhì)調(diào)查監(jiān)測(cè)所, 長(zhǎng)沙 410129)

0 引 言

高密度電法具有數(shù)據(jù)采集高效、 涵蓋信息量大等優(yōu)點(diǎn), 已被廣泛應(yīng)用于水文、 工程和環(huán)境等領(lǐng)域, 并且在地下水資源、 巖溶、 滑坡、 采空區(qū)、 隱伏構(gòu)造及污染場(chǎng)地監(jiān)測(cè)等方面均得到了較好的應(yīng)用[1-8]。二維反演作為高密度電法數(shù)據(jù)處理的重要技術(shù)手段, 在很大程度上影響了高密度電法的勘探效果。對(duì)高密度電阻率數(shù)據(jù)反演之前需要設(shè)定反演參數(shù), 不同的反演參數(shù)對(duì)應(yīng)得到的反演結(jié)果常常存在差別, 因此如何根據(jù)觀測(cè)裝置、 數(shù)據(jù)質(zhì)量等因素選取合理的反演參數(shù)顯得尤為重要。目前針對(duì)反演參數(shù)的選取問(wèn)題已有相關(guān)研究, 黃真萍等[9]通過(guò)數(shù)值試驗(yàn)研究了高密度反演的最佳深度轉(zhuǎn)化因子; 柳建新等[10]討論了正則化參數(shù)對(duì)反演結(jié)果分辨率和反演穩(wěn)定性的影響; 余長(zhǎng)恒等[11]討論了反演過(guò)程中阻尼系數(shù)和正演網(wǎng)格大小對(duì)突出異常體的作用。

為全面分析反演參數(shù)對(duì)反演結(jié)果的影響, 結(jié)合文獻(xiàn)[12-13], 本課題組編寫(xiě)了高密度電阻率二維反演程序, 對(duì)高密度電阻率二維反演參數(shù)的選取問(wèn)題展開(kāi)試驗(yàn)研究。通過(guò)設(shè)置不同的反演參數(shù), 包括正反演網(wǎng)格大小、 先驗(yàn)約束條件、 雅可比矩陣計(jì)算方法、 目標(biāo)函數(shù)的范數(shù)、 反演迭代次數(shù)、 初始阻尼因子、 初始模型等, 對(duì)比分析不同參數(shù)對(duì)電阻率二維反演結(jié)果的影響, 并針對(duì)裝置和數(shù)據(jù)特點(diǎn)總結(jié)出反演參數(shù)的優(yōu)化選擇, 可以為高密度電阻率二維反演處理提供參考。

1 高密度電阻率法二維正反演基本原理

1.1 高密度電阻率二維正演

高密度電法的二維正演問(wèn)題歸結(jié)為求解波數(shù)域點(diǎn)源穩(wěn)定電流場(chǎng)的邊值問(wèn)題

(1)

其中:V為波數(shù)域電位;λ為波數(shù);σ為電導(dǎo)率; 點(diǎn)源項(xiàng)f=-Iδ(A)/2;C=K1(λR)cos(R,n)/K0(λR),K0和K1分別為第二類(lèi)0階和1階修正貝塞爾函數(shù)。

利用泛函分析將波數(shù)域點(diǎn)源穩(wěn)定電流場(chǎng)的邊值問(wèn)題(式(1))轉(zhuǎn)化為等價(jià)的變分問(wèn)題

(2)

采用有限元法求解變分方程(2)[14], 可將變分問(wèn)題轉(zhuǎn)化為線性方程組

KV=f,

(3)

其中：K為對(duì)稱(chēng)稀疏數(shù)矩陣;V為網(wǎng)格節(jié)點(diǎn)處波數(shù)域電位構(gòu)成的列向量。

利用一維變帶寬存儲(chǔ)矩陣, 并用喬里斯基分解法求解方程, 即可得到波數(shù)域電位V, 然后通過(guò)傅氏逆變換將波數(shù)域電位轉(zhuǎn)換到空間域電位

(4)

其中:Wi是與波數(shù)λi對(duì)應(yīng)的傅氏逆變換權(quán)系數(shù);k為波數(shù)個(gè)數(shù)。最后根據(jù)觀測(cè)裝置將電位U轉(zhuǎn)化為視電阻率。

1.2 高密度電阻率二維反演

對(duì)數(shù)據(jù)空間和模型空間在混合范數(shù)下構(gòu)建目標(biāo)函數(shù)[15-16]

(5)

其中:da為實(shí)測(cè)數(shù)據(jù);dc(m)為模擬數(shù)據(jù);m為地電模型;mb為背景模型;C為光滑度矩陣;λ為阻尼因子;p和q為范數(shù)階, 通常取為1或2。

式(5)兩端對(duì)m求導(dǎo)并令其為0, 得

(6)

然后, 將dc(m)在初始模型m0處泰勒級(jí)數(shù)展開(kāi), 并略去二次以上的高次項(xiàng), 得

=dc(m0)+A(m-m0),

(7)

其中,A為偏導(dǎo)數(shù)矩陣。將方程(7)代入方程(6), 得到混合范數(shù)下的線性反演方程

λ·CTRmC(m-mb),

(8)

其中： Δd=da-dc(m0)為數(shù)據(jù)殘差向量; Δm=m-m0為模型修正量;Rd=|Δd-AΔm|p-2為數(shù)據(jù)加權(quán)矩陣;Rm=|C(m-mb)|q-2為模型加權(quán)矩陣。

利用共軛梯度法求解線性方程組(8), 即可得到模型參數(shù)修正量Δm, 再將其代入公式

m(k+1)=m(k)+Δm,

(9)

便得到新的地電模型, 重復(fù)上述迭代過(guò)程, 直到滿足反演終止條件為止。

2 反演參數(shù)的設(shè)計(jì)方法

以下分別對(duì)正反演網(wǎng)格剖分、 目標(biāo)函數(shù)的范數(shù)、 初始模型、 偏導(dǎo)數(shù)矩陣的結(jié)合方式、 反演迭代次數(shù)、 先驗(yàn)約束方式、 阻尼因子等參數(shù)的設(shè)計(jì)方法進(jìn)行分析討論。

2.1 正反演網(wǎng)格剖分方法

在高密度電法反演時(shí), 研究區(qū)域的網(wǎng)格剖分將涉及橫向網(wǎng)格剖分、 縱向網(wǎng)格剖分、 反演深度, 其中橫向網(wǎng)格剖分通過(guò)橫向單位電極距剖分的網(wǎng)格數(shù)調(diào)整, 縱向網(wǎng)格剖分由縱向首層網(wǎng)格厚度、 縱向?qū)雍耠S深度增加的遞增系數(shù)調(diào)整, 具體為:

橫向正演網(wǎng)格:ehs=a/c,

其中:c=1, 2, 3, 4;a為相鄰電極間距。

橫向反演網(wǎng)格:ehm=ehs×h,h=1,2,3,4。

縱向正演和反演網(wǎng)格: 縱向正、 反演網(wǎng)格單元厚度一致, 網(wǎng)格尺度采用逐漸遞增方式

evs(i)=evm(i)=evs(i-1)×f,i=1,2,…,n。

(10)

其中:f∈[1.1, 2]為層厚度調(diào)整系數(shù);evs(0)=a×g為縱向首層網(wǎng)格厚度;g∈[0.1, 0.5]為首層網(wǎng)格調(diào)整系數(shù)。

反演深度:D=k×ABmax/2,

其中:ABmax/2為最大視深度;k∈[0.25, 0.8]為深度調(diào)整系數(shù)。

2.2 目標(biāo)函數(shù)的范數(shù)

根據(jù)線性反演方程(8), 可以采用4種方式構(gòu)建目標(biāo)函數(shù), 即： ①數(shù)據(jù)空間p和模型空間q均為2范數(shù); ②數(shù)據(jù)空間p為2范數(shù), 模型空間q為1范數(shù); ③數(shù)據(jù)空間p為1范數(shù), 模型空間q為2范數(shù); ④數(shù)據(jù)空間p和模型空間q均為1范數(shù)。

圖1 正演網(wǎng)格示意圖

2.3 初始模型

由于廣義線性反演結(jié)果很大程度上受初始模型影響, 在高密度電阻率二維反演時(shí), 可以給定均勻或非均勻初始模型(標(biāo)識(shí)IM=0為均勻,IM=1非均勻)。

對(duì)于均勻模型, 可將實(shí)測(cè)數(shù)據(jù)取平均值作為初始模型, 即

(11)

其中：ρa(bǔ)i為實(shí)測(cè)視電阻率,i=1, 2, …,n；n為數(shù)據(jù)個(gè)數(shù)。

對(duì)于非均勻初始模型, 可利用反距離加權(quán)插值公式

(12)

根據(jù)實(shí)測(cè)數(shù)據(jù)的視記錄點(diǎn)對(duì)地下模型進(jìn)行二維插值, 其中di為已知點(diǎn)到待插點(diǎn)的距離,m為選取的臨近數(shù)據(jù)點(diǎn)的個(gè)數(shù)。

2.4 阻尼因子

誤差收斂曲線通常在前幾次迭代下降較快而后緩慢下降, 考慮采用函數(shù)

λ(k)=a·k-2+b

(13)

構(gòu)造阻尼λ序列。其中:k為迭代序號(hào);a和b為待求系數(shù)。首先給定反演迭代次數(shù)nmax, 并根據(jù)數(shù)據(jù)所含噪音情況選擇相對(duì)合理的初始阻尼λmax(通常在0.1～1內(nèi)選擇), 為方便起見(jiàn), 設(shè)最小阻尼為λmin=λmax/10, 這樣即可確定a、b。若取nmax=5,λmax=0.5,λmin=0.05, 則根據(jù)式(13)可得λ序列為0.5、 0.152 9、 0.089、 0.066、 0.056、 0.05。

2.5 反演迭代次數(shù)

選取反演迭代次數(shù)nmax(如5～15)作為反演的終止條件, 若僅以擬合差作為終止條件， 實(shí)際中有時(shí)難以滿足。

2.6 先驗(yàn)約束

為減少反演產(chǎn)生的假異常, 反演時(shí)需要對(duì)模型空間施加先驗(yàn)約束, 主要包括光滑約束SC(選擇0或1)和背景約束B(niǎo)C(選擇0或1)。

2.7 偏導(dǎo)數(shù)矩陣的結(jié)合方式

在高密度電阻率反演中, 常采用以下兩種計(jì)算偏導(dǎo)數(shù)矩陣的方法。

①互換原理: 視電阻率對(duì)模擬電阻率的偏導(dǎo)數(shù)計(jì)算可以歸結(jié)為網(wǎng)格節(jié)點(diǎn)電位對(duì)模擬電阻率的偏導(dǎo)數(shù)計(jì)算。互換原理表明, 在某節(jié)點(diǎn)iA供電時(shí)jM點(diǎn)的電位等價(jià)于在jM點(diǎn)供電時(shí)iA點(diǎn)的電位, ?V/?ρ為所有網(wǎng)格節(jié)點(diǎn)的電位的線性組合, 經(jīng)過(guò)變換求出偏導(dǎo)數(shù)矩陣。

②擬牛頓法: 采用Broyden秩一校正公式近似計(jì)算雅可比矩陣

(14)

其中：Bk為第k次迭代的雅可比矩陣;sk為第k次迭代模型參數(shù)的改正量;yk為第k次和第k+1次迭代模擬視電阻率的差。

擬牛頓法計(jì)算偏導(dǎo)數(shù)矩陣的效率遠(yuǎn)遠(yuǎn)高于互換原理法, 本文將兩者結(jié)合起來(lái), 即前nc次迭代采用互換原理, 后nb=nmax-nc次迭代采用擬牛頓法, 在保證不降低分辨率的情況下提高反演速度。

3 反演試驗(yàn)結(jié)果分析

3.1 反演試驗(yàn)斷面介紹

選取兩種高密度電阻率法觀測(cè)裝置的實(shí)測(cè)斷面作為試驗(yàn)斷面。試驗(yàn)斷面一為浙江某地溫納Beta裝置采集的高密度電阻率數(shù)據(jù)(圖2a), 點(diǎn)距2 m, 電極60個(gè), 測(cè)點(diǎn)570個(gè), 最大視深度為56 m, 測(cè)線地形起伏。根據(jù)已知地質(zhì)信息, 該地區(qū)的地質(zhì)情況為: 里程0～50 m主要以砂頁(yè)巖為主, 里程75～85 m主要以灰?guī)r為主, 里程90～100 m、 深10～20 m處存在溶洞。

試驗(yàn)斷面二為湖南某地溫納施倫貝爾裝置采集的高密度電阻率數(shù)據(jù)(圖2b), 點(diǎn)距5 m, 電極60個(gè), 測(cè)點(diǎn)722個(gè), 測(cè)線地形水平。已有鉆孔資料顯示: 在730 m處0～9.8 m為黏土, 9.8～15.2 m為白云質(zhì)灰?guī)r, 15.2～18.6 m為含礫長(zhǎng)黏土, 18.6～72.8 m為白云質(zhì)灰?guī)r, 其中67.3～73.6 m為溶洞。

3.2 不同反演參數(shù)對(duì)反演結(jié)果的影響分析

3.2.1 正反演網(wǎng)格 采用試驗(yàn)斷面一作為測(cè)試對(duì)象, 反演時(shí)只改變縱向網(wǎng)格剖分參數(shù), 其他反演參數(shù)如表1所示。通過(guò)增大層厚度調(diào)整系數(shù)f, 正演網(wǎng)格尺寸增大, 深部的分辨率略有降低, 但不影響異常位置, 網(wǎng)格越小均方差越小且收斂速度快, 耗時(shí)長(zhǎng)(圖3a—c); 通過(guò)增大深度調(diào)整系數(shù)k, 反演深度越大, 越能反映深層電阻率情況以解決邊界問(wèn)題, 但耗時(shí)長(zhǎng)(圖3d—f), 如果超出測(cè)深范圍, 易造成虛假異常(此處f、k對(duì)應(yīng)2.1節(jié)所述, 下文相同)。采用稀疏網(wǎng)格反演耗時(shí)短, 可壓制噪聲; 精細(xì)網(wǎng)格降低圍巖對(duì)異常體的影響, 有利于分辨異常界面, 突出自身異常。

表1 設(shè)置不同正反演網(wǎng)格剖分參數(shù)時(shí)的參數(shù)設(shè)置

圖3 改變縱向網(wǎng)格調(diào)整系數(shù)時(shí)試驗(yàn)斷面一的反演結(jié)果

3.2.2 混合范數(shù) 反演時(shí)數(shù)據(jù)空間和模型空間選取不同的范數(shù)組合構(gòu)建目標(biāo)函數(shù), 其他反演參數(shù)設(shè)置如表2所示。對(duì)試驗(yàn)斷面一進(jìn)行反演試算, 反演結(jié)果如圖4所示。不同范數(shù)組合的反演結(jié)果形態(tài)類(lèi)似, 僅局部區(qū)域存在少量差異： 數(shù)據(jù)空間選擇1范數(shù)時(shí), 突出淺部高阻, 模型空間選擇1范數(shù)時(shí), 深部低阻范圍大； 當(dāng)數(shù)據(jù)空間和模型空間均選擇2范數(shù)時(shí), 限制了深部低阻范圍, 壓制了淺部高阻。范數(shù)的選擇主要與數(shù)據(jù)噪聲的分布規(guī)律, 即高斯分布或正態(tài)分布有關(guān)。

圖4 數(shù)據(jù)空間與模型空間采用不同范數(shù)組合時(shí)試驗(yàn)斷面一的反演結(jié)果

表2 改變范數(shù)組合反演時(shí)的參數(shù)設(shè)置

3.2.3 初始模型 反演時(shí)分別選取均勻和非均勻初始模型對(duì)試驗(yàn)斷面二進(jìn)行反演試算, 其他反演參數(shù)設(shè)置見(jiàn)表3。具體反演結(jié)果如圖5所示, 反演結(jié)果對(duì)初始模型的依賴(lài)性較大, 均勻初始模型的反演結(jié)果不能反映出溶洞的存在(圖5a), 而非均勻初始模型的反演結(jié)果能較好地反映異常體的存在。對(duì)于擬斷面圖中深部顯示的弱異常信息, 通過(guò)給定較優(yōu)的初始模型可以提高反演結(jié)果的分辨率。

圖5 采用不同初始模型時(shí)試驗(yàn)斷面二的反演結(jié)果

表3 給定不同初始模型反演時(shí)的參數(shù)設(shè)置

3.2.4 阻尼因子 對(duì)試驗(yàn)斷面二選取不同的阻尼因子進(jìn)行反演試算, 其他反演參數(shù)設(shè)置見(jiàn)表4。初始阻尼因子分別選取0.1和0.5時(shí)反演結(jié)果分別如圖6a、 6b所示, 反演結(jié)果受阻尼因子的影響較大： 初始阻尼因子為0.1時(shí)， 深部溶洞異常在反演剖面中得到較好的體現(xiàn)； 而初始阻尼因子為0.5時(shí)， 深部溶洞異常已被平滑掉。實(shí)際中由于數(shù)據(jù)所含噪聲通常不同, 很難根據(jù)噪聲水平選取能使模型得到最佳分辨又可使誤差收斂曲線穩(wěn)步下降的最佳阻尼因子, 但可以通過(guò)給定不同的阻尼因子進(jìn)行反演試驗(yàn), 從中選取符合地質(zhì)規(guī)律的反演結(jié)果。

圖6 設(shè)置不同阻尼因子時(shí)試驗(yàn)斷面二的反演結(jié)果

3.2.5 反演迭代次數(shù) 對(duì)試驗(yàn)斷面二選取不同的迭代次數(shù)進(jìn)行反演試算, 其他反演參數(shù)設(shè)置見(jiàn)表5。迭代4次和8次的反演結(jié)果分別如圖7a、 7b所示, 當(dāng)反演迭代4次時(shí)(耗時(shí)60 s), 溶洞未能反映出來(lái), 當(dāng)?shù)?次時(shí)(耗時(shí)87 s), 溶洞已得到較好的分辨。因此, 反演迭代次數(shù)也是影響反演結(jié)果的一個(gè)因素, 迭代次數(shù)不足, 反演結(jié)果可能未達(dá)到最佳逼近, 迭代次數(shù)過(guò)多, 會(huì)增加反演耗時(shí), 實(shí)際中反演迭代6～10次誤差收斂曲線已基本趨于穩(wěn)定。

表5 改變反演迭代次數(shù)時(shí)的參數(shù)設(shè)置

圖7 設(shè)置不同迭代次數(shù)時(shí)試驗(yàn)斷面二的反演結(jié)果

3.2.6 先驗(yàn)約束 對(duì)試驗(yàn)斷面二選取不同的先驗(yàn)約束進(jìn)行反演試算, 其他反演參數(shù)設(shè)置見(jiàn)表6。施加背景約束和同時(shí)施加光滑與背景約束的反演結(jié)果分別如圖8a、 8b所示, 兩者反演異常形態(tài)大致類(lèi)似, 僅施加背景約束時(shí)電阻率等值線多呈鋸齒狀, 而同時(shí)施加光滑與背景約束時(shí)電阻率等值線較光滑。實(shí)際中選取背景約束更能凸顯弱異常信息, 而對(duì)于強(qiáng)異常信息同時(shí)施加背景和光滑約束的反演異常等值線將更為平滑流暢。

表6 改變先驗(yàn)約束反演時(shí)的參數(shù)設(shè)置

圖8 采用不同先驗(yàn)約束時(shí)試驗(yàn)斷面二的反演結(jié)果

3.2.7 互換原理與擬牛頓法結(jié)合計(jì)算雅可比矩陣 采用互換原理和擬牛頓法相結(jié)合的方式計(jì)算偏導(dǎo)數(shù)矩陣, 當(dāng)反演迭代次數(shù)設(shè)置為6次時(shí), 改變互換原理和擬牛頓法使用的次數(shù), 其他反演參數(shù)設(shè)置見(jiàn)表7, 對(duì)試驗(yàn)斷面一進(jìn)行反演試算, 結(jié)果分別如圖9a、 9b所示。隨著互換原理計(jì)算雅可比矩陣次數(shù)的增加, 斷面異常形態(tài)無(wú)明顯變化, 互換原理使用3次時(shí)反演耗時(shí)72 s, 互換原理使用5次時(shí)反演耗時(shí)132 s。這說(shuō)明在電阻率二維反演中, 互換原理計(jì)算偏導(dǎo)數(shù)矩陣比較耗時(shí), 而擬牛頓法則比較節(jié)省時(shí)間, 嘗試將互換原理與擬牛頓法結(jié)合計(jì)算偏導(dǎo)數(shù)矩陣?？紤]到反演迭代誤差收斂曲線通常僅前2～3次收斂較快, 在后續(xù)迭代中收斂曲線下降緩慢。因此, 可將前2～3次迭代采用互換原理法后續(xù)迭代采用擬牛頓法計(jì)算偏導(dǎo)數(shù)矩陣, 在保證反演穩(wěn)定和分辨率的情況下加快反演計(jì)算速度。

表7 采用不同方式計(jì)算偏導(dǎo)數(shù)矩陣時(shí)的反演參數(shù)設(shè)置

圖9 互換原理與擬牛頓法結(jié)合時(shí)試驗(yàn)斷面一的反演結(jié)果

4 結(jié) 論

(1)通過(guò)對(duì)不同反演參數(shù)進(jìn)行反演試驗(yàn), 正反演網(wǎng)格大小應(yīng)根據(jù)裝置特點(diǎn)選擇, 偶極裝置橫向分辨率高, 可設(shè)置橫向正反演網(wǎng)格為點(diǎn)距的一半或更小; 溫施裝置探測(cè)深度大, 設(shè)置縱向反演網(wǎng)格深度盡量大于最大視深度的一半。通常情況下, 橫向正反演單元設(shè)置為點(diǎn)距的一半, 縱向首層網(wǎng)格單元設(shè)置為點(diǎn)距的0.2倍, 縱向網(wǎng)格調(diào)整系數(shù)設(shè)置為1.2, 縱向反演網(wǎng)格深度調(diào)整系數(shù)在很大程度上取決于最大視深度AB/2和地層電性, 電性值越低深度調(diào)整系數(shù)越小, 反之越大, 一般情況下設(shè)置為0.5, 必要時(shí)加大。

(2)初始模型和阻尼因子在很大程度上影響反演結(jié)果的分辨率, 如果異常信息較弱, 則選擇非均勻初始模型和較小的初始阻尼因子。

(3)當(dāng)數(shù)據(jù)噪聲較小且服從拉普拉斯分布, 數(shù)據(jù)空間和模型空間均選擇2范數(shù), 反之當(dāng)數(shù)據(jù)噪聲較大且服從高斯分布時(shí), 數(shù)據(jù)空間或模型空間選擇L1范數(shù); 對(duì)于其他反演參數(shù), 通常設(shè)置為同時(shí)施加光滑和背景約束、 反演迭代次數(shù)為6(前兩次迭代采用互換原理、 后4次迭代采用擬牛頓法計(jì)算偏導(dǎo)數(shù)矩陣)。

(4)在處理和解釋實(shí)測(cè)高密度電法數(shù)據(jù)時(shí), 選擇合適反演參數(shù)能夠提高反演效率和反演結(jié)果的準(zhǔn)確性, 討論高密度電法二維反演參數(shù)是非常有必要的。實(shí)際中對(duì)于較為復(fù)雜的地質(zhì)問(wèn)題, 需要反復(fù)調(diào)整各參數(shù)進(jìn)行反演試驗(yàn), 從中選取符合地質(zhì)規(guī)律的反演結(jié)果。