基于光滑約束正則化的電磁波層析反演

林芳鵬，戴前偉，雷 軼

(1.中南大學 有色金屬成礦預測教育部重點實驗室，長沙 410083;2.廣西交通工程檢測有限公司，南寧 530011；3.中南大學 地球科學與信息物理學院，長沙 410083)

0 引言

層析成像技術(Computed Tomography，CT)是反演孔間電磁波數據的主要方法之一，其分為走時層析成像和衰減層析成像[1-3]。電磁波層析成像原理基于Radon變換，即電磁波初至時或振幅可用射線經過的各單元的慢度或衰減系數的線性積分表示。射線路徑與慢度或衰減系數有關，且射線覆蓋角度與網格剖分大小直接影響射線長度，因此層析反演的實質是解病態方程組[4]。

正則化可以有效地解決地球物理中的不適定性問題，通過將反問題的先驗信息和對解估計的預期，加入到反演計算中以消除方程病態性，從而保證計算的穩定性[4-6]。常用的正則化約束條件有3種：①最小模型約束；②最平緩模型約束；③最光滑模型約束[7-8]，筆者采用最光滑模型約束，即模型參數二階導數最小。正則化參數的選擇直接決定了正則化的有效性，常用的選擇方法分為2類：①已知原始數據的誤差水平時可以采用Morozov偏差原理、廣義偏差原理、Arcangeli準則等方法；②不能預先估計原始數據誤差水平時，擬最優準則、廣義交叉檢驗準則、L曲線法等都是有效的方法[8-11]。由于層析反演中初至時或振幅值的誤差水平都是未知的，故筆者基于L曲線法選取正則化參數。層析反演的常用算法有(反投影技術)BPT、(代數重建技術)ART、(聯合迭代重建技術)SIRT、共軛梯度法(CG)、最小二乘QR分解法(LSQR)等[6,12]。因為LSQR在解大型稀疏矩陣方程組時具有節約存儲空間、運算速度快的特點，所以與其他算法相比LSQR更適合解決病態問題[13]。

1 光滑約束正則化LSQR算法

層析成像求解的基本方程組為式(1)。

Gm=d

(1)

假設有n條射線，p個網格單元。G為n×p的系數矩陣，包含第ni條射線在第pj個網格單元中的射線長度；m為p×1的慢度矩陣或衰減系數矩陣；d為n×1的走時矩陣或振幅矩陣。

由于式(1)方程組具有病態性，為保證解的穩定性需要加入正則化：

‖Gm-d‖2+λ‖Dm‖2

(2)

式中：λ為正則化參數；D是光滑約束矩陣，即二階導數型正則化。式(2)等價于解最小二乘問題：

(3)

Ax=b

(4)

因為LSQR法具有節約存儲空間、計算速度快、數值穩定性好等特點，所以應用LSQR法解方程組式(4)。LSQR法遵循最小二乘的原則，在限定的解空間內尋求最佳解，它包含了Lanczos迭代和Givens正交變換兩個過程，其詳細的推導可參考文獻[14-15]。LSQR算法過程如下：

1)初始化。

2)循環計算k=1、2、3、…重復以下步驟。

3)Lanczos雙對角化。

βk+1uk+1=Avk-αkuk，

αk+1vk+1=ATuk+1-βk+1vk。

4)正交變換。

5)迭代計算。

6)收斂判別。若迭代次數達到最大值或方程組的解隨迭代次數增加無明顯變化時，停止迭代。

2 L曲線法

正則化參數λ的選擇直接影響層析成像的精度[16]。λ過小時強調數據擬合，數據噪音可能會被放大而出現過擬合的現象；λ過大時強調光滑度而忽視了數據擬合，數據噪音被忽略而出現過抑制的現象。

L曲線法基于偏差函數‖Gm-d‖和約束項‖Dm‖在對數坐標上的平面曲線，該方法通過求解(lg‖Gm-d‖，lg‖Dm‖)曲線的最大曲率點(拐角點)來確定最佳正則化參數[17-18]。求取最大曲率點的推導公式如下，令：ρ=lg‖Gm-d‖，η=lg‖Dm‖，那么正則化參數λ對應的點的曲率為：

(5)

其中：ρ'、ρ''、η'、η''分別表示ρ、η的一階和二階導數，且均是正則化參數λ的函數。式(5)最大值對應的正則化參數λ即為最佳正則化參數，此時的L曲線拐點是偏差函數和約束項的最佳平衡點。

3 計算實例

3.1 簡單模型

圖1為簡單模型示意圖，模型為5 m×10 m，背景的相對介電常數為20；異常體邊長為1 m，相對介電常數為6。數據采集為單發多收模式，發射天線置于X=0 m，垂直移動間隔為0.5 m；接收天線置于X=5 m，垂直移動間隔為0.25 m。計算實例均基于走時層析成像，采用直線射線路徑計算系數矩陣G，利用能量比法提取電磁波初至時。

圖1 簡單模型圖Fig.1 Simple model

圖2 簡單模型正則化參數的L曲線Fig.2 L-curve for computing the optimal regularization parameter in simple model

正則化參數λ選用0.01、0.05、0.5、1、2、5、10，代入正則化目標函數中得到7組不同的層析反演結果。求取每個正則化參數對應的偏差函數值lg‖Gm-d‖和約束項lg‖Dm‖，利用三次樣條擬合得到L曲線(如圖2)，拐點處正則化參數λ=0.5。

如圖3列舉了正則化參數λ為0.05、0.5、5時的三組層析反演結果，由圖3可知，不同的正則化參數得到的反演結果有明顯差別。圖3(a)為正則化參數λ=0.05時的電磁波速度剖面圖，由于正則化參數選取過小，背景電磁波速度出現波動，且異常體的特征反映較為模糊；圖3(b)為最佳正則化參數λ=0.5時的電磁波速度剖面，從圖中可清晰識別異常體的位置和大小，異常體的高速特性也得以準確體現；圖3(c)為正則化參數λ=5時的電磁波速度剖面，當正則化參數選取過大時，高速異常區過于光滑，導致層析成像失真。

3.2 斷層模型

圖4為斷層模型示意圖，模型為5 m×10 m，背景的相對介電常數為15；異常體的相對介電常數為25。數據采集以及反演參數設置同上。

正則化參數λ選用0.25、0.5、1、2、4、8、16，代入正則化目標函數中得到7組不同的層析反演結果，其對應的L曲線如圖5所示，拐點處正則化參數λ=8。

圖6為正則化參數λ為2、8、16時的三組層析反演結果。圖6(a)為正則化參數λ=2時的電磁波速度剖面圖，由于正則化參數選取過小，方程組的解穩定性較差，導致反演結果與實際模型誤差較大；圖6(b)為最佳正則化參數λ=8時的電磁波速度剖面圖，從圖6(b)中可清晰識別低速斷層的傾向、大小和位置；圖6(c)為正則化參數λ=16時的電磁波剖面圖，當正則化參數選取過大時，斷層異常區過于圓滑、數據擬合較差，導致難以從圖像中識別異常體特征。

圖4 斷層模型Fig.4 Fault model

圖5 斷層模型正則化參數L曲線Fig.5 L-curve for computing the optimal regularization parameter in fault model

圖6 斷層模型的層析反演剖面圖Fig.6 The results of tomography inversion in fault model(a) λ=2；(b) λ=8；(c) λ=16

4 總結

筆者提出了基于光滑約束正則化的LSQR層析反演算法，介紹了L曲線法選取最佳正則化參數。由簡單模型和斷層模型的走時層析成像結果可知：①基于光滑約束正則化的LSQR算法適用于層析成像計算，它具有穩定性和高效性；②利用L曲線的拐點可準確選擇最佳正則化參數，其對應的速度剖面圖能清晰反映異常體的特征。

層析成像的影響因素有多種，需要指出，本文中直線射線路徑以及射線覆蓋程度都會造成一定程度的誤差，導致電磁波速度剖面中異常體的位置、大小與實際模型有所不同。

參考文獻：

[1] 趙連鋒，朱介壽，嚴忠瓊，等．巖溶勘察層析成像試驗：成像效果分析[J]．物探化探計算技術，2002，24(4)：332-336．

ZHAO L F, ZHU J S, YAN Z Q, et al. A test of tomography in the cavern survey: imaging effect analysis [J]. Computing Techniques for Geophysical and Geochemical Exploration, 2002, 24(4):332-336. (In Chinese)

[2] 馮銳．電磁波層析成象[J]．地震學報，1997，19(5)：524 - 534．

FENG R.Electromagnetic wave tomography[J]. Journal of seismic，1997, 19(5):524-534. (In Chinese)

[3] 陳昌彥，沈小克，蘇兆鋒，等．電磁波層析成像技術在復雜地質邊坡工程勘察中的應用研究[J]．地球物理學進展，2012，27(2)：796-803．

CHEN C Y, SHEN X K, SHU Z F, et al. Application of electromagnetic wave tomography on geotechnical investigation of the complex geological slope engineering [J]. Progress in Geophysics, 2012, 27(2): 796-803. (In Chinese)

[4] 成谷, 張寶金, 馬在田，等. 淺談反射地震走時層析中的正則化[J]. 地球物理學進展, 2007, 22(6):1715-1721.

CHENG G, ZHANG J B, MA Z T, et al. Regularization in traveltime tomography for reflection seismic data [J]. Progress in Geophysics, 2007, 22(6): 1715-1721. (In Chinese)

[5] 董一飛, 張致付. 正則化方法在無線電波層析中的應用[J]. CT理論與應用研究, 2016, 25(3):287-298.

DONG Y F, ZHANG Z F. Application of regularization method in radio wave tomography [J]. CT Theory and Applications, 2016, 25(3): 287-298. (In Chinese)

[6] 楊文采. 地球物理反演和反射地震走時層析成像[M].北京:地質出版社, 1989.

YANG W C. Geophysical inversion and seismic tomography [M] Beijing: Geological Publishing House , 1989. (In Chinese)

[7] 柳建新，童孝忠，郭榮文，等. 大地電磁測深勘探：資料處理反演與解釋[M]. 北京：科學出版社，2012.

LIU J X, TONG X Z, GUO R W, et al. Magnetotelluric sounding: data inversion and interpretation [M]. Beijing: Science Press, 2012. (In Chinese)

[8] 張文權, 翁愛華. 地面核磁共振正則化反演方法研究[J]. 吉林大學學報(地球科學版), 2007, 37(4):809-813.

ZHANG W Q, WENG A H. On regularization inversion method of surface nuclear magnetic resonance data [J]. Journal of Jilin University (Earth Science Edition), 2007, 37(4): 809-813. (In Chinese)

[9] 王化祥, 何永勃, 朱學明，等. 基于L曲線法的電容層析成像正則化參數優化[J]. 天津大學學報(自然科學與工程技術版), 2006, 39(3):306-309.

WANG H X, HE Y B, ZHU X M, et al. Regularization Parameter optimum of electrical capacitance tomography based on L-curve method [J]. Journal of Tianjin University, 2006, 39(3):306-309. (In Chinese)

[10] 馮德山, 王珣. 大地電磁雙二次插值FEM正演及最小二乘正則化聯合反演[J]. 中國有色金屬學報, 2013(9):2524-2531.

FENG D S, WANG X. Magnetotelluric finite element method forward based on biquadratic interpolation and least squares regularization joint inversion [J]. The Chinese Journey of Nonferrous Metals, 2013, 23(9): 2524-2531. (In Chinese)

[11] 馬濤, 陳龍偉, 吳美平,等. 基于L曲線法的位場向下延拓正則化參數選擇[J]. 地球物理學進展, 2013, 28(5):2485-2494.

MA T, CHEN L W, WU M P, et al. The selection of regularization parameter in downward continuation of potential field based on L-curve method[J]. Progress in Geophysics, 2013, 28(5): 2485-2494. (In Chinese)

[12] 楊薇, 劉四新, 馮彥謙. 跨孔層析成像LSQR算法研究[J]. 物探與化探, 2008, 32(2):199-202.

YANG W, LIU S X, FENG Y Q, et al. A study of the LSQR algorithm for cross-hole tomography [J]. Geophysical and Geochemical Exploration, 2008, 32(2):199-202. (In Chinese)

[13] 劉勁松, 劉福田, 劉俊,等. 地震層析成像LSQR算法的并行化[J]. 地球物理學報, 2006, 49(2):540-545.

LIU J S, LIU F T, LIU J, et al. Parallel LSQR algorithms used in seismic tomography [J]. Chinese Journal of Geophysics, 2006 , 49(2): 540-545.(In Chinese)

[14] PAIGE C C, SAUNDERS M A. LSQR: An algorithm for sparse linear equations and sparse least squares[J]. Acm Transactions on Mathematical Software, 1982, 8(1):43-71.

[15] 張東, 喬友峰, 姜麟舜,等. 地震層析成像中LSQR算法的快速求解[J]. 物探化探計算技術, 2011, 33(6):632-635.

ZHANG D, QIAO Y F, JIANG L Y, et al. The fast algorithm of LSQR in seismic tomography [J]. Computing Techniques for Geophysical and Geochemical Exploration, 2011, 33(6):632-635. (In Chinese)

[16] 冉利民, 劉四新, 李玉喜,等. 影響跨孔雷達層析成像效果的幾個因素[J]. 吉林大學學報(地球科學版), 2013, 43(5):1672-1680.

RAN L M, LIU S X, LI Y X, et al. Several factor affecting cross-hole tomography [J]. Journal of Jilin University (Earth Science Edition), 2013, 43(5):1672-1680. (In Chinese)

[17] HANSEN P C. Analysis of discrete Ill-posed problems by means of the L-curve [J]. Siam Review, 1992, 34(4):561-580.

[18] HANSEN P C, O'LEARY D P. The use of the L-curve in the regularization of discrete ill-posed problems[J]. Siam Journal on Scientific Computing, 1993, 14(6):1487-1503.