于明羽,吳遵紅,譚 凱,彭代誠(chéng),熊絢兮,劉江平

(1.中國(guó)地質(zhì)大學(xué) 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430074;2.湖北特種設(shè)備檢驗(yàn)檢測(cè)研究院，湖北 武漢 430074;3.中國(guó)冶金地質(zhì)總局 中南地質(zhì)調(diào)查院,湖北 武漢 430081;4.中國(guó)地質(zhì)大學(xué) 地球物理空間與信息學(xué)院，湖北 武漢 430074)

1 引 言

隧道工程是復(fù)雜地質(zhì)地形區(qū)域交通基礎(chǔ)設(shè)施建設(shè)最優(yōu)選方案。隧道施工環(huán)境圍巖地質(zhì)情況復(fù)雜多變，經(jīng)常會(huì)因遇到巖溶、斷裂破碎帶等不良地質(zhì)體，造成塌方、突水等各種地質(zhì)災(zāi)害，因此準(zhǔn)確預(yù)報(bào)掌子面前方地質(zhì)情況，避免突發(fā)隧道地質(zhì)災(zāi)害損失，已經(jīng)成為國(guó)內(nèi)外隧道工程界關(guān)注的問題[1,2]。隧道超前探測(cè)方法，是一種全空間條件下在掌子面或迎頭后方布置觀測(cè)系統(tǒng)超前探測(cè)前方隱伏地質(zhì)災(zāi)害的地球物理技術(shù)，目前隧道施工中多采用地球物理超前預(yù)報(bào)方法進(jìn)行掌子面前方地質(zhì)預(yù)報(bào)，其中，地震波法隧道超前探測(cè)方法以不良地質(zhì)體與圍巖的波速差異為基礎(chǔ)進(jìn)行隧道地質(zhì)超前預(yù)報(bào)，具有較長(zhǎng)的探測(cè)距離及較好的地質(zhì)結(jié)構(gòu)識(shí)別效果，是隧道內(nèi)常用的超前探測(cè)地球物理方法之一[3-5]。但由于隧道地震方法超前預(yù)報(bào)觀測(cè)系統(tǒng)受隧道環(huán)境限制，采集數(shù)據(jù)量無法滿足對(duì)隧道地質(zhì)體波速計(jì)算精度的高要求，超前預(yù)報(bào)數(shù)據(jù)處理成像結(jié)果對(duì)異常體分辨能力，無法達(dá)到對(duì)不良地質(zhì)體精細(xì)探測(cè)的需求，只適用于簡(jiǎn)單地質(zhì)情況[6]，因此亟需適用于隧道超前地質(zhì)預(yù)報(bào)的精細(xì)探測(cè)手段來改變這一現(xiàn)狀。

全波形反演在地震勘探領(lǐng)域的發(fā)展，為研究者們提供了新的思路。全波形反演方法是一種基于全波場(chǎng)信息實(shí)現(xiàn)介質(zhì)參數(shù)反演的高分辨率成像方法，20世紀(jì)80年代Tarantola和Lailly將時(shí)間域全波形反演的理念引入至地震勘探領(lǐng)域[7,8]，相比常規(guī)利用單一反射波或初至波數(shù)據(jù)信息的處理方法進(jìn)行物性參數(shù)反演成像，全波形反演充分使用了全波場(chǎng)的信息，達(dá)到了高精度、高分辨率的反演成像效果，因而備受推崇。全波形反演由于非線性和周期跳躍現(xiàn)象，導(dǎo)致目標(biāo)函數(shù)存在多極值，使得反演對(duì)初始模型有嚴(yán)重依賴性，為此早期研究者們提出了對(duì)地震數(shù)據(jù)依次濾波的時(shí)間域多尺度全波形反演方法[9]來降低反演結(jié)果對(duì)初始模型的依賴。20世紀(jì)90年代，Pratt等將全波形反演的理論擴(kuò)展到頻率域，稱為頻率域全波形反演方法，頻率域全波形反演中低頻成分的反演結(jié)果可以直接作為高頻成分的初始模型直接達(dá)到多尺度反演的效果，降低對(duì)初始模型的依賴性，這一優(yōu)勢(shì)使得頻率域全波形反演在地震勘探領(lǐng)域的研究廣泛發(fā)展[10,11]。

因此近年來越來越多研究者開始嘗試將全波形反演的方法應(yīng)用于隧道空間，包括地質(zhì)雷達(dá)的隧道全波形反演、地震方法的隧道前方巖體速度反演以及隧道內(nèi)部破碎位置的全波形反演[12-14]?？紤]到近年來數(shù)學(xué)方法與地球探測(cè)理論的結(jié)合以及計(jì)算機(jī)技術(shù)的發(fā)展，更多最優(yōu)化方法被應(yīng)用于全波形反演的目標(biāo)函數(shù)極值求取中，以期得到更精確高效的反演結(jié)果，從早期Tarantola等研究者們采用的梯度下降法，以及在Mora等討論預(yù)條件處理策略的基礎(chǔ)上，Santosa采用牛頓類方法全波形反演，而近來研究者們又將擬牛頓類方法應(yīng)用到全波形反演中[15-17]。因此將全波形反演的方法應(yīng)用于隧道空間不可避免地涉及到反演優(yōu)化算法的選擇，鑒于本文將目前幾類常用方法中具有代表性的最速下降法、截?cái)嗯ｎD法以及l(fā)-BFGS方法應(yīng)用到隧道環(huán)境下的頻率域聲波全波形反演，并對(duì)比不同方法的反演迭代收斂速率、反演效果以及時(shí)間消耗，以期達(dá)到為后續(xù)研究者們提供合理建議的目的。

2 頻率域全波形反演原理及方法

2.1 頻率域全波形正演及目標(biāo)函數(shù)

二維頻率域聲波波動(dòng)方程公式如下[18]：

=-s(x,z,ω)

(1)

式(1)中：k(x,z)為體變模量，ρ(x,z)為密度，d(x,z,ω)和s(x,z,ω)分別為波場(chǎng)壓力函數(shù)與震源函數(shù)，(x,z)為空間位置向量，Δ為拉普拉斯算子。

波場(chǎng)壓力函數(shù)d(x,z,ω)與震源函數(shù)s(x,z,ω)的關(guān)系是線性的，因此離散化的二維聲波方程等式可以簡(jiǎn)化為如下的大型稀疏線性方程組：

A(ω)?d(x,z,ω)=s(x,z,ω)

(2)

式(2)中：d是波場(chǎng)矩陣，s為震源矩陣，A是關(guān)于頻率和介質(zhì)屬性的阻抗矩陣，模型空間網(wǎng)格數(shù)為(Nx,Nz)，d和s均表示(Nx,Nz)的向量，那么A(ω)是一個(gè)(Nx,Nz)×(Nx,Nz)的有限差分算子矩陣，可采用LU分解法對(duì)其進(jìn)行求解[19]。

頻率域全波形反演的目標(biāo)函數(shù)選擇為觀測(cè)波場(chǎng)與計(jì)算波場(chǎng)差值的最小二范數(shù)，具體表示為：

(3)

式(3)中：f(d)表示波場(chǎng)殘差的最小二范數(shù)，Δd表示觀測(cè)波場(chǎng)與計(jì)算波場(chǎng)之間的波場(chǎng)殘差，Nω表示反演頻率個(gè)數(shù)，Ns表示震源個(gè)數(shù)，Nr表示接收點(diǎn)個(gè)數(shù)，dobs(s,r,ω)為實(shí)際觀測(cè)值，dcal(s,r,ω)為計(jì)算波場(chǎng)。

2.2 頻率域全波形反演優(yōu)化方法

反演過程即是目標(biāo)函數(shù)梯度尋優(yōu)過程，由于全波形反演問題的高度非線性，一般采用逐步線性化的方法來減小地震記錄和模型參數(shù)之間的非線性關(guān)系[20]。在全波形反演求解極小化目標(biāo)函數(shù)的優(yōu)化迭代中，對(duì)目標(biāo)函數(shù)使用局部最優(yōu)化算法求取極值，進(jìn)行目標(biāo)函數(shù)梯度、二階導(dǎo)數(shù)等信息尋優(yōu)來解決非線性極小化問題。目標(biāo)函數(shù)梯度尋優(yōu)中多采用梯度類方法、擬牛頓類方法以及牛頓類等方法，常見的梯度類方法有最速下降法、非線性共軛梯度法；牛頓類方法主要包括牛頓法、高斯牛頓法和截?cái)嗯ｎD法；擬牛頓法則以DFP算法和l-BFGS為主要代表[21,22]。

全波形反演目標(biāo)函數(shù)最優(yōu)化的標(biāo)準(zhǔn)格式及給定初始模型x0模型情況下模型更新迭代如下：

mind∈Cf(d)dk+1=dk+αΔdk

(4)

式(4)中,C為模型的約束集合，α為迭代步長(zhǎng)，Δdk即為下降方向。

梯度下降法是目前最為常用的局部最優(yōu)化方法，梯度下降法以最速下降法(STD)為代表，采用當(dāng)前位置負(fù)梯度方向作為搜索方向，該方向?yàn)楫?dāng)前位置的最快下降方向。最速下降法計(jì)算較為簡(jiǎn)單，即：

Δdk=-f(dk)

(5)

Δf(dk)為目標(biāo)函數(shù)f(d)在dk處的梯度。

牛頓法通過構(gòu)建目標(biāo)函數(shù)對(duì)模型的二階導(dǎo)數(shù)Hessian 矩陣或者Hessian矩陣與向量的乘積來實(shí)現(xiàn)牛頓方程的線性求解[23,24]。以截?cái)嗯ｎD法(TRN)為例：

H(dk)Δdk=-f(dk)

(6)

梯度f(dk)為目標(biāo)函數(shù)f(d)在dk處的梯度，H(dk)為海森算子H(dk)=-2f(dk)。在數(shù)學(xué)算法上牛頓類方法為二階收斂，梯度下降法屬于一階收斂，牛頓法需要在每次迭代中重新求取海森算子的逆，計(jì)算步驟比梯度類方法多。

擬牛頓類方法以l-BFGS為例，相比直接的梯度計(jì)算法，l-BFGS法通過相鄰迭代次數(shù)的模型向量差和梯度向量差構(gòu)建Hessian矩陣的逆，從而實(shí)現(xiàn)對(duì)梯度的預(yù)處理作用，本質(zhì)思想是使用正定矩陣來近似Hessian矩陣的逆，簡(jiǎn)化運(yùn)算[25,26]:

Δdk=-Qκf(dk)

(7)

3 不同反演優(yōu)化方法在隧道全波形反演中的反演效果

3.1 隧道模型設(shè)計(jì)與反演參數(shù)選擇

地震隧道超前探測(cè)方法是在隧道掌子面及邊墻內(nèi)布置激發(fā)孔，在孔內(nèi)通過人工激發(fā)地震波向隧道圍巖中傳播，當(dāng)圍巖波阻抗發(fā)生變化時(shí)，一部分地震波會(huì)被反射回來，另一部分地震波將會(huì)繼續(xù)向前傳播，反射回的地震波由接收器接收并傳回電腦主機(jī)形成地震波記錄，處理后的記錄可預(yù)報(bào)隧道掌子面前方的地質(zhì)變化，為隧道施工提供可靠的地質(zhì)資料[27]。為方便討論反演結(jié)果及相關(guān)影響因素，本文設(shè)計(jì)圖1所示的隧道異常體模型和U型觀測(cè)系統(tǒng)，圖1(a)為隧道地質(zhì)模型，模型長(zhǎng)200 m，寬 30 m, 掌子面寬12 m，網(wǎng)格步長(zhǎng)1 m；根據(jù)實(shí)際條件，已開挖隧道空間填充空氣，隧道圍巖速度4 000 m/s, 在掌子面前方14 m，46 m，84 m深度分別設(shè)計(jì)三個(gè)半徑為3 m的圓形低速異常體，速度3 500 m/s，隧道地質(zhì)體內(nèi)采用均勻密度。圖1(b)為U型觀測(cè)系統(tǒng)，包含掌子面和垂直于掌子面的已開挖側(cè)邊，具體的觀測(cè)系統(tǒng)為：炮點(diǎn)從51～101 m，間距2 m，兩側(cè)邊各25炮，掌子面5炮，共55炮；檢波點(diǎn)從51～101 m，間距1 m，兩側(cè)共102道，掌子面檢波點(diǎn)13道，共115道，每次放炮時(shí)全排列接收。

圖1 隧道異常體模型及隧道觀測(cè)系統(tǒng)示意圖

全波形反演總體上是一個(gè)極小化目標(biāo)函數(shù)的優(yōu)化迭代求解過程，由于反演的非線性和周期跳躍現(xiàn)象，目標(biāo)函數(shù)存在多極值，為了提高全波形反演對(duì)初值的穩(wěn)健性，本文在隧道空間的全波形反演中采用頻率域多尺度的反演方法(圖2)。頻率域多尺度反演方法先從低頻出發(fā)開始反演，然后再逐級(jí)提高頻率進(jìn)行反演，由于對(duì)低頻數(shù)據(jù)反演陷入局部極值導(dǎo)致不收斂的概率要小，因此可以采用低頻反演結(jié)果得出相對(duì)較好的初始模型，以此提高反演過程的收斂。在進(jìn)行隧道頻率域聲波全波形反演中，以異常體模型背景速度作為反演的初始速度模型，以圖1(b)所示U型觀測(cè)系統(tǒng)作為模型波場(chǎng)計(jì)算的觀測(cè)系統(tǒng)，以200 Hz雷克子波為主頻進(jìn)行頻率域正演模擬對(duì)隧道掌子面前方低速異常體模型進(jìn)行頻率域聲波全波形反演。在多尺度全波形反演中，對(duì)垂直波數(shù)覆蓋完成度越高的頻率組，多尺度反演的效果就會(huì)越好，而垂直波數(shù)覆蓋與最大偏移距相關(guān)，在隧道空間偏移距限制的情況下，隧道頻率組無法完成垂直波數(shù)覆蓋，因此多尺度的反演效果與頻率組個(gè)數(shù)是正相關(guān)的關(guān)系[28]，頻率組內(nèi)頻率數(shù)量越多，反演效果越好。為對(duì)比不同反演優(yōu)化算法，參考地震子波頻譜范圍設(shè)計(jì)[10 Hz,500 Hz]內(nèi)自低頻到高頻頻率組個(gè)數(shù)分別為15，18和21的三組頻率(圖2b)，對(duì)每組頻率分別進(jìn)行三種不同反演優(yōu)化方法下的隧道低速異常體模型頻率域聲波全波形反演，子波頻譜和三個(gè)頻率組選擇如圖2所示。

圖2 隧道異常體正演子波頻譜及全波形反演頻率選擇示意圖

3.2 不同反演優(yōu)化方法的反演結(jié)果對(duì)比

3.2.1 不同反演優(yōu)化方法第一組頻率反演結(jié)果對(duì)比

分別采用STD、TRN以及l(fā)-BFGS方法對(duì)相同初始模型進(jìn)行圖2(b)所示第一組頻率的相同迭代次數(shù)隧道空間頻率域聲波全波形反演，經(jīng)過依次對(duì)第一組頻率的15個(gè)頻率進(jìn)行反演，最終三種反演優(yōu)化方法分別得到如圖3(c)～圖3(e)所示的速度反演結(jié)果。

圖3 隧道異常體反演模型及不同反演優(yōu)化方法計(jì)算結(jié)果

可見經(jīng)過由低至高的頻率依次反演后，三種反演算法都可以得到隧道掌子面前方低速異常體的速度成像結(jié)果，表明全波形反演方法在隧道空間是可以得到速度成像的。圖3(c)、圖3(d)以及圖3(e)反演成像結(jié)果表明,三類反演方法都能對(duì)異常地質(zhì)體成像，對(duì)淺部異常體成像效果優(yōu)于深部，隨著深度加深，異常體成像效果變差，對(duì)淺部異常體反演結(jié)果的對(duì)比顯示TRN以及l(fā)-BFGS方法的反演成像分辨率優(yōu)于STD方法。為進(jìn)一步對(duì)比三種反演優(yōu)化方法的反演效果，繪制反演過程迭代收斂曲線以及過反演結(jié)果模型中線的速度剖面曲線如圖4所示。

圖4 三種反演方法迭代誤差圖及不同反演方法結(jié)果速度剖面對(duì)比

由圖4(a)可見，不同方法反演的反演迭代收斂過程并非始終最速收斂，存在均方根誤差值的局部跳動(dòng)，參考反演過程中可能存在陷入局部極值的情況，這一現(xiàn)象反映了反演過程中搜索方向與實(shí)際誤差最小值方向會(huì)存在一定差異，同時(shí)不同反演優(yōu)化方法經(jīng)過20次迭代后的誤差值差異較小，其中l(wèi)-BFGS反演算法相對(duì)收斂速度較快，經(jīng)過20次迭代后誤差值最小。由圖4(b)可見，反演結(jié)果速度曲線中圍巖范圍速度出現(xiàn)起伏，異常體位置速度變化相對(duì)明顯，同時(shí)TRN和l-BFGS反演算法得到的異常體速度值準(zhǔn)確度高于STD反演算法結(jié)果。對(duì)反演結(jié)果剖面速度做速度平均，得到結(jié)果如圖5。

STD、TRN和l-BFGS三種反演算法反演結(jié)果速度與真實(shí)模型速度相對(duì)誤差分別為2.78 %，2.55 %，2.31 %，由圖5可見三種反演方法得到的圍巖平均速度誤差較小，異常體的平均速度誤差相對(duì)較大，其中l(wèi)-BFGS反演算法速度結(jié)果最為接近真實(shí)模型速度，綜合收斂速度分析可以得出l-BFGS方法在本模型的本次反演結(jié)果中表現(xiàn)最好。

圖5 三種不同方法反演結(jié)果平均速度剖面對(duì)比

3.2.2 不同反演優(yōu)化方法第二組頻率反演結(jié)果對(duì)比

分別采用STD、TRN以及l(fā)-BFGS方法對(duì)相同初始模型進(jìn)行圖2(b)第二組頻率的相同迭代次數(shù)隧道空間頻率域聲波全波形反演，經(jīng)過依次對(duì)第二組頻率的18個(gè)頻率進(jìn)行反演，最終三種反演優(yōu)化方法分別得到如圖6(c)～圖6(e)所示速度反演結(jié)果。

圖6 隧道異常體反演模型及不同方法反演結(jié)果

對(duì)比第一組頻率反演結(jié)果，可見第二組頻率反演結(jié)果整體上淺部異常體成像分辨率變高，深部異常體成像效果也出現(xiàn)明顯提高。由圖6(c)、圖6(d)以及圖6(e)的反演成像結(jié)果可見，三類反演方法對(duì)異常體深度的位置反演都較為準(zhǔn)確，在異常體位置與臨近位置的圍巖都有明顯速度差異，隨著深度加深，異常體成像效果變差，整體反演成像結(jié)果顯示TRN以及l(fā)-BFGS方法的反演成像分辨率優(yōu)于STD方法。為進(jìn)一步對(duì)比三種反演優(yōu)化方法的反演效果，繪制反演過程迭代收斂曲線以及過反演結(jié)果模型中線的速度剖面曲線見圖7。

圖7 三種反演方法迭代誤差圖及不同反演方法結(jié)果速度剖面對(duì)比

圖7(a)對(duì)比可見TRN以及l(fā)-BFGS的收斂曲線相對(duì)收斂速度更快，TRN以及l(fā)-BFGS方法的最終結(jié)果模型的誤差值小于STD方法的結(jié)果誤差值，l-BFGS與TRN算法反演最終誤差值近似。參考圖7(b)剖面速度曲線對(duì)比可見,背景圍巖速度存在起伏，異常體位置速度下降明顯，同時(shí)TRN和l-BFGS反演算法結(jié)果值準(zhǔn)確度高于STD反演算法結(jié)果。對(duì)反演結(jié)果剖面速度做速度平均得到結(jié)果如圖8所示。

圖8 三種不同方法反演結(jié)果平均速度剖面對(duì)比.

STD、TRN和l-BFGS三個(gè)反演算法反演結(jié)果速度與真實(shí)模型速度相對(duì)誤差分別為2.45 %，2.27 %，2.23 %，TRN和l-BFGS方法的反演結(jié)果值差異不大，l-BFGS方法在本模型的本次反演結(jié)果中表現(xiàn)相對(duì)更好。

3.2.3 不同反演優(yōu)化方法第三組頻率反演結(jié)果對(duì)比

分別采用STD、TRN以及l(fā)-BFGS方法對(duì)相同初始模型進(jìn)行圖2(b)第三組頻率的相同迭代次數(shù)隧道空間頻率域聲波全波形反演，經(jīng)過依次對(duì)第三組頻率的21個(gè)頻率進(jìn)行反演，最終三種反演優(yōu)化方法分別得到如圖9(c)～圖9(e)所示的速度反演結(jié)果。

可見隨著頻率組個(gè)數(shù)的增加，隧道環(huán)境頻率域聲波全波形反演方法得到了隧道前方地質(zhì)體較高分辨率的速度成像結(jié)果，體現(xiàn)了隧道環(huán)境下全波形反演方法高分辨率成像的優(yōu)勢(shì)。由圖9(c)、圖9(d)以及圖9(e)的反演成像結(jié)果可見，三類反演方法對(duì)異常體深度的位置反演都較為準(zhǔn)確，在異常體位置速度與圍巖速度有明顯差異，且TRN和l-BFGS優(yōu)化反演算法成像結(jié)果分辨率相比STD明顯更高。為進(jìn)一步對(duì)比三種反演優(yōu)化方法的反演效果，繪制反演過程迭代收斂曲線以及過反演結(jié)果模型中線的速度剖面曲線見圖8。

圖9 隧道異常體反演模型及不同方法反演結(jié)果

經(jīng)過對(duì)比反演迭代收斂曲線圖10(a)可見，STD方法收斂曲線跳動(dòng)較為明顯，TRN以及l(fā)-BFGS的收斂曲線相對(duì)收斂速度更快，TRN以及l(fā)-BFGS反演算法最終結(jié)果模型的誤差值小于STD方法結(jié)果誤差值，l-BFGS與TRN算法反演最終誤差值近似，但收斂速度相對(duì)較快。參考圖10(b)過掌子面中線的速度剖面對(duì)比可見結(jié)果背景圍巖速度存在波動(dòng)，反演速度在異常體位置準(zhǔn)確下降，模型淺部的反演結(jié)果精度高于深部，同時(shí)TRN和l-BFGS優(yōu)化反演算法反演得到的速度結(jié)果值準(zhǔn)確度高于STD優(yōu)化反演算法。對(duì)反演結(jié)果剖面速度做速度平均得到的結(jié)果如圖11所示。

圖10 三種反演方法迭代誤差圖及不同反演方法結(jié)果速度剖面對(duì)比

圖11 三種不同方法反演結(jié)果平均速度剖面對(duì)比

STD、TRN和l-BFGS優(yōu)化反演算法反演得到的速度反演結(jié)果相對(duì)誤差值分別為2.25 %，2.07 %，2.01 %，TRN和l-BFGS方法的反演結(jié)果值相差不大，l-BFGS方法在本模型的本次反演中相對(duì)反演精度更高。

3.2.4 不同反演優(yōu)化方法相鄰頻率組反演結(jié)果對(duì)比

對(duì)于主頻較高的頻率域多尺度全波形反演，為達(dá)到較好的反演效果，通常會(huì)選擇相對(duì)較多的頻率組來滿足波束覆蓋率，本文三種不同頻率組的反演對(duì)比也表明了同一反演算法下頻率組個(gè)數(shù)多的反演效果相對(duì)較好。為進(jìn)一步對(duì)比不同反演算法的效果，本文采用頻率個(gè)數(shù)不同的不同反演算法進(jìn)行對(duì)比，將第二組18個(gè)頻率下采用STD反演算法得到的模型速度剖面結(jié)果與第一組15個(gè)頻率下l-BFGS算法和TRN算法反演的模型速度剖面分別進(jìn)行對(duì)比，結(jié)果見圖12。

圖12 第二組頻率反演結(jié)果與第一組對(duì)比

對(duì)比圖12(a)可見，第一組15個(gè)頻率下l-BFGS算法反演得到的結(jié)果模型速度剖面對(duì)異常體的反演精度比第二組18個(gè)頻率下采用STD反演算法得到的模型速度剖面高，這表明了l-BFGS算法反演可以在頻率組頻率個(gè)數(shù)較少的情況下得到好于STD反演算法在頻率個(gè)數(shù)更多條件下的反演結(jié)果。圖12(b)顯示15個(gè)頻率下TRN算法反演得到的結(jié)果模型速度剖面精度與第二組18個(gè)頻率下采用STD反演算法得到的模型速度近似，圖12的對(duì)比結(jié)果側(cè)面表明了反演算法的選擇對(duì)反演結(jié)果精度的影響，即在其他參數(shù)相同的情況下，l-BFGS算法和TRN算法以更少的頻率個(gè)數(shù)進(jìn)行反演比STD算法以更多的頻率個(gè)數(shù)反演的反演結(jié)果表現(xiàn)更好。

將第三組21個(gè)頻率下采用STD反演算法得到的模型速度剖面結(jié)果與第二組18個(gè)頻率下l-BFGS算法和TRN算法反演的模型速度剖面分別進(jìn)行對(duì)比，結(jié)果見圖13。

圖13 第三組頻率反演結(jié)果與第二組對(duì)比

對(duì)比圖13(a)可見，第二組18個(gè)頻率下l-BFGS算法反演得到的結(jié)果模型速度剖面比第三組21個(gè)頻率下采用STD算法得到的模型速度剖面結(jié)果反演精度更高，這表明了l-BFGS算法反演可以在頻率組頻率個(gè)數(shù)較少的情況下得到明顯好于STD算法在頻率個(gè)數(shù)更多條件下的反演結(jié)果。圖13(b)顯示18個(gè)頻率下TRN算法反演得到的結(jié)果模型速度剖面與第三組21個(gè)頻率下采用STD反演算法得到的模型剖面速度誤差更小，這一對(duì)比進(jìn)一步驗(yàn)證了在其他參數(shù)相同的情況下，采用l-BFGS算法和TRN算法以更少的頻率個(gè)數(shù)進(jìn)行反演可以得到比STD算法更多頻率個(gè)數(shù)反演情況下的分辨率更高的速度反演結(jié)果。

通過對(duì)不同頻率組進(jìn)行隧道空間頻率域聲波全波形反演三種優(yōu)化算法的結(jié)果對(duì)比，綜合可得二階收斂的牛頓類方法TRN和擬牛頓類方法l-BFGS在隧道聲波頻率域全波形反演中得到的速度結(jié)果值準(zhǔn)確度明顯高于梯度類下降法的STD，同時(shí)TRN和l-BFGS優(yōu)化反演算法的隧道聲波頻率域全波形反演結(jié)果精度近似，l-BFGS優(yōu)化反演方法得到的反演結(jié)果精度相對(duì)更高；在反演迭代誤差收斂過程中TRN以及l(fā)-BFGS方法的最終結(jié)果模型的誤差值小于STD方法的結(jié)果誤差值，l-BFGS優(yōu)化反演算法的最終結(jié)果模型誤差值與TRN方法最終誤差值近似，但l-BFGS優(yōu)化反演算法的收斂速度更快。

3.3 不同反演優(yōu)化方法的并行運(yùn)算消耗對(duì)比

不同于常規(guī)地震波形反演方法只需要單一波形的數(shù)據(jù)存儲(chǔ)，全波形反演需要存儲(chǔ)數(shù)量較多的不同時(shí)刻正向傳播波場(chǎng)快照，此外二階計(jì)算的截?cái)嗯ｎD法還需要另外存儲(chǔ)反向傳播的波場(chǎng)快照，這些均涉及到大量?jī)?nèi)存開銷，因此全波形反演方法的波場(chǎng)模擬和反演對(duì)內(nèi)存的需求導(dǎo)致利用單個(gè)進(jìn)程計(jì)算已變得越來越困難[29,30]。隨著計(jì)算機(jī)技術(shù)的發(fā)展，MPI并行編程技術(shù)(分布式內(nèi)存多線程同時(shí)工作)逐漸在全波形反演中應(yīng)用開來，這一方法的使用，極大地提高了計(jì)算效率，減少了計(jì)算時(shí)間，使得全波形反演向多維及大數(shù)據(jù)量的應(yīng)用研究發(fā)展成為可能[31,32]。本文在進(jìn)行基于隧道空間頻率域聲波全波形反演的反演優(yōu)化方法對(duì)比時(shí)，將MPI并行全波形反演算法應(yīng)用于全波形反演的正演波場(chǎng)計(jì)算和反演梯度計(jì)算過程中，采用四線進(jìn)程同時(shí)工作的方式完成了三種不同反演方法的隧道空間全波形反演計(jì)算。包含有并行算法的基于隧道空間頻率域聲波全波形反演計(jì)算流程見圖14。

圖14 基于并行算法的隧道空間頻率域聲波全波形反演流程

在相同內(nèi)存條件下對(duì)相同的模型和數(shù)據(jù)采用三組頻率進(jìn)行反演，分別對(duì)三種優(yōu)化算法程序計(jì)算時(shí)長(zhǎng)對(duì)比，計(jì)算過程中STD、TRN以及l(fā)-BFGS方法反演程序計(jì)算總時(shí)長(zhǎng)具體對(duì)比曲線見圖15。

圖15 三種反演方法并行計(jì)算耗時(shí)以及計(jì)算總耗時(shí)對(duì)比

其中TRN優(yōu)化反演算法的計(jì)算時(shí)長(zhǎng)最長(zhǎng)，體現(xiàn)了二階算法在計(jì)算和存儲(chǔ)方面的更多消耗，其次是l-BFGS算法，l-BFGS算法的計(jì)算時(shí)長(zhǎng)相比TRN明顯減少，三組頻率計(jì)算時(shí)長(zhǎng)平均只相當(dāng)于TRN算法消耗時(shí)長(zhǎng)的61.2 %，而STD算法的計(jì)算時(shí)長(zhǎng)最小，平均只相當(dāng)于TRN算法消耗時(shí)長(zhǎng)的43.6 %。相同內(nèi)存下計(jì)算時(shí)長(zhǎng)的對(duì)比表明TRN優(yōu)化反演算法的計(jì)算時(shí)長(zhǎng)最長(zhǎng)，l-BFGS算法的計(jì)算時(shí)長(zhǎng)相比TRN明顯減少，而STD算法的計(jì)算時(shí)長(zhǎng)最小，節(jié)約了超過一半的時(shí)間。

4 結(jié)論與建議

本文通過將分屬于梯度類反演方法、牛頓類反演方法以及擬牛頓類反演方法的三種不同優(yōu)化方法應(yīng)用于隧道環(huán)境下的二維頻率域聲波全波形反演，對(duì)比反演效果、反演過程的迭代收斂速度以及反演過程中的時(shí)間消耗得出如下結(jié)論：

1)經(jīng)過對(duì)隧道超前探測(cè)地震觀測(cè)系統(tǒng)下的低速異常體隧道模型進(jìn)行三種不同優(yōu)化反演方法的頻率域聲波全波形反演，在頻率組選擇合適情況下三種方法均可得到隧道前方地質(zhì)體的速度成像結(jié)果，證明了隧道環(huán)境下全波形反演方法是可行的，同時(shí)三類反演方法都對(duì)異常體位置較敏感，反演結(jié)果顯示異常體位置速度與圍巖速度有明顯差異。對(duì)比可見TRN方法和l-BFGS方法反演結(jié)果成像分辨率明顯高于STD，采用TRN和l-BFGS優(yōu)化反演算法可以在頻率數(shù)量較少的情況下得到與STD反演算法較多頻率個(gè)數(shù)反演相近的結(jié)果，佐證了反演算法的選擇對(duì)反演結(jié)果的影響。

2)在隧道異常體反演迭代收斂過程中可見，TRN和l-BFGS方法比STD方法具有更快的收斂速度，同時(shí)可以得到更小的最終誤差結(jié)果，在相同的反演迭代次數(shù)中TRN和l-BFGS方法表現(xiàn)更為優(yōu)秀，l-BFGS方法的最終結(jié)果模型的誤差值與TRN方法最終誤差值近似，但l-BFGS方法收斂速度更快。

3)相同內(nèi)存下計(jì)算時(shí)長(zhǎng)的對(duì)比表明TRN優(yōu)化反演算法的計(jì)算時(shí)長(zhǎng)最長(zhǎng)，l-BFGS算法的計(jì)算時(shí)長(zhǎng)次之，STD算法的計(jì)算時(shí)長(zhǎng)最小，可見二階算法時(shí)間消耗更大。同時(shí)l-BFGS算法的計(jì)算時(shí)長(zhǎng)比TRN明顯減少，只相當(dāng)于TRN算法消耗時(shí)長(zhǎng)的61.2 %，而STD算法的計(jì)算時(shí)長(zhǎng)只相當(dāng)于TRN算法消耗時(shí)長(zhǎng)的43.6 %，節(jié)約了超過一半的時(shí)間。

本文綜合反演結(jié)果的準(zhǔn)確程度以及反演迭代過程中的收斂速度和成像精度，在相同的初始模型和計(jì)算參數(shù)下，得出l-BFGS優(yōu)化反演算法表現(xiàn)更好，TRN算法其次，但STD算法和l-BFGS算法相比TRN算法具有較為優(yōu)越的計(jì)算速度，基于全波形反演巨大的計(jì)算消耗和內(nèi)存需求，以及隧道超前探測(cè)的高分辨率成像效果需求，本文建議使用l-BFGS優(yōu)化反演算法來進(jìn)行隧道空間的頻率域全波形反演，希望這一結(jié)論可為隧道空間的全波形反演后續(xù)研究和發(fā)展提供借鑒。

致 謝

感謝對(duì)本文做出貢獻(xiàn)的Seiscope課題組開發(fā)的全波形反演算法以及對(duì)本文提出寶貴意見的課題組成員們。