多先驗信息約束的三維電阻率反演方法

郭來功 戴廣龍 楊本才 薛俊華 陳本良

(①安徽理工大學電信學院，安徽淮南 232001； ②安徽理工大學能源學院，安徽淮南 232001；③深部煤炭開采與環境保護國家重點實驗室，安徽淮南 232001)

1 引言

直流電法是一種重要的地球物理勘探方法，在礦井突水、地質滑坡、隧道施工等巖土工程領域，地下水勘探、環境地質監測等環境工程領域，以及礦產資源調查等資源勘探領域均有廣泛應用[1-3]。直流電法三維成像存在固有的多解性問題，會出現異常體定位不準確的現象。另外，反演數據量大、求解效率不高，其在時效性要求較高的應用中易受到限制，因此提高反演精度和計算效率是直流電三維勘探需要解決的主要問題[4,5]。

在三維反演中引入更多的先驗信息進行約束，是提高成像效果的主要手段之一。劉征宇等[6]提出幾種有效方法，一是引入距離加權約束，構建距離加權因子分布模型，對于不同的電源點距離，施加指數規律的權重因子，對不同單元賦予差異權重，反演的異常體形態得到有效改善；二是引入不等式約束，將介質電阻率的取值范圍作為先驗約束條件，提高反演精度[7]；三是引入參考模型約束，將已經獲取的異常體結構信息作為先驗信息施加到約束條件中，有效去除反演成像中的假異常[8]。劉鑫明等[9]采用光滑約束構造三維反演函數中的粗糙度矩陣，采用合理的網格剖分，在注漿檢測中取得良好的效果。吳小平等[10]采用不完全Gauss-Newton法，應用非結構網格實現任意偶極—偶極視電阻率的三維反演，較好地展現了地下真實電阻率結構。Johnson等[11]在正演建模時通過埋設導電件，建立精準模型，提高反演成像效果。Johnson 等[12]還采用并行分布式計算方法，提高了反演效率。Farquharson[13]針對最小構造反演，提出一種分段常數的地質模型構造方法，采用L1范數準則下的迭代重加權最小二乘方法，成像結果與實際地質結構高度接近。在反演方法上，基于全局優化的模擬退火方法、人工神經網絡方法等受到計算機計算能力的限制，目前還無法全面應用到實踐，當前的工程應用以基于最小二乘準則的線性反演為主[14-18]。

實際應用中，仍有諸多問題需要解決。關于反演算法，當目標函數確定后，先驗信息無法靈活地施加到約束中，例如引入距離加權約束后，無法再施加不等式約束信息；另外，反演耗時過長，普通計算機無法實現快速反演。本文采用分布式并行計算方法[19]，在最小二乘準則的基礎上，在重加權函數中增加結構度量、權重均值和標準偏差等參數，對模型的目標區域施加不同約束，并通過合成數據進行直流電三維反演算例分析。

2 目標函數的建立

三維有限元模型空間中，模型參數用向量m表示，模型響應為f，測量數據用列向量d表示。反演過程中，網格固定不變，其目的是確定模型參數。反演通過構造一個包含模型參數和測量數據的目標函數，使其最小化。根據Farquharson[13]提出的一般準則，目標函數Φ可表示為

Φ=Φd(ud)+βΦm(um)

(1)

式中：下標“m”和“d”分別代表模型和數據；Φd(ud)是度量模型響應數據與測量數據d之間擬合差的函數；Φm(um)是度量模型參數變化的函數；β是正則化參數，控制擬合數據和模型向量變化對目標函數的影響。其中ud的表達式為

ud=Wd(d-f)

(2)

式中Wd是數據加權對角矩陣，其元素是測量數據中誤差(噪聲)的標準偏差σi的倒數，即1/σi。如果測量時沒有考慮誤差，則需對誤差做合理的估計[20]。假設誤差由測量電壓誤差ΔU和一個固定百分比誤差(q%)構成，則第i個數據的標準偏差為

(3)

式中：Ui為第i個測量電壓； ΔUi為第i個數據的測量電壓誤差。目標函數Φ中模型參數向量um可寫作

um=Wm(m-mref)

(4)

式中：Wm是模型參數權重矩陣；mref表示參考模型。在一種簡單的情況下，考慮三維空間的各向異性，有

(5)

式中：Cx、Cy和Cz分別是x、y、z方向的平滑矩陣；αx、αy和αz分別是x、y、z方向的相對權重。平滑矩陣C的一般形式是一階差分矩陣[21]。

3 基于迭代重加權最小二乘法的三維反演

采用最小二乘法時，反演方程的求解采用將非線性問題線性化的方法，經過多次迭代修正模型參數，使估計解不斷逼近真實解。迭代重加權的核心思想是通過計算一系列權系數并不斷更新，得到目標函數的最優解。假設在第n次迭代時，得到最佳參數模型m(n)，使目標函數最小化。由式(1)，第n次迭代的目標函數可以表示為

(6)

設模型參數向量的變化量為Δm=m(n)-m(n-1)，則第n次迭代時，得到模型m(n-1)的正演結果

f(n)≈f(n-1)+JΔm

(7)

式中J為偏導數的靈敏度雅克比矩陣，其元素為

(8)

表示第j個模型參數的變化導致的第i個模型響應的變化，因此，由式(2)和式(4)可得

ud=Wd[d-f(n-1)-JΔm]

(9)

um=Wm[m(n-1)+Δm-mref]

(10)

式(6)是關于Δm的函數，最小化式(6)，令其偏導為零，結合式(9)、式(10)，應用Gauss-Newton法求解，可得

(11)

式中：R是對角矩陣。Rd的元素為

Rd,i,i=

(12)

式中：p表示所使用范數的階次；ε是一個極小的數，用以保證di-fi-1-Ji,jΔmi→0； 矩陣Rd是非奇異的；Rm與Rd具有相似的形式，其元素為

Rm,i,i=

(13)

式(11)寫成通用形式

(14)

在求解地球物理反演問題時，以Lp范數(p≥1)作為測度進行解的估計，根據觀測數據和模型參數的統計特征，使用不同的范數。如果實際地下電阻率平滑地漸變，測量誤差服從高斯分布，L2范數是最合適的方法，它對每一個元素進行平方，將向量的L2范數最小化，在數據擬合時，剔除相鄰單元的異常值，空間區域之間無明顯邊界，排除了空間結構的突變。但在實際地質結構中，如煤層、金屬礦山等的邊界會出現明顯的突變，此時，需要采用L1范數求解，它對模型電阻率值的絕對值變化最小化，對誤差大的數據不敏感，保留了相鄰元素之間的異常值。

4 約束條件

對于迭代重加權最小二乘法，每次迭代的權函數均不相同，為了在重加權函數中施加先驗約束信息，定義結構度量方程[19]

X=|mi-mt|

(15)

式中，與L1范數不同，mt在不同約束信息下代表不同的模型信息，例如相鄰單元的電導率、常數電導率等。引入關于X的重加權函數Wf，根據不同的約束要求，Wf可以用不同的函數替代。本文選擇互補誤差函數Wf1或誤差函數Wf2

(16)

(17)

式中：M是重加權的權重均值；S是權重標準偏差。重加權函數是一個基于正態的累積分布函數，其目的是在X變化時，施加不同的權重。Wf1、Wf2隨X的變化如圖1所示。

圖1 權重函數曲線

(1)L2范數平滑約束：mt代表相鄰單元的電導率，選擇Wf1，取一個較大的M值和較小的S值，當每個單元與相鄰單元的電導率差異減小且小于M+2S時，則加權函數施加權重；當電導率差異小于M-2S時，則施加滿權，使X最小化。此時的約束鼓勵兩個相鄰單元的電導率接近，等效于L2范數施加在模型空間。

(2)L1范數塊狀約束：mt代表相鄰單元的電導率，選擇Wf1，取一個較小的M值和較大的S值，此時，當每個單元與相鄰單元的電導率之差增加時，施加在相鄰單元的電導率平滑約束值減少； 當二者的差異大于M+2S時，平滑約束被移除。這樣，相鄰單元的電導率差異被保留，兩個單元的電導率值之間形成明顯邊界，等效于L1范數施加在模型空間。

(3)電導率極值約束：通過地質勘探等手段獲取地下電導率最大值和最小值的合理范圍時，反演可以施加電導率極值約束。mt代表電導率的極限值，選擇Wf1，設定M=0和一個較小的S值。此時將每個單元的電導率與mt的差值作為約束，當X小于M-2S時，權重函數提供滿權以最小化X；當X大于M+2S時，權重函數為0，移除約束，這使得反演結果中的元素電導率的最小值為mt。同樣，如果選擇Wf2，則反演結果中的元素電導率的最大值為mt。

在已知勘探區域的電阻率極大值和極小值、區域斷面位置和地下空間結構等信息的條件下，可將反演區域分為若干子區域，通過修改重加權函數的參數，增加、減少或改變每個子區域的約束信息，提高反演精度。

5 反演方法

共軛梯度法是一種常用的地球物理反演優化算法，在求解式(14)時，采用并行共軛梯度最小二乘法(Parallel Conjugate-Gradient Least-Squares,PCGLS)[12]。該方法將處理器分為主處理器和多個從處理器，反演時由主處理器分配任務給不同的從處理器，自身也同時執行計算數據量較小的任務，待從處理器完成計算任務后，主處理器再分配下一次計算任務，主、從處理器協同工作，達到提高計算效率的目的，IRLS反演流程見圖2。令式(14)的右側為b，設PCGLS的內迭代次數為j=0，1，2，…，n，λ(j)為Δm的修正因子，γ(j)是b的修正因子，數據擬合的條件是測量數據與模擬數據之間的均方誤差小于設定的容許值，算法收斂的條件是滿足數據擬合，或者內部迭代次數j達到預設值。

圖2 IRLS法迭代流程

6 合成數據反演實例

6.1 三維模型

為了評價最小二乘法三維電阻率反演效果，利用合成數據進行三維反演。應用有限元法，由TetGen軟件建立四面體非結構網格三維模型，并采用局部加密的方式，減少網格數量，提高網格精度。三維模型如圖3a所示，四面體模型由三角形網格構成，有限元空間定義為800m×800m×500m，監測區尺寸為16m×12m×10m，內部設有兩個2m×2m×2m的低阻異常體，其中低阻塊1的電導率為0.020S/m，低阻塊2的電導率為0.200S/m，圍巖電導率為0.002S/m，低阻塊埋深均為1m。測區布置3條測線，共48個電極，電極距為1.0m，測線間距為5.5m，電極分布和剖面位置見圖3b和圖3c。生成的合成數據引入5%高斯噪聲。

圖3 合成模型(a)三維模型示意圖； (b)模型水平截面圖； (c)模型y=0切面圖

為了對比混合約束條件下先驗信息的效果，分別采用平滑約束、塊狀約束和混合約束對該模型進行反演，比較不同約束信息下的反演結果。計算機CPU為 Intel i5-3470，主頻為3.2G，內存為16GB，模型網格信息見表1。

有限元正演過程采用預條件共軛梯度法計算剛度矩陣，求解空間的電場分布[7,23]。

表1 模型參數信息表

6.2 不同先驗信息約束的反演

6.2.1 平滑約束反演

對測區使用平滑約束反演。結構度量中，mt代表相鄰單元的電導率，選擇Wf1，設M=10，標準偏差S=0.1，測區及低阻塊采用同樣的約束規則。反演結果如圖4所示，可見在平滑約束下，反演結果中低阻異常體的位置基本正確，但反演精度較低，電導率從低阻的中心位置向圍巖逐漸降低，地質結構具有平滑過渡的特征；在低阻塊2中，反演低阻區域形態與模型不完全一致，低阻區向測線位置延伸。

圖4 平滑約束反演結果(a)三維切片； (b)y=0切片； (c)y=1m切片； (d)y=2m切片

6.2.2 塊狀約束反演

對測區及低阻塊施加塊狀約束反演(非連續地質結構反演)，仍采用函數Wf1，結構度量中mt代表相鄰單元的電導率，約束信息的施加通過減小權重函數的均值實現。取值M=0.1、S=0.01時的反演結果如圖5所示，可以看出，M值減小后，反演時會移除更多的相似約束，導致相鄰單元的電導率呈現較大差異，反演結果中低阻異常體位置基本準確，但精度不夠。與地質原型相比，形態出現較大偏差，低阻塊2表現得尤為明顯，同時低阻塊1和低阻塊2向測線位置偏移。

從圖4和圖5的反演結果可以看出，單純采用平滑約束(L2范數)或者塊狀約束(L1范數)的反演精度均不夠理想。在工程實踐中，由于測量誤差和人工操作的影響，反演結果的精度會更低。

6.2.3 混合約束1反演

對測區和低阻塊同時施加平滑約束、最大值約束、最小值約束進行反演(稱為混合約束1)，對測區和低阻塊分別施加約束值?？紤]實際地質結構中，同一塊地質單元電導率可能會有變化，甚至可能出現數倍的差異，因此設定某區域電導率最大值和最小值時，保留較大的余量。對每個區域同時施加三個約束條件，以測區為例，施加約束信息為：①對于平滑約束，mt代表相鄰單元的電導率，選擇Wf1，設定M=10，S=0.1； ②對于最大電導率值約束，設定mt=0.01S/m，選擇Wf1，M=0，S=0.1；③對于最小電導率值約束，設定mt=0.0004S/m，選擇Wf2，M=0，S=0.1。對低阻塊1和低阻塊2，同樣施加平滑約束和電導率極值約束信息，其中低阻塊1最大電導率值mt=0.1S/m，最小電導率值mt=0.004S/m；對于低阻塊2，設定最大電導率值mt=1.0S/m，最小電導率值mt=0.04S/m。三種約束信息同時施加的反演結果如圖6所示，可以看出，施加多個約束信息后，反演結果的精度有較大的提高，反演的測區電導率值較接近于理論模型，其中最大電導率σmax=0.03705S/m，最小電導率σmin=0.001209S/m，低阻塊的位置形態與理論模型基本一致。需要說明的是，對低阻塊2施加電導率極值約束(0.04S/m≤σ≤1.00S/m)后，實際上該區域電導率接近0.04S/m，在圖例中接近紅色區域，因此視覺上像是一個紅色方塊。

圖5 非連續地質結構約束反演結果(a)三維切片； (b)y=0切片； (c)y=1m切片； (d)y=2m切片

圖6 混合約束1反演結果(a)三維切片； (b)y=0切片； (c)y=1m切片； (d)y=2m切片

6.2.4 混合約束2反演

對低阻塊1施加平滑約束和電導率極值約束，對低阻塊2施加塊狀約束和電導率極值約束(稱為混合約束2)，方法同混合約束1，電導率極值不變。反演結果如圖7所示，可以看出，對低阻塊2的約束信息改為塊狀約束后，成像質量較好，低阻區域的位置和形態與模型高度一致。整個測區最大電導率σmax=0.04149S/m，最小電導率σmin=0.001688S/m。由于在低阻塊2中施加了塊狀約束和電導率極值約束，反演結果中低阻塊2邊界處的電導率允許出現不連續。而在低阻塊2內部，電導率極值的約束條件為0.04S/m≤σ≤1.0S/m，反演得到該區域電導率為0.04～0.04149S/m，因此低阻塊2呈現單獨的塊狀形態。與低阻塊1所施加的平滑約束相比，在不連續地質結構中，塊狀約束可以更好地反映地質體的結構變化，這在礦產資源的勘探和開采監測中尤其重要。

圖7 混合約束2反演結果(a)三維切片； (b)y=0m切片； (c)y=0.99m切片； (d)y=1.1m切片

6.3 不同約束下PCGLS法計算效率

施加不同約束信息時，PCGLS法反演性能數據統計見表2。合成算例中PCGLS算法計算時間在數分鐘至數十分鐘以內，與文獻[16](約160min)、文獻[5](約86min)的方法相比，計算速度得到較大提高。不同約束信息對迭代次數和運行時間有較大的影響，從表2可見，混合約束2反演效果優于混合約束1，且迭代次數和運行時間更少，因此在施加反演約束信息時，應根據地質條件合理施加。圖8給出了迭代過程中數據失配數量變化，迭代中第一個值為起始模型的數據失配數量。迭代過程中，最大數據失配數量為54個，不超過總數據量的2.6%。擬合均方根誤差(RMSE)為最后一次迭代后的結果，平滑約束和混合約束1由于采用了平滑約束，地質結構更加光滑，反演結果中(自然對數)電導率數據與測量數據偏差較大，均方根誤差大于1.5；而采用塊狀約束時，反演允許相鄰單元的電導率出現較大差異，保留較多的原始數據，因此擬合均方根誤差相對較小。

表2 不同約束信息的迭代次數和運行時間統計表

圖8 迭代過程中的數據失配數量

7 結論

基于迭代重加權最小二乘法實現三維電阻率反演時，將數據加權矩陣中的重加權函數修改為誤差函數和互補誤差函數的形式，通過對重加權函數參數做簡單改動，達到不同先驗信息約束的效果；反演過程中，主、從處理器多線程協同工作，共同完成反演迭代流程。合成數據的成像結果證明該方法的有效性，與單一信息約束的反演對比說明，合理施加多信息約束可以提高成像結果中異常體位置和形態的準確性；與同類型算法相比，反演耗時較短，有助于提高三維電阻率層析成像的時效。

結合其他勘探手段，獲得更多的地質結構信息，是提高直流電法三維成像效果的重要途徑，對所述方法進行工程驗證，并進一步挖掘測量數據的價值，是后續研究的重點。