基于快速步進(jìn)法的三維反射波追蹤

姜莉

摘 要： 引入了一種改進(jìn)的快速步進(jìn)算法，對三維復(fù)雜介質(zhì)中反射波進(jìn)行波前追蹤，并利用QT和OpenGL相結(jié)合的方法對速度場以及射線路徑進(jìn)行可視化模擬。通過模型仿真，驗(yàn)證了該方法的效率、精度以及可靠性。

關(guān)鍵詞： 快速步進(jìn)法； 反射波； 射線追蹤； 波前追蹤

中圖分類號： TN911.7?34； P314.3 文獻(xiàn)標(biāo)識碼： A 文章編號： 1004?373X（2014）11?0132?03

Abstract： In the paper， an improved fast stepping algorithm is introduced to make wavefront tracing of the reflected wave in 3D heterogeneous media， and perform the visualization simulation of velocity field and ray path by using QT and OpenGL. The efficiency， accuracy and reliability of the newalgorithm were verified by model simulation.

Keywords： fast stepping method； reflected wave； ray tracing； wavefront tracing

費(fèi)馬原理指出，地震波沿著所需時間最短的路徑傳播，在高頻近似條件下，地震波場的主能量沿射線軌跡傳播。地震射線追蹤的目的就是尋找震源到各個接收點(diǎn)的地震波傳播的最短傳播路徑，從而研究地震波的傳播規(guī)律，對地層探測及層析成像具有重要的意義。針對傳統(tǒng)射線追蹤方法僅能求出地震波的初至旅行時的缺點(diǎn)，本文在快速步進(jìn)算法的基礎(chǔ)上給出一種多級的方法，對反射波進(jìn)行射線追蹤。該方法將波前到達(dá)的每一層分割為一個獨(dú)立計算的區(qū)域，把波前到達(dá)的每層界面重新作為波源重新計算，從而逐層推進(jìn)式地追蹤反射波或透射波的波前。

1 算法實(shí)現(xiàn)

快速步進(jìn)算法的主要思想就是通過系統(tǒng)地用上風(fēng)格式構(gòu)建旅行時場來解程函方程。程函方程的上風(fēng)差分格式實(shí)際上意味著信息傳播是從走時的小值到大值。因此，快速步進(jìn)算法重點(diǎn)在通過從[T（x，y，z）]的小值向外擴(kuò)展來構(gòu)建方程的解[1]。

1.1 旅行時場方程

在3D旅行時計算時，程函方程為：

1.2 模型建立

把目標(biāo)區(qū)域離散為網(wǎng)格，將沒被界面穿過的網(wǎng)格稱作規(guī)則網(wǎng)格，每個網(wǎng)格頂點(diǎn)稱作規(guī)則結(jié)點(diǎn)，如圖1中大實(shí)心點(diǎn)所示。在規(guī)則網(wǎng)格部分，每一個網(wǎng)格都是具有八個網(wǎng)格頂點(diǎn)的區(qū)域，同時每個規(guī)則結(jié)點(diǎn)屬于八個網(wǎng)格。如果規(guī)則網(wǎng)格被界面切割，稱作不規(guī)則網(wǎng)格，網(wǎng)格邊與界面相交，產(chǎn)生的交點(diǎn)稱作不規(guī)則網(wǎng)格結(jié)點(diǎn)，如圖1中小實(shí)心點(diǎn)所示，每個交點(diǎn)屬于四個網(wǎng)格。如果要更新某一點(diǎn)的走時，認(rèn)為每個網(wǎng)格相互獨(dú)立。

1.3 算法實(shí)現(xiàn)

如圖2所示，為了加快計算，把目標(biāo)區(qū)域的所有網(wǎng)格結(jié)點(diǎn)劃分為三個子集，首先把具有初始旅行時值的點(diǎn)記為確認(rèn)值點(diǎn)（實(shí)心黑點(diǎn)），這些點(diǎn)的旅行時在計算中不再改變；然后把與確認(rèn)值點(diǎn)相鄰一個網(wǎng)格的點(diǎn)記為近點(diǎn)（實(shí)心綠點(diǎn)），近點(diǎn)需要計算一個旅行時的值，但它不一定是最小值，是可改變的；其他點(diǎn)記為遠(yuǎn)點(diǎn)（空心點(diǎn)），遠(yuǎn)點(diǎn)的旅行時值還沒有計算，需要逐次推進(jìn)計算。

計算最小旅行時值的過程如下：首先在近點(diǎn)中搜索旅行時最小的點(diǎn)，把這個點(diǎn)從近點(diǎn)子集中去除，加入確認(rèn)值點(diǎn)子集，然后把此點(diǎn)的所有非確認(rèn)值的臨近點(diǎn)歸入近點(diǎn)子集，如果臨近點(diǎn)屬于遠(yuǎn)點(diǎn)子集，則從中去除；接著重新計算該點(diǎn)臨近點(diǎn)的旅行時值，如果比原值大，則保留原值，如果小，則取而代之。重復(fù)前面的步驟，直到確認(rèn)值點(diǎn)集包含整個目標(biāo)區(qū)域。本算法利用堆排序技術(shù)，可有效地在近點(diǎn)集合中搜索最小元素，使算法加快速度。

1.4 反射波前追蹤

對反射波射線追蹤時如果直接利用快速步進(jìn)法計算時，會存在以下問題：使用快速步進(jìn)方法的有限差分解必須服從絕對初至?xí)r間。對于反射波，波前至少穿過計算區(qū)域兩次，因此只能通過多級方法來解決。多級方法的原理如圖3所示。

第一步，波前傳播到反射界面，然后將該區(qū)域進(jìn)行復(fù)制。

第二步，將快速步進(jìn)方法用于復(fù)制的區(qū)域。第二步開始時，窄帶圍繞著界面，這樣做并不破壞因果性，因?yàn)槌踔练瓷洳ㄇ爸坏竭_(dá)每個界面結(jié)點(diǎn)一次。

1.5 走時更新

1.5.1 規(guī)則網(wǎng)格

對于規(guī)則網(wǎng)格點(diǎn)的走時更新，用上風(fēng)格式近似表示式（1）的梯度函數(shù)，可得到如下的表達(dá)式：

1.5.2 不規(guī)則網(wǎng)格

如果被更新結(jié)點(diǎn)處于不規(guī)則網(wǎng)格中，使用相鄰的不規(guī)則點(diǎn)的走時來更新某一結(jié)點(diǎn)時，假設(shè)波前近似為平面。當(dāng)使用不規(guī)則網(wǎng)格中的確認(rèn)值點(diǎn)更新時，此網(wǎng)格單元中的確認(rèn)值點(diǎn)數(shù)可能是一個，兩個，三個或大于三個。如果由三個確認(rèn)值點(diǎn)的信息來更新此點(diǎn)走時，這三點(diǎn)可以惟一確定一個波陣面，準(zhǔn)確計算該點(diǎn)走時。如果僅有一個或兩個點(diǎn)是確認(rèn)值點(diǎn)，必須增加額外的約束信息，找到波陣面的解。當(dāng)確認(rèn)值點(diǎn)多于三個時，將產(chǎn)生冗余信息。下面分別考慮這四種情況。

（1） 由三個確認(rèn)值點(diǎn)更新

[x，][y，][z]分別為[x，][y]和[z]三個方向的單位矢量。由以上的計算公式可以看出，每個網(wǎng)格的慢度矢量在計算旅行時的過程中，就可求出，不用特別計算，從而可節(jié)省計算量。根據(jù)射線理論，射線與慢度矢量的方向一致，射線在模型中某一點(diǎn)的傳播方向，用該點(diǎn)周圍網(wǎng)格點(diǎn)的慢度矢量插值求出，在本文中采用雙線性插值。

2 實(shí)驗(yàn)及仿真

本文采用層狀介質(zhì)模型進(jìn)行實(shí)驗(yàn)仿真，模型為一個兩層的三維介質(zhì)模型。如圖4所示，上層用單一顏色表示其傳播速度為常數(shù)，下層的速度隨深度變化，顏色越深表示傳播速度越快。目標(biāo)區(qū)域長寬深分別為2 000 m，2 000 m和1 000 m，在[x，][y]和[z]方向分別劃分為81，81和41個網(wǎng)格，網(wǎng)格間距在各個方向上都是25 m，炮點(diǎn)位于地面上（0，1 000，1 000）坐標(biāo)處。圖4是模型速度場的剖面圖，拖動右邊滑塊可以控制各方向剖面的移動。圖5是射線及界面模擬圖，拖動滑塊可以對圖像進(jìn)行旋轉(zhuǎn)。從圖中可以看出，模擬結(jié)果與理論結(jié)果一致，快速步進(jìn)法確實(shí)是一種快速穩(wěn)定，適用于實(shí)際的方法。

3 結(jié) 論

本文采用Fortran和QT類庫以及OpenGL相結(jié)合，的方法很好地實(shí)現(xiàn)了快速步進(jìn)法射線追蹤，不僅利用前人已有的成果進(jìn)一步降低開發(fā)難度，提高開發(fā)效率。同時易于向Unix移植，以適應(yīng)一些用戶的特殊需要。由模型實(shí)驗(yàn)看，所得結(jié)果比較客觀地反映出了實(shí)際情況。所以，此法是一種很好的三維正演方法。

參考文獻(xiàn)

[1] KIM Seongjai. 3?Deikonalsolvers： first?arrival traveltimes [J]. Geophysical， 2002， 67（4）： 1225?1231.

[2] 許琨，吳律，王妙月.改進(jìn)Moser法射線追蹤[J].地球物理學(xué)進(jìn)展，1998（4）：60?64.

[3] 張云姝，顧漢明，師學(xué)明.基于MOSER曲射線追蹤的SIRT聲波層析成像[J].工程地球物理學(xué)報，2005（3）：167?169.

[4] 張霖斌，劉迎曦，趙振峰，等.有限差分法射線追蹤[J].石油地球物理勘探，1993，42（6）：198?200.

[5] 張美根，賈豫葛，王妙月，等.界面二次源波前擴(kuò)展法全局最小走時射線追蹤技術(shù)[J].地球物理學(xué)報，2006，49（4）：1169?1175.

[6] 李艷.地震勘探中二維波動方程的小波數(shù)值正演模擬[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2002.

[7] 周竹生，劉喜亮，熊孝雨.彈性介質(zhì)中瑞雷面波有限差分法正演模擬[J].地球物理學(xué)報，2007，50（2）：567?573.

[8] 劉恩儒，岳建華，劉彥.具有離散裂縫空間分布的二維固體中地震波傳播的有限差分模擬[J].地球物理學(xué)報，2006，49（1）：180?188.

[9] 朱亞平，楊慧珠.OpenGL技術(shù)在地震數(shù)據(jù)可視化中的應(yīng)用[J].石油地球物理勘探，2000，35（4）：403?414.

[10] WOO Mason， NEIDER Jaekie， DAVIS Tom， et al. OpenGL編程權(quán)威指南[M].北京：中國電力出版社，2001.

[11] 王博樺，陳振一，王關(guān)橋，等.地震信號檢測新機(jī)理研究[J].現(xiàn)代電子技術(shù)，2013，36（21）：111?112.