一種以逆時(shí)偏移算法為引擎的波動(dòng)方程初至走時(shí)層析方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地震層析成像技術(shù)領(lǐng)域的層析速度反演技術(shù)，尤其涉及一種基于雙程 聲波方程的初至走時(shí)層析速度反演方法。
【背景技術(shù)】
[0002] 根據(jù)利用的數(shù)據(jù)信息不同，地震層析成像技術(shù)主要可分為走時(shí)層析，振幅層析及 波形層析。其中，走時(shí)層析方法最為穩(wěn)健。按正過(guò)程不同，走時(shí)層析主要可分為射線層析， 胖射線層析，菲涅爾體層析和波動(dòng)方程層析。
[0003] 射線層析（Bishop 等，1985 ;Tomographic determination of velocity and depth in laterally varying media, Geophysics, 50, 903-923)基于高頻射線理論，對(duì)初始 模型精度的要求較低，但反演精度不高。同時(shí)，由于層析矩陣十分稀疏，矩陣的零空間比較 大，因而求解時(shí)收斂較慢，且反演結(jié)果受射線照明的影響嚴(yán)重(}111和1&11^;[111?)￥1(311，2012 ;八 sensitivity controllable target-oriented tomography algorithm, 82nd SEG Annual Meeting，Expanded Abstracts，1-5)0
[0004] 考慮到地震波的傳播不僅受到射線路徑上的速度的影響，也受到射線周圍速度結(jié) 構(gòu)的影響。胖射線層析（Vasco 等，1995 ;Beyond ray tomography:Wavepaths and Fresnel volumes, Geophysics, 60 (6)，1790-1804)將射線加寬，其物理背景是波傳播是有一定寬度 的，但沒(méi)有嚴(yán)格的數(shù)學(xué)物理推導(dǎo)。此外，射線寬度的選擇也有技巧，如固定寬度（Michelena 和 Harris,1991 ;Michelena R. J. , and J.M.Harris, 1991, Tomographic traveltime inversion using natural pixels)或隨著射線傳播距離增大而逐漸增大射線寬度（Xu 等，2006 ;Enhanced tomography resolution by a fat ray technique. 76th SEG Annual Meeting, Expanded Abstracts, 3354-3358)。但這些做法不能定量描述波動(dòng)現(xiàn)象。
[0005] 相比而言，菲涅爾帶（體）層析（Yomogida, I"2 ;Fresnel zone inversion for lateral heterogeneities in the earth,Pure and Applied Geophysics, 138 (3)，391-406)在改善射線層析矩陣稀疏性的同時(shí)，還考慮了地震波在第 一菲涅爾帶（體）內(nèi)的有限頻效應(yīng)。但單頻菲涅爾帶只有在常速介質(zhì)中存在解析表達(dá)式 (Cerveny 和 Soares, 1992 ;Fresnel volume ray tracing, Geophysics, 57,902-915) 〇 在 變速介質(zhì)中，可以用常速介質(zhì)中主頻對(duì)應(yīng)的菲涅爾帶來(lái)近似有限頻帶地震波的菲涅爾帶 寬度（劉玉柱等，2009 ;初至波菲涅耳體地震層析成像，2009, 52 (9) ,2310-2320);或計(jì)算 炮點(diǎn)和檢波點(diǎn)的走時(shí)場(chǎng)，利用單頻菲涅爾帶的公式計(jì)算菲涅爾帶的范圍（Watanabe等， 1999 ；Seismic traveltime tomography using Fresnel volume approach:69thAnnual International Meeting, SEG, ExpandedAbstracts, 1402 - 1405)。但這些做法只能近似描 述變速介質(zhì)中的帶限地震波的菲涅爾帶（體）范圍。
[0006] 波動(dòng)方程走時(shí)層析沒(méi)有引入高頻近似假設(shè)，在理論上比射線類的層析具有更高 的反演分辨率。在一定條件下（波動(dòng)方程線性化近似，如一階Born/Rytov近似等）能 夠描述波場(chǎng)（或相位）擾動(dòng)與模型擾動(dòng)的定量關(guān)系。Woodward(1992 ;Wave_equation tomography:Geophysics, 57, 15-26)引入了 一 階 Born 及 Rytov 近似下"波路徑"的概 念（即波形與相位敏感度核函數(shù)）。Luo 和 Schuster (1991 ;Wave_equation traveltime inversion. Geophysics, 56, 645 - 653)利用互相關(guān)函數(shù)測(cè)量觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的時(shí)差， 并給出了互相關(guān)準(zhǔn)則下走時(shí)擾動(dòng)對(duì)模型擾動(dòng)的Fr6chet導(dǎo)數(shù)（基于一階Born近似）。 Marquering 等（1999 ;Three_dimensional sensitivity kernels for finite-frequency traveltimes: the banana-doughnut paradox. Geophys. J. Int.，I37, 8〇5_81?5)則給出了擾 動(dòng)場(chǎng)與走時(shí)擾動(dòng)的關(guān)系，進(jìn)而可以導(dǎo)出一階Born近似下走時(shí)敏感度核函數(shù)的表達(dá)。
[0007] 考慮到一階Rytov近似更適用于描述波場(chǎng)的前向散射現(xiàn)象，且適用條件較Born近 似更為廣泛。因此有必要發(fā)展一套在一階Rytov近似下以波動(dòng)方程為傳播算子的層析速度 反演技術(shù)。

【發(fā)明內(nèi)容】

[0008] 本發(fā)明的目的是針對(duì)上述現(xiàn)有技術(shù)存在的不足，提出了一種以逆時(shí)偏移算法為引 擎的波動(dòng)方程初至走時(shí)層析方法。該方法是一種高精度的背景速度反演方法，能夠反演大 尺度的宏觀背景速度場(chǎng)。此外，由于正問(wèn)題采用聲波波動(dòng)方程，沒(méi)有引入高頻近似，在理論 上比射線類的層析具有更高的反演分辨率，因而可以定量的描述有限頻帶地震波的傳播效 應(yīng)，對(duì)小尺度的速度異常有更好的反演能力。
[0009] 本發(fā)明從波動(dòng)方程線性化出發(fā)，并考慮地震波場(chǎng)的前向散射近似（一階Rytov近 似），導(dǎo)出了基于波動(dòng)方程的走時(shí)敏感度核函數(shù)在時(shí)空域的顯式計(jì)算方法，其具體算法與著 名的逆時(shí)偏移方法相似。而走時(shí)敏感度核函數(shù)就是層析線性方程組中的稀疏矩陣，通過(guò)求 解層析線性方程組，進(jìn)行速度更新，可以實(shí)現(xiàn)波動(dòng)方程初至走時(shí)層析方法。
[0010] 本發(fā)明的核心在于基于波動(dòng)方程的走時(shí)敏感度核函數(shù)的構(gòu)造，技術(shù)方案包括以下 步驟：
[0011] 步驟1 :輸入地震子波，觀測(cè)系統(tǒng)，深度域初始速度場(chǎng)，觀測(cè)地震記錄的初至走時(shí)， 以及走時(shí)殘差范數(shù)的最小值；
[0012] 步驟2 :利用聲波波動(dòng)方程正演方法，計(jì)算初始速度場(chǎng)中正演的模擬地震記錄，求 解如下波動(dòng)方程；
[0013] 步驟3 :計(jì)算模擬地震記錄的初至波到達(dá)時(shí)；
[0014] 步驟4:計(jì)算初至波走時(shí)殘差，判斷走時(shí)殘差的范數(shù)是否小于預(yù)先設(shè)定的精度要 求，若滿足要求，跳轉(zhuǎn)至步驟8 ;否則，執(zhí)行下一步；
[0015] 步驟5 :計(jì)算波動(dòng)方程走時(shí)敏感度矩陣；
[0016] 步驟6 :求解層析線性方程組；
[0017] 步驟7 :更新速度模型，進(jìn)行下一次迭代反演，跳轉(zhuǎn)回步驟2 ;
[0018] 步驟8 :輸出反演得到的速度模型。
[0019] 上述技術(shù)方案進(jìn)一步優(yōu)化為：
[0020] 步驟1 :輸入地震子波f(t)，觀測(cè)系統(tǒng)，深度域初始速度場(chǎng)v(x)，觀測(cè)地震記錄的 初至走時(shí)T°bs (\，xs)，走時(shí)殘差L2范數(shù)的最小值ε。
[0021] 步驟2 :利用聲波波動(dòng)方程正演方法，計(jì)算初始速度場(chǎng)中正演的模擬地震記錄 ural (xr, 11 xs)，求解如下波動(dòng)方程：
[0023] 其中，&和\分別代表觀測(cè)系統(tǒng)中檢波器和震源的坐標(biāo)，t代表時(shí)間，Λ代表拉普 拉斯算子。
[0024] 步驟3 :計(jì)算模擬地震記錄的初至波到達(dá)時(shí)ral (x" xs)。
[0025] 步驟4 :計(jì)算初至波走時(shí)殘差Λ T (\，xs)，按照如下公式：
[0026] Δ T (xr，xs) = Tobs (xr，xs) -Tcal (xr，xs)。 （8)
[0027] 其中，T°bs(\，xs)代表輸入的觀測(cè)數(shù)據(jù)的初至波到達(dá)時(shí)。
[0028] 同時(shí)，判斷走時(shí)殘差的L2范數(shù)是否小于預(yù)先設(shè)定的精度要求（ )〇
[0029] 若滿足要求，跳轉(zhuǎn)至步驟8。否則，執(zhí)行下一步。
[0030] 步驟5 :計(jì)算波動(dòng)方程走時(shí)敏感度矩陣K(x | xs)，按照如下公式：
[0032] 其中，uQ(x|xs ;t)為震源在\處的地震子波產(chǎn)生的地震波場(chǎng)，X為地下空間任意一 點(diǎn)的坐標(biāo)。上標(biāo)?代表時(shí)間導(dǎo)數(shù)。P(x|\;t)為檢波器處逆時(shí)傳播的地震波場(chǎng)，表示為如下 公式
[0033] p (x | xr ；t) = g〇 (x, -t I xr, 0) *u〇 (xr | xs ；t) 〇 (10)
[0034] 其中，g。(x, _t I 0)為聲波波動(dòng)方程的格林函數(shù)，代表在震源Xl^處，從時(shí)間t = 0 開始的逆時(shí)傳播過(guò)程，UQ(X」XS ;t)為正演波場(chǎng)uQ(x|xs ;t)在檢波器位置Χρ處的地震記錄， 即
[0036] 顯然，波動(dòng)方程走時(shí)敏感度矩陣K(x|xpXs)的計(jì)算與業(yè)內(nèi)著名的波動(dòng)方程逆時(shí)偏 移有相類似的計(jì)算流程，即：震源處正演的模擬波場(chǎng)與檢波器處逆時(shí)傳播的地震記錄對(duì)時(shí) 間的積分。不同之處在于，公式（9)中需要的是正演的模擬波場(chǎng)的時(shí)間導(dǎo)數(shù)，并且需要能量