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

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

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

0 引 言

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

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

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

1.1 高密度電阻率二維正演

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

(1)

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

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

(2)

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

KV=f,

(3)

其中：K為對稱稀疏數矩陣;V為網格節點處波數域電位構成的列向量。

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

(4)

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

1.2 高密度電阻率二維反演

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

(5)

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

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

(6)

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

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

(7)

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

λ·CTRmC(m-mb),

(8)

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

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

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

(9)

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

2 反演參數的設計方法

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

2.1 正反演網格剖分方法

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

橫向正演網格:ehs=a/c,

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

橫向反演網格:ehm=ehs×h,h=1,2,3,4。

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

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

(10)

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

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

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

2.2 目標函數的范數

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

圖1 正演網格示意圖

2.3 初始模型

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

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

(11)

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

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

(12)

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

2.4 阻尼因子

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

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

(13)

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

2.5 反演迭代次數

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

2.6 先驗約束

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

2.7 偏導數矩陣的結合方式

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

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

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

(14)

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

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

3 反演試驗結果分析

3.1 反演試驗斷面介紹

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

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

3.2 不同反演參數對反演結果的影響分析

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

表1 設置不同正反演網格剖分參數時的參數設置

圖3 改變縱向網格調整系數時試驗斷面一的反演結果

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

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

表2 改變范數組合反演時的參數設置

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

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

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

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

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

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

表5 改變反演迭代次數時的參數設置

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

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

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

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

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

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

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

4 結 論

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

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

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

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