不同反演參數(shù)對高密度電阻率二維反演結果的影響

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

(1.中南大學 a.地球科學與信息物理學院; b.有色資源與地質災害探查湖南省重點實驗室, 長沙 410083;2.湖南省水文地質環(huán)境地質調查監(jiān)測所, 長沙 410129)

0 引 言

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

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

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

1.1 高密度電阻率二維正演

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

(1)

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

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

(2)

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

KV=f,

(3)

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

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

(4)

其中:Wi是與波數(shù)λi對應的傅氏逆變換權系數(shù);k為波數(shù)個數(shù)。最后根據(jù)觀測裝置將電位U轉化為視電阻率。

1.2 高密度電阻率二維反演

對數(shù)據(jù)空間和模型空間在混合范數(shù)下構建目標函數(shù)[15-16]

(5)

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

式(5)兩端對m求導并令其為0, 得

(6)

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

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

(7)

其中,A為偏導數(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ù)加權矩陣;Rm=|C(m-mb)|q-2為模型加權矩陣。

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

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

(9)

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

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

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

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

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

橫向正演網(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]為層厚度調整系數(shù);evs(0)=a×g為縱向首層網(wǎng)格厚度;g∈[0.1, 0.5]為首層網(wǎng)格調整系數(shù)。

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

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

2.2 目標函數(shù)的范數(shù)

根據(jù)線性反演方程(8), 可以采用4種方式構建目標函數(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 初始模型

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

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

(11)

其中：ρai為實測視電阻率,i=1, 2, …,n；n為數(shù)據(jù)個數(shù)。

對于非均勻初始模型, 可利用反距離加權插值公式

(12)

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

2.4 阻尼因子

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

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

(13)

構造阻尼λ序列。其中:k為迭代序號;a和b為待求系數(shù)。首先給定反演迭代次數(shù)nmax, 并根據(jù)數(shù)據(jù)所含噪音情況選擇相對合理的初始阻尼λmax(通常在0.1～1內選擇), 為方便起見, 設最小阻尼為λ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)作為反演的終止條件, 若僅以擬合差作為終止條件， 實際中有時難以滿足。

2.6 先驗約束

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

2.7 偏導數(shù)矩陣的結合方式

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

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

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

(14)

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

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

3 反演試驗結果分析

3.1 反演試驗斷面介紹

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

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

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

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

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

圖3 改變縱向網(wǎng)格調整系數(shù)時試驗斷面一的反演結果

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

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

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

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

圖5 采用不同初始模型時試驗斷面二的反演結果

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

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

圖6 設置不同阻尼因子時試驗斷面二的反演結果

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

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

圖7 設置不同迭代次數(shù)時試驗斷面二的反演結果

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

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

圖8 采用不同先驗約束時試驗斷面二的反演結果

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

表7 采用不同方式計算偏導數(shù)矩陣時的反演參數(shù)設置

圖9 互換原理與擬牛頓法結合時試驗斷面一的反演結果

4 結 論

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

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

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

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