林芳鵬，戴前偉，雷 軼

(1.中南大學(xué) 有色金屬成礦預(yù)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083;2.廣西交通工程檢測(cè)有限公司，南寧 530011；3.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

0 引言

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

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

1 光滑約束正則化LSQR算法

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

Gm=d

(1)

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

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

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

(2)

式中：λ為正則化參數(shù)；D是光滑約束矩陣，即二階導(dǎo)數(shù)型正則化。式(2)等價(jià)于解最小二乘問(wèn)題：

(3)

Ax=b

(4)

因?yàn)長(zhǎng)SQR法具有節(jié)約存儲(chǔ)空間、計(jì)算速度快、數(shù)值穩(wěn)定性好等特點(diǎn)，所以應(yīng)用LSQR法解方程組式(4)。LSQR法遵循最小二乘的原則，在限定的解空間內(nèi)尋求最佳解，它包含了Lanczos迭代和Givens正交變換兩個(gè)過(guò)程，其詳細(xì)的推導(dǎo)可參考文獻(xiàn)[14-15]。LSQR算法過(guò)程如下：

1)初始化。

2)循環(huán)計(jì)算k=1、2、3、…重復(fù)以下步驟。

3)Lanczos雙對(duì)角化。

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

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

4)正交變換。

5)迭代計(jì)算。

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

2 L曲線法

正則化參數(shù)λ的選擇直接影響層析成像的精度[16]。λ過(guò)小時(shí)強(qiáng)調(diào)數(shù)據(jù)擬合，數(shù)據(jù)噪音可能會(huì)被放大而出現(xiàn)過(guò)擬合的現(xiàn)象；λ過(guò)大時(shí)強(qiáng)調(diào)光滑度而忽視了數(shù)據(jù)擬合，數(shù)據(jù)噪音被忽略而出現(xiàn)過(guò)抑制的現(xiàn)象。

L曲線法基于偏差函數(shù)‖Gm-d‖和約束項(xiàng)‖Dm‖在對(duì)數(shù)坐標(biāo)上的平面曲線，該方法通過(guò)求解(lg‖Gm-d‖，lg‖Dm‖)曲線的最大曲率點(diǎn)(拐角點(diǎn))來(lái)確定最佳正則化參數(shù)[17-18]。求取最大曲率點(diǎn)的推導(dǎo)公式如下，令：ρ=lg‖Gm-d‖，η=lg‖Dm‖，那么正則化參數(shù)λ對(duì)應(yīng)的點(diǎn)的曲率為：

(5)

其中：ρ'、ρ''、η'、η''分別表示ρ、η的一階和二階導(dǎo)數(shù)，且均是正則化參數(shù)λ的函數(shù)。式(5)最大值對(duì)應(yīng)的正則化參數(shù)λ即為最佳正則化參數(shù)，此時(shí)的L曲線拐點(diǎn)是偏差函數(shù)和約束項(xiàng)的最佳平衡點(diǎn)。

3 計(jì)算實(shí)例

3.1 簡(jiǎn)單模型

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

圖1 簡(jiǎn)單模型圖Fig.1 Simple model

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

正則化參數(shù)λ選用0.01、0.05、0.5、1、2、5、10，代入正則化目標(biāo)函數(shù)中得到7組不同的層析反演結(jié)果。求取每個(gè)正則化參數(shù)對(duì)應(yīng)的偏差函數(shù)值lg‖Gm-d‖和約束項(xiàng)lg‖Dm‖，利用三次樣條擬合得到L曲線(如圖2)，拐點(diǎn)處正則化參數(shù)λ=0.5。

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

3.2 斷層模型

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

正則化參數(shù)λ選用0.25、0.5、1、2、4、8、16，代入正則化目標(biāo)函數(shù)中得到7組不同的層析反演結(jié)果，其對(duì)應(yīng)的L曲線如圖5所示，拐點(diǎn)處正則化參數(shù)λ=8。

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

圖4 斷層模型Fig.4 Fault model

圖5 斷層模型正則化參數(shù)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 總結(jié)

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

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

參考文獻(xiàn)：

[1] 趙連鋒，朱介壽，嚴(yán)忠瓊，等．巖溶勘察層析成像試驗(yàn)：成像效果分析[J]．物探化探計(jì)算技術(shù)，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]．地震學(xué)報(bào)，1997，19(5)：524 - 534．

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

[3] 陳昌彥，沈小克，蘇兆鋒，等．電磁波層析成像技術(shù)在復(fù)雜地質(zhì)邊坡工程勘察中的應(yīng)用研究[J]．地球物理學(xué)進(jìn)展，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] 成谷, 張寶金, 馬在田，等. 淺談反射地震走時(shí)層析中的正則化[J]. 地球物理學(xué)進(jìn)展, 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] 董一飛, 張致付. 正則化方法在無(wú)線電波層析中的應(yīng)用[J]. CT理論與應(yīng)用研究, 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] 楊文采. 地球物理反演和反射地震走時(shí)層析成像[M].北京:地質(zhì)出版社, 1989.

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

[7] 柳建新，童孝忠，郭榮文，等. 大地電磁測(cè)深勘探：資料處理反演與解釋[M]. 北京：科學(xué)出版社，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] 張文權(quán), 翁愛(ài)華. 地面核磁共振正則化反演方法研究[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 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] 王化祥, 何永勃, 朱學(xué)明，等. 基于L曲線法的電容層析成像正則化參數(shù)優(yōu)化[J]. 天津大學(xué)學(xué)報(bào)(自然科學(xué)與工程技術(shù)版), 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正演及最小二乘正則化聯(lián)合反演[J]. 中國(guó)有色金屬學(xué)報(bào), 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曲線法的位場(chǎng)向下延拓正則化參數(shù)選擇[J]. 地球物理學(xué)進(jìn)展, 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]. 地球物理學(xué)報(bào), 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]. 物探化探計(jì)算技術(shù), 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] 冉利民, 劉四新, 李玉喜,等. 影響跨孔雷達(dá)層析成像效果的幾個(gè)因素[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 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.