基于偽解析法的TTI介質純qP波地震波場正演模擬

張奎濤,廖家榮,顧漢明,孫瑛瑩,陳懌旸,王凱

(1.江西省交通設計研究院有限責任公司,江西 南昌 330052; 2.中國地質大學(武漢) 地球物理與空間信息學院,湖北 武漢 430074)

0 引言

實際地下介質充滿流體、裂縫及裂隙,因而表現出各向異性特征,這種現象也在實驗室和波場測量中被觀測到[1],使得各向異性成為地震波場數值模擬中越來越不容忽略的因素。

進行彈性波數值模擬時,得到的混合波場通常同時含有縱波和橫波,而在進行地震偏移成像,尤其是彈性波逆時偏移成像時,須將縱、橫波解耦,以得到物理意義明確的成像剖面[2]。為此,Alkhalifah[3-4]提出了聲學近似假設下的擬聲波方程,人為設定沿對稱軸方向的橫波速度為零,推導出標量形式的四階微分方程。為了簡化方程,提高計算效率,Du等[5]和Zhou等[6]通過不同的降階方法,推導出不同形式的二階耦合qP波微分方程。Duveneck等[7]從胡克定律和運動方程出發,推導得到另一種形式的擬聲波方程,且其波場具有明確的物理意義。程玖兵等[8-9]推導得到各向異性介質的偽純模式波動方程,使其在運動學上同彈性波動方程等價。

上述各方法均存在擬聲波方程的突出缺點,即當模型參數中有變化較大的傾角和方位角時,會出現不穩定的現象[10]、殘留著橫波假象[11],或只有在Thomsen參數滿足ε≥δ時正演模擬才是穩定的[12]。為了解決此類問題,有學者另辟蹊徑,從qP波與qSV波的頻散關系出發,直接解耦以構建純qP波控制方程。要實現qP與qSV波的完全分離,關鍵在于導出簡單、靈活的純qP波控制方程和設計有效的擬微分求解算法。例如,Du等[5]和Zhang等[13]基于弱各向異性近似和平方根近似,推導出一種純qP波解耦方程,其表現為時間—波數域形式,能較準確地描述TTI介質中的運動學特征;Zhan等[14]推導了一種新的解耦方程,將P波與SV波完全分離開來,并成功將其應用于逆時偏移中;Xu等[15]提出了一種純qP波橢圓微分方程,通過將微分方程分解成一個標量算子和一個微分算子,并用橢圓分解方法修正標量算子,達到校正振幅誤差的目的;胡書華等[16]則將Xu等[15]二階微分方程做降階處理,通過Lebedev交錯網格,得到TTI介質下純qP波的穩定傳播;張慶朝等[17]利用近似相速度公式,導出TTI介質三維qP波頻散方程及波動方程,并使用偽譜法做正演模擬;顧漢明等[18]將各向異性介質推廣到黏聲各向異性介質,通過Low-rank方法,實現了黏聲各向異性介質純qP波正演模擬。

純qP波波動方程中往往存在著擬微分算子,譜方法能夠較好地進行求解。偽譜法(pseudo spectral method,PSM)是一種早期發展出來的譜方法[19],其基本思想是通過傅里葉變換方法計算波場對空間的導數,在時間域通過有限差分法計算波場對時間的導數。正因如此,使得偽譜法具有非常高的空間精度,沒有誤差累積效應和空間頻散,是一種較為理想的數值模擬方法。由于偽譜法計算時間導數時采用的是有限差分,使得其具有時間頻散,為了解決這一問題,Etgen等[20]提出了偽解析算法,其基本思想是波數域推導時間誤差補償因子并作用于拉普拉斯算子,從而補償在偽譜法中因時間導數的二階差分造成的誤差,達到進一步提高數值模擬精度的目的;Chu等[21]在Etgen提出的偽解析法的基礎上,提出了歸一化偽拉普拉斯算子的偽解析法,隨后,張衡等[22]將這一方法運用于解耦的TTI介質波動方程中,得到了高精度的地震波場,在時間和空間上能夠同時達到Nyquist頻率。

本文在前人的研究基礎上,將qP波擬微分方程變換到空間—波數域,并通過坐標變換,推導了時間域TTI介質二階純qP波波動方程,引入偽解析算法,實現了基于偽解析法的TTI介質純qP波地震波場正演模擬。通過簡單和復雜模型測試,驗證了本文方法的正確性及對復雜介質的適用性。

1 時間域TTI介質二階純qP波波動方程

為了推導本文的TTI介質純qP波方程,首先從如下VTI介質qP波標量方程

(1)

出發[23]。式中,Q為擬微分算子,其具體表達式為

(2)

式中:t為時間坐標;x、y、z為空間坐標;p為地震波場;v0為qP波垂直速度;ε、δ為Thomsen參數。

式(1)控制著VTI介質中qP波的傳播,其中所有參數都是隨空間變化的,且與混合波場中縱波有著一樣的頻散關系,即相位相同,能夠準確地表征其運動學特征。

將式(1)變換到頻率—波數域中有

(3)

將式(3)改寫為

(4)

(5)

將式(4)重新改寫為

(6)

記橢圓微分算子Se為

(7)

則將式(6)變換到時間域有

(8)

二維情況下有

(9)

(10)

式中:θ為TTI介質中的極化角;Gx、Gz分別為x、z方向經過坐標旋轉后新的一階微分算子;Hx、Hz分別為x、z方向經過坐標旋轉后新的二階微分算子。

將式(10)代入式(9)中即可得到二維TTI介質純qP波方程:

(11)

式(11)克服了擬聲波方程的局限性,消除了偽橫波干擾,不受模型參數限制且地震波場能穩定傳播,能夠準確地表征地震波場的運動學特征。

偽譜法可以很好求解式(11),即有

(12)

(13)

(14)

(15)

2 偽解析法

在用偽譜法進行數值模擬時,雖然在空間方向能達到一個波長兩個采樣點的Nyquist譜精度,但波場對時間的導數仍用有限差分方法計算。因此,在時間方向上仍存在時間頻散。而另一種譜方法,即偽解析方法(pseudo analytical method,PAM)能夠較好地解決上述問題,并且相比于偽譜法,精度能夠進一步提高,使得空間和時間方向上均無頻散。其主要思想是在波數域通過修正拉普拉斯算子來達到對時間方向上的有限差分離散誤差進行補償的目的。

對于二階聲波方程,有

(16)

其二階時間離散形式有

本文提出一種抗盲檢測直擴隱蔽信號設計方法，提出基于數據分級的大信號掩蓋技術，給出了波形參數設計方案，并對大信號掩蓋下的機密信號解調BER的理論值進行推導.最后仿真驗證了所提方案可實現機密信號的抗盲檢測.同時在保證機密信號抗盲檢測能力的情況下，解調損失可控制在接受范圍內.未來工作將分析初始相位、擴頻碼等參數對信號抗截獲的影響.

(17)

對上式進行傅里葉變換并利用時移定理有

(18)

化簡上式并代入ω=v0|k|有

(19)

因此,式(17)的偽解析法為

(20)

定義歸一化偽拉普拉斯(NPL)算子[21]有

(21)

則式(17)的NPL偽解析法為

(22)

故,式(11)的偽解析法格式為

(23)

偽解析法能夠準確地補償時間方向二階差分所造成的誤差,得益于其一個重要特性,即歸一化偽拉普拉斯算子的空間導數項隨速度的變化緩慢,使得其能夠較為準確地補償時間方向上的誤差。為了說明其空間導數項隨速度變化緩慢的特點,采用速度1 500 m/s和3 000 m/s作為對比,時間步長為0.001 s,選擇波數范圍為[0,0.21]。圖1為不同速度時的歸一化偽拉普拉斯(NPL)算子;圖2為基于歸一化偽拉普拉斯算子的各類空間二階導數,包括x方向二階導數、z方向二階導數、xz方向混合二階導數及擬微分算子;圖3為基于Hess復雜介質模型的歸一化偽拉普拉斯算子。從圖1～3中可以發現,歸一化偽拉普拉斯算子及各類空間二階導數均隨速度變化緩慢,能夠適應較復雜介質。

圖1 歸一化偽拉普拉斯(NPL)算子Fig.1 Normalized pseudo-Laplace(NPL) operator

圖2 基于歸一化偽拉普拉斯(NPL)算子的空間二階導數Fig.2 A second derivative of space based on normalized pseudo-Laplace(NPL) operators

3 模型試算

3.1 均勻模型

為了驗證本文構建的二維TTI介質純qP波方程不受模型參數的限制,且無偽橫波干擾,設計了尺寸為2 000 m×2 000 m的均勻介質模型。其空間網格大小為5 m×5 m,時間步長為1 ms,總時間采樣點數為1 000,記錄長度為1 s;縱波震源位于模型中心,采用主頻為25 Hz的Ricker子波震源;邊界條件為海綿邊界條件;接收排列范圍是0～2 000 m,排列深度為750 m,道間距為10 m,共201道接收。其均勻介質彈性參數見表1。

表1 均勻介質模型參數Table 1 Homogeneous media model parameters

圖4上部為聲學TTI介質常規擬聲波方程分別在介質類型Ⅲ(圖4a)、Ⅳ(圖4b)、Ⅴ(圖4c)時400 ms時刻波場快照,圖4下部為聲學TTI介質純qP波方程分別在介質類型Ⅲ(圖4d)、Ⅳ(圖4e)、Ⅴ(圖4f)時400 ms時刻波場快照。從圖中可見常規擬聲波方程(即直接將橫波速度設為0)在ε>δ時雖能穩定傳播,但會出現菱形假象,嚴重干擾了有效波;在ε<δ時會出現不穩定現象;只有在ε=δ時才能穩定傳播且無假象,這樣就極大地限制了其應用。而本文構建純qP波方程(圖4下部),無論是ε<δ、ε>δ還是ε=δ時,都能穩定地傳播且不受橫波干擾,表明本文方法突破了該限制,具有更強適用性。

a、d—對應介質類型Ⅲ; b、e—對應介質類型Ⅳ; c、f—對應介質類型Ⅴ

圖5為純qP波方程分別在介質類型Ⅰ、Ⅱ、Ⅲ時400 ms時刻波場快照。對比VTI(圖5a)、HTI(圖5b)和TTI(圖5c)3種介質結果發現,通過引入坐標旋轉,使得TTI介質具有更強適用性,更能反映實際地下真實介質的情況。

a—介質類型Ⅰ; b—介質類型Ⅱ; c—介質類型Ⅲ

為了說明偽解析法的計算精度,設計了大小為2 000 m×2 000 m的均勻介質TTI模型,空間網格間距為5 m×5 m,時間步長1 ms,采用60 Hz主頻的Ricker子波,采用海綿邊界條件,縱波震源位于模型中心,模型參數為縱波波速為1 700 m/s,ε=0.25,δ=0.1,θ=45°,由于60 Hz主頻的Ricker子波最大頻率約為150 Hz,最大縱波速度為1 700 m/s,空間網格間距為5 m,因此一個波長內的采樣點數約為1700/150/5=2.2個,基本達到了Nyquist采樣定理極限精度的要求。

圖6為使用不同方法得到的350 ms時刻純qP波波場快照,其中6a為有限差分方法、6b為偽譜法和6c為偽解析法。對比這3種方法可以看到,有限差分方法出現明顯的頻散現象,這是有限差分在時間和空間方向的誤差造成的;而偽譜法雖然在空間方向上能夠得到Nyquist譜精度,但由于時間方向上的二階近似誤差,使得其仍出現了較為明顯的頻散;而圖6c采用偽解析方法在計算空間導數的同時,對時間二階近似誤差進行了補償,使得其能夠在時間和空間上同時達到Nyquist譜精度,說明了相比于有限差分方法和偽譜法,偽解析方法具有更高的計算精度。

a—有限差分法(FDM); b—偽譜法(PSM); c—偽解析法(PAM)

3.2 復雜BP模型

為檢驗所提方法對傾角劇變復雜介質的穩定性,截取一段傾角變化范圍是-55°～55°的二維BP2007聲學TTI介質模型進行測試。該模型規模為17 km×11.25 km,空間網格單元為6.25 m×6.25 m,時間步長為5 ms,總時間采樣點為1 600,記錄長度為8 s,縱波震源位于(68.5 km,0),采用主頻為20 Hz的Ricker子波震源,采用海綿邊界條件,接收排列范圍是0～17 km,排列深度為0,接收點間距為25 m,共681道接收,模型參數見圖7。

a—縱波速度v0;b—縱波各向異性參數ε;c—變異系數δ;d—極化角θ

圖8為采用二維BP2007模型得到的TTI介質純qP波4 s時刻波場快照及其地震記錄。從波場快照(圖8a)可見,地震波場能穩定地傳播,說明所提方法能適應傾角劇烈變化的情形,具有較高穩定性;同時,在地震記錄(圖8b)中,地震記錄和波場快照信噪比高,無頻散,波場反射信息清楚且豐富,繞射波明顯,反映出偽解析法能夠通過補償時間步長上的誤差有效地提高計算精度,壓制頻散。

圖8 二維BP2007聲學TTI介質4 s時刻波場快照(a)及地震記錄(b)Fig.8 Two-dimensional BP2007 acoustic TTI medium 4 s time wave field snapshot(a) and seismic record(b)

4 結論

本文在前人的研究基礎上,將qP波擬微分方程變換到空間—波數域,并通過坐標變換,推導了時間域TTI介質二階純qP波波動方程,引入偽解析算法,實現了基于偽解析法的TTI介質純qP波地震波場正演模擬。通過多種模型的測試,得出以下認識和結論:

1)克服了擬聲波方程的局限性,消除了偽橫波干擾,不受模型參數限制且地震波場能穩定傳播,具有較為廣泛的應用前景;

2)與其他方法相比(有限差分法及偽譜法),偽解析法能有效提高數值模擬精度;

3)簡單及復雜模型測試驗證了本文方法的正確性與適用性,能夠適應復雜介質,且具有更高的計算精度。