復雜近地表速度廣義近似反演方法研究

李 哲,程丹丹

(中國石油大學(華東)，山東青島266580)

隨著我國油氣勘探的重點目標由東部地區(qū)向西部和南方探區(qū)轉(zhuǎn)移，復雜近地表成為制約地震勘探效果的關鍵因素。在以勝利油田和大慶油田為主體的東部探區(qū)，近地表相對比較簡單，具有典型的低速帶、降速帶及正常速度帶3層結(jié)構(gòu)特征。此類近地表結(jié)構(gòu)中地震波速度向下逐層遞增，一般采用折射波法[1]來獲取近地表結(jié)構(gòu)。在復雜近地表探區(qū)，地形起伏劇烈，低速帶與降速帶不明顯，同時表層橫向巖性與速度變化快，常用的近地表調(diào)查和速度建模方法效果較差，不能滿足地震勘探的要求。

地震波旅行時層析是一種利用地震波在介質(zhì)中傳播時間來反演地震波所經(jīng)過路徑上介質(zhì)速度的方法。該方法對地下介質(zhì)的速度分層結(jié)構(gòu)沒有特殊要求，并能適用于橫向速度變化劇烈的地下介質(zhì)，因而成為最有效的近地表速度反演方法。國內(nèi)外眾多學者對旅行時層析反演方法進行了研究。自1984年第54屆SEG年會上設置了地震層析成像研究的專題之后，地震層析成像研究在地震勘探領域得以快速推進。經(jīng)過近30年的發(fā)展，國內(nèi)外關于地震層析成像方法的研究取得了豐碩的成果，主要有透射波層析[2-4]、反射波層析[5-8]、折射波層析[9-10]、井間層析[11-13]、多波聯(lián)合層析[14-16]、初至波層析[17-19]、波形層析[20-21]、起伏地表層析[22-23]等。這些工作也極大地推動了我國地震層析成像方法研究和實際應用的步伐。

對于復雜近地表探區(qū)，旅行時速度反演在地震勘探中主要有3個用途：一是靜校正量的計算，只要靜校正量合適，不要求速度值與地層厚度正確，相當于求解半確定性問題;二是用于觀測系統(tǒng)參數(shù)的確定，要求大尺度范圍內(nèi)一致，主要涉及震源激發(fā)深度的確定;三是偏移速度模型的建立，要求小尺度上的一致性。本文針對復雜近地表旅行時層析成像的數(shù)據(jù)特征，設計和推導了基于廣義逆的層析反演近似算法，并提出了4種針對性的約束方法。通過約束廣義近似旅行時層析速度反演建立合適的近地表速度模型，用于復雜近地表區(qū)地震資料的層析靜校正，恢復了炮集記錄上反射波同相軸的雙曲線形態(tài)。

1 復雜近地表探區(qū)的特點

1.1 復雜近地表探區(qū)的地質(zhì)特點

復雜近地表是指地形起伏大、速度變化迅速、巖性變化劇烈、厚度橫向差異大的淺表結(jié)構(gòu)。其主要地質(zhì)特點體現(xiàn)在以下方面。

1) 地表高程差別大。例如南方海相碳酸鹽巖探區(qū)，一般為大山區(qū)，地表高差往往在四五百米以上，鄂西、鎮(zhèn)巴和黑池梁等探區(qū)局部高程差可達1000多米[24];在西部黃土塬探區(qū)，地形起伏高差也有100多米。

2) 表層速度變化迅速。有的地表速度只有幾百米每秒，而不遠處就能達到2000m/s以上;在一些灰?guī)r出露地區(qū)，已經(jīng)不存在傳統(tǒng)意義上的低、降速帶，表層速度已經(jīng)達到3000～4000m/s。

3) 表層巖性變化劇烈。復雜近地表受到了大自然環(huán)境變化的影響，是大氣圈滲入巖石圈而形成的特殊地質(zhì)結(jié)構(gòu)，其巖性特征、層理結(jié)構(gòu)、膠結(jié)程度等受大自然環(huán)境作用而發(fā)生了改變，橫向均質(zhì)性極差。例如我國西部典型的山地、沙漠、戈壁、礫石覆蓋區(qū)、黃土覆蓋區(qū)等探區(qū)，橫向上巖性極不穩(wěn)定。正是這種橫向上巖性的變化，結(jié)合差異風化作用，使得低、降速帶厚度橫向差異明顯。

1.2 復雜近地表條件下地震資料特點

復雜近地表條件嚴重制約地震資料質(zhì)量，一方面是通過震源激發(fā)特性影響地震資料的幅頻特性和信噪比;另一方面是接收到的地震波先后兩次經(jīng)過復雜近地表傳播，其運動學與動力學特性受到明顯影響。由于復雜近地表介質(zhì)的非完全彈性性質(zhì)與橫向非均質(zhì)性，地震波發(fā)生吸收、衰減、散射和漫射等現(xiàn)象，不僅使地震波場變得復雜，也嚴重衰減有效波地震能量，降低了信噪比。

復雜近地表條件下，在共炮點道集上直達波同相軸不再是直線，反射波時距曲線也明顯偏離雙曲線特征(圖1)。這類地震記錄的初至波拾取存在2個問題：①存在極低信噪比的地震道，初至波無法識別與拾取，減少了可用于旅行時層析的數(shù)據(jù);②由于地震資料信噪比低，以及波形發(fā)生畸變等原因，拾取到的初至時間存在較大的誤差，這將嚴重影響表層旅行時層析速度反演的精度。因此，研究設計層析反演方法要考慮實際旅行時數(shù)據(jù)誤差。

圖1 起伏地表區(qū)典型炮集記錄

為此，在復雜近地表旅行時層析反演中，需要通過旅行時數(shù)據(jù)的統(tǒng)計分析確保用于速度反演的旅行時數(shù)據(jù)正確，并充分利用已知信息來約束速度反演，以彌補旅行時數(shù)據(jù)的不足。

2 復雜近地表旅行時層析廣義近似方法

旅行時層析成像的基本數(shù)學原理是Radon變換[25]。在旅行時層析中將地下介質(zhì)按一定大小進行網(wǎng)格化，用vj(j=1,2,…,M)表示第j個網(wǎng)格的速度。根據(jù)地震波傳播時間與速度、地震波傳播路徑的關系，得到的表達式為

(1)

式中：ti為地震波沿第i(i=1,2，…，N)條射線傳播的時間;lij為第i條射線在第j個網(wǎng)格中的長度。公式(1)中速度vj與旅行時ti之間的關系不是線性關系，用速度vj的倒數(shù)sj來替換公式(1)中的1/vj，從而得到用慢度表示的方程：

(2)

公式(2)形式上旅行時ti與未知量慢度sj已經(jīng)成為線性關系。但旅行時層析反演依然不是一個簡單的線性反演，因為這個反演中地震波傳播的路徑也不是已知的。實際上，這個反演不僅要反演出速度分布，而且必須同時反演出射線分布。由于地震波傳播路徑由速度分布與費馬原理共同決定，旅行時層析速度反演實際上是一個復雜的非線性反演問題。

2.1 廣義反演及其近似表達式

用列向量d=(t1,t2，…,tN)T表示有效的旅行時數(shù)據(jù)，列向量s=(s1,s2，…，sM)T表示需要求取的慢度;將每條射線在網(wǎng)格中的長度用下面的射線矩陣G來表示

(3)

則公式(2)變?yōu)?/p>

(4)

射線矩陣G不是方陣，公式(4)不能直接求解?？紤]到射線矩陣G與s的依賴關系，假設初始慢度s0，用于獲取射線矩陣G，并計算出地震波在慢度s0下傳播的時間d0，即

(5)

進而得到

(6)

式中：Δs表示s與s0的差值；Δd表示d與d0的差值。將公式(5)代入公式(6)得到

(7)

求解公式(7)得到修正的慢度。由于射線參數(shù)矩陣G是根據(jù)初始慢度s0計算的，與實際射線矩陣有差異，因而求取的慢度和實際慢度有差異。為此，需要用修正后的慢度作為初始慢度，重新進行計算。這樣就形成了旅行時層析速度反演的迭代算法。

由于射線矩陣G不是正定矩陣，公式(7)不能直接求解?？梢岳米钚《朔╗27]得到慢度修正量Δs的廣義反演形式：

(8)

公式(8)的求解需要計算M×M階的逆矩陣，M為慢度修正量的未知數(shù)數(shù)目，通常比較大。公式(8)難以直接計算，同時參數(shù)矩陣G與實際射線矩陣有差異，公式(8)精確計算的結(jié)果肯定是錯誤的，為此考慮近似求解。

用J表示正定矩陣GTG，其元素Jkm是經(jīng)過第k網(wǎng)格和第m網(wǎng)格的射線長度數(shù)據(jù)的互相關，即

(9)

公式(9)右邊被求和項只有當一條射線同時通過第k網(wǎng)格和第m網(wǎng)格才不為零，其它項均為零。即：當k≠m時，Jkm接近于零;當k=m時，被求和項為網(wǎng)格的自相關，當一個網(wǎng)格沒有射線經(jīng)過時，Jkk=0。剔除沒有射線經(jīng)過的網(wǎng)格，則可以得到一個主對角線元素絕對占優(yōu)的對角陣J。因此反演中只考慮J的主對角線元素，得到公式(8)的近似公式

(10)

公式(10)不需要求解大型矩陣，在近似計算中只需計算主對角線元素和射線矩陣的乘法，因此算法是穩(wěn)定的。

2.2 復雜近地表初始模型建立

建立初始模型是旅行時層析反演的關鍵。常用的建模方法有2種。

1) 建立與已知信息相一致的初始模型。在地震勘探中往往利用微測井、小折射和巖性探測技術(shù)等獲取探區(qū)內(nèi)某些關鍵點上的實測表層速度與巖性信息。利用這些已知點的地層速度分布，結(jié)合地形起伏情況，可以建立起一個比較接近實際的速度模型?；谠摮跏寄Ｐ?，只需要反演出較小的慢度修正量。

2) 按旅行時層析反演算法的要求來建立初始模型。旅行時反演中只有那些有射線經(jīng)過的網(wǎng)格才能被反演出來。對于近地表來說，如果一開始建立的模型淺層速度比深層速度大，那么地震波射線就只會在淺表，下部介質(zhì)里沒有射線，自然就反演不出來。因此，應建立速度隨深度逐漸增加的模型，讓射線能夠下傳到下部介質(zhì)，符合大多數(shù)情況，即表層速度基本結(jié)構(gòu)是速度隨深度增加的。

無論采用哪種建模方法，都應在反演中利用微測井、小折射和巖性探測技術(shù)等提供的實測速度信息作為反演的約束條件，以減少反演的多解性。

2.3 旅行時拾取誤差的影響

第一個地震波到達接收點，接收點開始振動的時刻為起跳時間，這個時間就是旅行時反演所需要的初至時間。由于存在背景噪聲，當?shù)卣鸩ㄒ鸬恼駝有∮诒尘霸肼晻r，無法確定哪個時刻是起跳時間。通常在技術(shù)上總是拾取地震記錄上第一個能量明顯高于背景振動的到達時間，作為初至波時間。

在復雜近地表條件下，地震記錄變得十分復雜，無論是直達波還是反射波，其同相軸均偏離了正常的形態(tài)且伴隨著低信噪比，進一步增加了旅行時拾取的難度和可信度。

在近地表旅行時層析反演中炮檢距相對較小，因而旅行時拾取的相關誤差很大。由于討論相對誤差對速度反演的影響比較困難，假設所有旅行時均具有相同的絕對誤差σ，則根據(jù)公式(10)進行反演，估算慢度誤差列向量為

(11)

公式(11)是旅行時數(shù)據(jù)誤差傳遞方程。假設射線在每個網(wǎng)格中的長度相同，則速度反演結(jié)果將被放大N/L倍(L為射線的平均長度)。因此，反演中應該考慮旅行時數(shù)據(jù)的誤差，避免將旅行時誤差反演到速度模型中。

2.4 低、降速帶的確定

旅行時層析反演獲得近地表速度分布后，需要進一步確定低、降速帶。根據(jù)公式(10)可知，沒有射線經(jīng)過的網(wǎng)格不參加反演，即這些網(wǎng)格的速度不正確。實際應用時首先根據(jù)射線分布確定速度反演的可靠范圍，在該范圍內(nèi)分析低、降速帶的分布。

利用微測井、小折射以及巖性探測等技術(shù)提供的局部近地表可靠速度分層信息，對近地表旅行時層析反演資料進行一致性檢測，進而標定速度范圍，同時結(jié)合層析射線的密度進行低、降速帶劃分。

3 實際地震資料試處理

對南方海相碳酸鹽巖某復雜地表區(qū)實際二維地震資料進行旅行時層析廣義近似反演試處理。該區(qū)地形起伏達到500m，巖石出露地表，最低速度在1400m/s左右，屬于典型的復雜近地表結(jié)構(gòu)。根據(jù)工區(qū)的基本情況，考慮到地表速度較高，建立如圖2 所示的速度隨深度逐漸增加的初始模型。

對該二維測線的所有炮集數(shù)據(jù)進行初至波拾取，并繪制初至波旅行時與炮檢距的交會圖，結(jié)果如圖3a所示。觀察圖3a可知，炮檢距在4000m以上時存在散點，即可能存在誤差。近地表初始模型的深度為1000m，而大炮檢距的地震波向下的穿透深度超過了1000m，為此切除了炮檢距大于2500m的初至時間，結(jié)果如圖3b所示。

圖2 南方復雜地表區(qū)實際二維測線初始速度模型

圖3 實際二維地震資料拾取的初至波旅行時隨炮檢距的分布a 所有初至波旅行時; b 小炮檢距初至波旅行時

圖4為直接利用圖3b所示小炮檢距初至波旅行時數(shù)據(jù)進行無約束廣義近似層析反演給出的射線分布和速度分布，可以看出剖面兩邊由于邊界原因，速度與初始模型接近，反演不充分;中部射線比較松散，淺表射線比較集中;從速度分布上看，沒有清晰地反演出低、降速帶。

圖4 無約束廣義近似層析反演的射線分布(a)與速度分布(b)

通過對實際資料的分析，根據(jù)工區(qū)特點引入4種約束條件：①速度區(qū)間的控制，根據(jù)已知的先驗信息設定一個速度區(qū)間，防止速度過小和過大;②利用小折射和微測井資料提供的先驗信息進行約束;③速度奇異點約束，利用中值濾波去除速度奇異值;④利用附近網(wǎng)格修正量均值修正無射線通過網(wǎng)格。圖5為仍然采用圖3b中的初至波旅行時數(shù)據(jù)進行增加約束后的廣義近似層析反演給出的射線分布與速度模型。從射線分布(圖5a)上看，對低、降速帶有較好的覆蓋，近地表速度模型(圖5b)的層狀結(jié)構(gòu)非常明顯。

利用約束廣義近似層析反演出的合理的近地表速度模型(圖5)，以3000m/s為低降速帶門檻值，得到低、降速帶底界，進而計算靜校正量，對實際二維地震資料的炮集記錄進行靜校正處理。圖6a 為原始炮集記錄，直達波和反射波同相軸畸變明顯;圖6b是應用靜校正量后的結(jié)果，反射波同相軸基本上恢復了雙曲線的特性。

圖5 約束廣義近似層析反演的射線分布(a)與速度分布(b)

圖6 復雜地表區(qū)實際二維地震資料單炮記錄靜校正前(a)、后(b)對比

4 結(jié)束語

復雜近地表條件下地震波場復雜，地震資料信噪比低，拾取的旅行時誤差比較大。我們提出的廣義近似層析反演方法是一種考慮旅行時拾取誤差的層析成像方法，能夠獲得較穩(wěn)定的近地表速度模型。針對旅行時誤差與速度模型誤差的關系，推導了慢度反演誤差的計算公式，可進一步研究將其應用于反演精度控制的方法。

針對復雜地表區(qū)實際地震資料的無約束廣義近似層析反演結(jié)果存在對近地表的分層不明顯、射線完全集中在淺表等問題，研究并提出了4種約束方法，使得約束后反演的射線能較好地覆蓋低降速帶，速度分布清晰地展示了近地表地層結(jié)構(gòu)。將約束廣義近似層析反演出的速度模型用于復雜地表區(qū)實際地震資料的靜校正，恢復了反射波同相軸的雙曲線形態(tài)。

致謝:中國石油大學(華東)地球科學與技術(shù)學院地球物理系莫言工作組為本方法的研究提供了所有資料，課題組老師和同學對本文的研究提出了寶貴意見，在此表示衷心感謝！

參 考 文 獻

[1] 徐凌,崔興福,齊莉,等.復雜近地表層析反演與波場延拓聯(lián)合基準面校正[J].石油物探,2007,46(3)：226-230

Xu L,Cui X F,Qi L,et al.Combined datum correction of tomography inversion with wavefield continuation under complex near surface condition[J].Geophysical Prospecting for Petroleum,2007,46(3):226-230

[2] Bording R P,Gersztenkorn A,Lines L R,et al.Applications of seismic travel-time tomography[J].Geophysical Journal of the Royal Astronomical Society,1987,90(2):285-303

[3] Gelbke C,Miranda F,Sattel G.Results of a seismic transmission tomography survey at the Grimsel Rock Laboratory[J].Log Analysist,1989,30(4):243-260

[4] 張建中,丁興號.一種2.5維井間透射波層析成像方法[J].聲學學報,2007,32(1):91-96

Zhang J Z,Ding X H.An approach to 2.5-D crosswell seismic transmission tomography[J].Acta Acustica,2007,32(1):91-96

[5] Zhou C,Crawley S,Whitmore D,et al.Seismic reflection tomography with 3D RTM angle gathers[EB/OL].(2013-06-10)[2013-12-30].http:∥earthdoc.eage.org/publication/publicationdetails/?publication=69508

[6] Stork C.Reflection tomography in the postmigrated domain[J].Geophysics,1992,57(5):680-692

[7] 井西利,楊長春,王世清.一種改進的地震反射層析成像方法[J].地球物理學報,2007,50(6):1831-1836

Jing X L,Yang C C,Wang S Q.A improved seismic reflection tomographic method[J].Chinese Journal of Geophysics,2007,50(6):1831-1836

[8] Brzostowski M A,McMechan G A.3-D tomographic imaging of near-surface seismic velocity and attenuation[J].Geophysics,1992,57(3):396-403

[9] Mari J L O,Mendes M.3D Near surface imaging thanks to refraction tomography and delay method[EB/OL].(2010-06-14)[2013-12-30].http:∥www.earthdoc.org/publication/publicationdetails/?publication=39624

[10] Afanasiev M,Pratt R G,Kamei R,et al.Crosshole Tomography,Waveform Inversion and Anisotropy-A Combined Approach Making Use of Simulated Annealing[EB/OL].(2012-06-04)[2013-12-30].http:∥www.earthdoc.org/publication/publicationdetails/?publication=59614

[11] 左建軍,林松輝,孔慶豐,等.井間地震直達波和反射波聯(lián)合層析成像及應用[J].石油地球物理勘探,2011,46(2):226-231

Zuo J J,Lin S H,Kong Q F,et al.Crosswell seismic direct wave and reflected wave joint tomographic imaging and Application[J].Oil Geophysical Prospecting,2011,46(2):226-231

[12] 瞿辰,楊文采,于常青.高分辨率過井地震層析成像的級聯(lián)算法及其應用[J].物探與化探,2010,34(6):814-820

Huo C,Yang W C,Yu C Q.High resolution over the cascade algorithm well seismic tomography and its application[J].Geophysical and Geochemical Exploration,2010,34(6):814-820

[13] 李培明,閆智慧,郭明杰,等.聯(lián)合應用初至和淺層反射的二維可形變層析靜校正[J].石油地球物理勘探,2010,45(5):647-654

Li P M,Yan Z H,Guo M J.The combined application of first break and shallow reflection of 2D deformable tomographic static correction[J].Oil Geophysical Prospecting,2010,45(5):647-654

[14] 張樹林,朱介壽,賀振華.井間地震和逆VSP 聯(lián)合層析成像[J].石油地球物理勘探,1993,28(5):577-583

Zhang S L,Zhu J S,He Z H.Joint tomography of interborehole seismic data and inverse VSP data[J].Oil Geophysical Prospecting,1993,28(5):577-583

[15] Meléndez A,Korenaga J,Sallarǒs V,et al.TOMO3D-A new 3-d joint refraction and reflection travel-time tomography code for active-source seismic data[EB/OL].(2013-06-04)[2013-12-30].http:∥www.researchgate.net/publication/258621526_tomo3d_a_new_3-D_joint_refraction_and_reflection_travel-time_tomography_code_for_active-source_seismic_data

[16] Noble M,Thierry P,Taillandier C,et al.High-performance 3D first-arrival traveltime tomography[J].The Leading Edge,2010,29(1):86-93

[17] Li S,Fomel S,Vladimirsky A.Prestack first-break traveltime tomography using the double-square-root eikonal equation[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012,1-6

[18] 袁剛,馮心遠,蔣波,等.約束層析反演及其在地震速度計算中的應用[J].石油物探,2013,52(1):55-59

Yuan G,Feng X Y,Jiang B,et al.Constrained tomography inversion and its application in seismic velocity computation[J].Geophysical Prospecting for Petroleum,2013,52(1):55-59

[19] Smithyman B,Pratt R G,Hayles J,et al.Detecting near-surface objects with seismic waveform tomography[J].Geophysics,2009,74(6):WCC119-WCC127

[20] Sheng J,Leeds A,Buddensiek M,et al.Early arrival waveform tomography on near-surface refraction data[J].Geophysics,2006,71(4):U47-U57

[21] Zhang K,Li Z C,Zeng T S,et al.Tomographic velocity inversion for ADCIGs in areas with a rugged surface[J].Applied Geophysics,2012,9(3):313-318

[22] 張兵,徐兆濤,王華忠,等.山前帶地震數(shù)據(jù)共成像點道集層析速度反演建模方法研究[J].石油物探,2012,51(6):590-597

Zhang B,Xu Z T,Wang H Z,et al.Common imaging gathers tomography velocity inversion and model building in foothill area[J].Geophysical Prospecting for Petroleum,2012,51(6):590-597

[23] 董良國,郭曉玲,吳曉豐,等.起伏地表彈性波傳播有限差分法數(shù)值模擬[J].天然氣工業(yè),2007,27(10):38-41

Dong L G,Guo X L,Wu X F,et al.Finite difference numerical simulation for the elastic wave propagation in rugged topography[J].Natural Gas Industry,2007,27(10):38-41

[24] Helgason S.The Radon transform[M].2nd ed.USA:Birkhauser Boston Inc,1999:1-200

[26] Fitzgibbon A,Pilu M,Fisher R B.Direct least square fitting of ellipses[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,1999,21(5):476-480