基于程函方程的克?；舴蛉矫嬗蚱瞥上穹椒?/h1>
2020-07-26 05:25智敏朱建剛
物探與化探 2020年4期
智敏,朱建剛
(中煤科工集團(tuán)西安研究院有限公司,陜西 西安 710077)
0 引言
克希霍夫積分偏移是一種基于波動(dòng)方程積分解的偏移成像方法，與其他偏移方法相比，該方法能夠適應(yīng)非規(guī)則觀測(cè)系統(tǒng)，具有較高的計(jì)算效率[1]。目前主流處理系統(tǒng)的克?；舴蚱颇K主要針對(duì)于地面地震數(shù)據(jù)成像，要求炮檢點(diǎn)均位于地表，波場(chǎng)信息來自地下半空間，屬于半空間成像方法。對(duì)于某些特殊地震觀測(cè)系統(tǒng)，如煤礦井下槽波觀測(cè)系統(tǒng)，地震信號(hào)可能來自以炮檢點(diǎn)為焦點(diǎn)的完整橢圓上，屬于全平面域偏移成像，主流處理軟件并不適用。王季、姬廣忠等基于直射線層析重構(gòu)的槽波速度場(chǎng)，采用直射線克?；舴蚍e分偏移方法對(duì)井下反射槽波數(shù)據(jù)進(jìn)行了偏移成像，取得了一定的效果[2-3]，該方法與地面地震直射線偏移成像一樣，將繞射點(diǎn)與炮檢點(diǎn)之間的射線路徑視為直線，沒有考慮速度場(chǎng)彎化引起的射線彎曲及旅行時(shí)變化，造成成像點(diǎn)與原繞射點(diǎn)不重合，降低了成像分辨率。
針對(duì)直射線層析成像速度的這一缺陷，本文采用彎曲射線層析成像方法反演成像區(qū)域速度。對(duì)于彎曲射線層析成像，射線追蹤的精度直接影響著成像精度。常用的射線追蹤方法有試射法、旅行時(shí)線性插值法(LTI)、快速步進(jìn)法(FMM)、有限差分求解程函方程法。試射法對(duì)于簡(jiǎn)單速度模型具有較高精度，但不能處理折射波問題，對(duì)于復(fù)雜構(gòu)造常存在追蹤盲區(qū)，并且計(jì)算效率低下[4]。LTI、FMM算法及有限差分法均首先計(jì)算從震源傳至各節(jié)點(diǎn)的最小旅行時(shí)，并基于旅行時(shí)間采用反向追蹤技術(shù)求取各檢波點(diǎn)到炮點(diǎn)的射線路徑。在計(jì)算最小旅行時(shí)時(shí)，LTI算法常采用按列(行)掃描的方式計(jì)算全局最小旅行時(shí)[5]，該方法未能考慮來自各個(gè)方向的射線，對(duì)于橫向速度變化較大的介質(zhì)，計(jì)算結(jié)果往往并非全局最小，葉佩等采用了全方位循環(huán)法求得了網(wǎng)格節(jié)點(diǎn)的全局最小旅行時(shí)[6]，張建中等采用Dijkstra算法逐點(diǎn)計(jì)算了各節(jié)點(diǎn)的最小旅行時(shí)[7]；FMM算法通過上風(fēng)差分近似和窄帶技術(shù)近似模擬地震波前的傳播，具有計(jì)算精度高，運(yùn)算速度快、無條件穩(wěn)定的優(yōu)點(diǎn)，可以計(jì)算網(wǎng)格節(jié)的全局最小旅行時(shí)，但該方法近源區(qū)域網(wǎng)格節(jié)點(diǎn)旅行計(jì)算誤差較大，并對(duì)后續(xù)節(jié)點(diǎn)的計(jì)算精度有影響，需要對(duì)近源區(qū)的網(wǎng)格進(jìn)行行加密處理，以提高近源區(qū)旅行時(shí)的計(jì)算精度[8]；有限差分求解程函方程法是一種計(jì)算精度高、運(yùn)算量小的旅行時(shí)計(jì)算方法[9]，對(duì)于速度變化大的介質(zhì)，采用傳統(tǒng)的波前矩形擴(kuò)展方可能因違返波傳播的因果關(guān)系導(dǎo)致計(jì)算不穩(wěn)定，潘紀(jì)順等通過在算法中額外考慮首波、體波和散射波對(duì)旅行時(shí)影響的方式有效地解決了這一問題[10]。本文采用有限差分法求解程函方程的方法求取成像平面內(nèi)各節(jié)點(diǎn)的初至波旅行時(shí)場(chǎng)，算法上借鑒了潘紀(jì)順的處理方法，波前擴(kuò)展方法采用Disjkstra算法逐點(diǎn)計(jì)算各節(jié)點(diǎn)的全局最小旅行，射線反向追蹤采用了Aldridge提出的最速下降法求取射線路徑的方法[11]，并基于全局最小旅行時(shí)對(duì)二維復(fù)雜速度模型進(jìn)行全平面域克?；舴蚩臻g域偏移成像，通過對(duì)偏移結(jié)果的對(duì)比分析，認(rèn)為該方法可以較真實(shí)地反映異常速度體的形態(tài)大小和位置，成像結(jié)果分辨率較高。
1 方法原理
與常規(guī)地面地震克?；舴蚍e分法疊前深度偏移類似，全平面域克?；舴蚍e分疊前偏移成像方法由偏移速度求取和克?；舴蚍e分偏移兩部分組成，兩者在實(shí)現(xiàn)過程中均涉及成像區(qū)域旅行時(shí)場(chǎng)計(jì)算，本文采用有限差分求解程函方程的方法計(jì)算節(jié)點(diǎn)旅行時(shí)。
1.1 有限差分程函方程初至旅行時(shí)計(jì)算
旅行時(shí)計(jì)算是積分法偏移成像的重要組成部分，其計(jì)算精度直接影響著成像效果，在高頻近似下，平面波波前的傳播符合程函方程式(1)[9]：
(1)
式中,t為初至?xí)r間；s(x,y)為平面模型的慢度場(chǎng)。
采用有限差分近似方法求解式(1)，將速度模型離散成規(guī)則的矩形網(wǎng)格(圖1所示)，X、Y方向網(wǎng)格邊長(zhǎng)分別為Δx、Δy，將每個(gè)網(wǎng)格單元內(nèi)的慢度視為常數(shù)，將方程中的兩個(gè)時(shí)間微分算子用有限差分算子近似，如式(2)～(3)所示[9,12]:
圖1 矩形網(wǎng)格剖分示意Fig.1 Diagram of rectangular mesh partition
(2)
(3)
將式(2)、(3)代入式(1)，經(jīng)推導(dǎo)可以得到式(4)[9]:
(4)
式中,A=Δx2-Δz2，B=2ΔxΔz，C=Δx2+Δz2，已知矩形網(wǎng)格3個(gè)節(jié)點(diǎn)旅行時(shí)t0、t1、t2可以利用式(4)得到旅行時(shí)t3，計(jì)算中常常會(huì)出現(xiàn)根號(hào)下為負(fù)值而出現(xiàn)計(jì)算不穩(wěn)定的情況，此時(shí)可將節(jié)點(diǎn)0、1、2視為散射源，按直射傳播方式，分別計(jì)算節(jié)點(diǎn)3的初至?xí)r間，結(jié)合節(jié)點(diǎn)3的已有旅行時(shí)間，取5個(gè)值里面的最小值作為節(jié)點(diǎn)3的全局最小旅行時(shí)。
在旅行時(shí)場(chǎng)外推過程中，筆者采用Dijkstra算法逐點(diǎn)計(jì)算所有節(jié)點(diǎn)的全局最小旅行時(shí)，將所有計(jì)算節(jié)點(diǎn)分成3個(gè)子集，分別用P、Q、R表示，其中P為尚未計(jì)算旅行時(shí)的節(jié)點(diǎn)集合(節(jié)點(diǎn)旅行為無窮大)，Q為已計(jì)算旅行時(shí)的節(jié)點(diǎn)集合(不確定是否為全局最小旅行時(shí))，R為已確定為全局最小旅行時(shí)節(jié)點(diǎn)的集合，采用如下方法對(duì)各節(jié)點(diǎn)的集合屬性進(jìn)行更新：
1)將所有節(jié)點(diǎn)置于P集合，且所有節(jié)點(diǎn)旅行時(shí)均置為無窮大，計(jì)算炮點(diǎn)所在網(wǎng)格4個(gè)節(jié)點(diǎn)的旅行時(shí)，將其置于Q集；
2)在Q集中選取旅行時(shí)最小的節(jié)點(diǎn)作為次級(jí)震源，并將其置于R集，計(jì)算該節(jié)點(diǎn)周圍8個(gè)節(jié)點(diǎn)的旅行時(shí)，若待計(jì)算節(jié)點(diǎn)屬于P集合，則將其置于Q集合，更新該節(jié)點(diǎn)旅行時(shí)；若節(jié)點(diǎn)屬于Q集合，則將計(jì)算結(jié)果與先前結(jié)果比較，取小值作為該節(jié)點(diǎn)旅行時(shí)；若節(jié)點(diǎn)屬于R集，則不改變?cè)担?/p>
3)重復(fù)步驟2)直至所有節(jié)點(diǎn)均屬于R集為止，即可求得成像區(qū)域內(nèi)所有節(jié)點(diǎn)的全局最小旅行時(shí)。
1.2 初至波層析成像反演速度場(chǎng)
偏移速度場(chǎng)的計(jì)算精度直接影響著偏移成像效果，是克?；舴蚍e分法偏移成像的關(guān)鍵，對(duì)于非地面地震觀測(cè)系統(tǒng)，采用地面地震觀測(cè)系統(tǒng)的偏移速度分析方法獲取偏移速度比較困難，通常采用初至波層析反演方法獲得偏移速度，其中直射線層析成像方法應(yīng)用比較普遍，如文獻(xiàn)[2-3]中的井下槽波繞射偏移成像，當(dāng)反演區(qū)域速度差異較小時(shí)，可近似認(rèn)為射線在介質(zhì)中沿直射傳播，如圖2所示，炮檢點(diǎn)間的射線路徑為直線，ai,j為第i條射線在第j個(gè)網(wǎng)格中的射線長(zhǎng)度，Sj為第j個(gè)網(wǎng)格的慢度，ti為第i條射線的初至?xí)r，利用直射線層析反演方法可以獲得較可靠的成像速度，當(dāng)成像區(qū)域速度變化較大時(shí)，由于實(shí)際射線路徑常繞過低速區(qū)，并向高速區(qū)集中，射線路徑發(fā)生較大彎曲，不滿足直射線假設(shè)條件，常導(dǎo)致高速異常體偏大，低速異常體偏小[13]。
圖2 離散網(wǎng)格中直射線層析反演示意Fig.2 Sketch of straight-ray tomography inversion in regular grids
為了克服直射線層析的這一缺陷，將成像區(qū)域進(jìn)行網(wǎng)格剖分后，采用有限差分法求得各成像網(wǎng)格節(jié)點(diǎn)的全局最小旅行時(shí)，應(yīng)用最速下降法反向求取炮檢點(diǎn)間的射線路徑[11]，在網(wǎng)格單元內(nèi)認(rèn)為射線路徑為直線段，這樣整體而言，炮檢點(diǎn)之間的射線路徑則是由一系列的直線段首尾依次連接形成的彎曲射線(圖3)，不妨設(shè)二維平面被剖分成N個(gè)網(wǎng)格，每個(gè)網(wǎng)格的慢度為sj，獲得了M條射線，假設(shè)第i條射線在第j個(gè)網(wǎng)格的射線長(zhǎng)度為ai,j，該射線的初至波旅行時(shí)為ti，則可以得到該射線的旅行時(shí)方程如式(5)所示:
圖3 離散網(wǎng)格中彎曲射線層析反演示意Fig.3 Sketch of curve-ray tomography inversion in regular grids
(5)
式中，ai,j為第i條射線在第j個(gè)網(wǎng)格中的射線長(zhǎng)度，Sj為第j個(gè)網(wǎng)格的慢度，ti為第i條射線的初至?xí)r，
將所有射線方程聯(lián)立可以得到線性方程組(6)：
A·S=T,
(6)
其中,A=(ai,j)M,N為射線路徑矩陣，S=(sj)N為慢度，列向量T=(ti)M為初至?xí)r間列向量，式中A通常為一大型稀疏病態(tài)矩陣，可采用ART、SVD、SIRT、LSQR、LSCG等算法對(duì)其進(jìn)行求解，本文采用LSCG算法對(duì)慢度場(chǎng)進(jìn)行更新，該算法是共軛梯度法的一種改進(jìn)算法，對(duì)于大型線性方程組Ax=b，可將其轉(zhuǎn)化為二次函數(shù)的極小點(diǎn)問題如式(7)所示[14]：
(7)
函數(shù)f(x)的梯度g(x)為：
g(x)=f(x)=ATAx-ATb。
(8)
對(duì)于任意非零向量p和實(shí)數(shù)t，有：
于是有：
(10)
考慮到p的隨機(jī)性，必須有g(shù)(u)=0成立，因而u為方程組的解，對(duì)于共軛向量序列{p(k)}，假定x(0)是任一給定的初始向量，對(duì)于k=0，1，2，3，…，以x(k)作為起點(diǎn)，沿方向向量p(k)求函數(shù)f(x)在直線x=x(k)+tp(k)上的極小點(diǎn)，可以得到：
x(k+1)=x(k)+αkp(k)，
(11)
s(k)=AT(b-Ax(k)) ，
(12)
(13)
以上3式中p(k)表示搜索方向，如果取p(0)=s(0)，則有：
p(k+1)=s(k+1)+βkp(k)，
(14)
(15)
式(11)～(15)構(gòu)成了LSCG算法的迭代格式，在迭代之前常給定初始解x(0)=0。從該迭代格式可知LSCG算法主要由大量Ax、ATx及矢量?jī)?nèi)積運(yùn)算構(gòu)成，為了提高計(jì)算速度降低內(nèi)存消耗，在程序設(shè)計(jì)中采用三元組存儲(chǔ)方式存儲(chǔ)距離陣A[15]。為了壓制迭代過程中出現(xiàn)的反演噪聲及野值，提高反演剖面的信噪比，采用均值濾波器對(duì)迭代的中間結(jié)果進(jìn)行了平滑處理[16]。
1.3 克?；舴蚍e分偏移
應(yīng)用有限差分法分別計(jì)算以炮點(diǎn)和檢波點(diǎn)作為激發(fā)點(diǎn)時(shí)各成像網(wǎng)格點(diǎn)的旅行時(shí)場(chǎng)后，可基于Kirchhof積分公式通過邊界積分法實(shí)現(xiàn)地震數(shù)據(jù)成像，對(duì)于二維全平面域偏移成像問題，其離散形式如式(16)所示[17]：
(16)
值得注意的是：由于在偏移過程中對(duì)數(shù)據(jù)已做擴(kuò)散補(bǔ)償，因此在預(yù)處理時(shí)不應(yīng)對(duì)數(shù)據(jù)再做擴(kuò)散補(bǔ)償處理。
2 數(shù)值模擬
2.1 模型建立
建立二維速度模型如圖4所示，該模型橫向長(zhǎng)1 100 m，深700 m，背景速度2 200 m·s-1，在模型上設(shè)置5個(gè)速度異常體，其中4個(gè)正方形速度異常體邊長(zhǎng)均為100 m，速度為3 000 m·s-1，模型中間長(zhǎng)條狀異常體長(zhǎng)500 m、寬10 m、速度1 700 m·s-1。在模型左側(cè)及頂部每隔20 m激發(fā)一單炮，地震子波為主頻60 Hz的Ricker子波，共激發(fā)了92張模擬記錄，并于模型頂部及右側(cè)每隔10 m安置檢波器接收地震信號(hào)，記錄長(zhǎng)度1 500 ms，采樣間隔0.5 ms；圖5為激發(fā)點(diǎn)位于(0 m，120 m)處的模擬單炮記錄。由該單炮記錄可以看出，對(duì)于圖4的簡(jiǎn)單模型，地震波場(chǎng)特征比較復(fù)雜，單炮記錄上發(fā)育有直達(dá)波、反射波、繞射波等。
圖4 模型及觀測(cè)系統(tǒng)Fig.4 Model and observation system
圖5 模擬單炮記錄Fig.5 Simulate single shot record
2.2 層析反演速度結(jié)構(gòu)
分別采用直射線和彎曲射線層析成像方法反演得到了成像區(qū)域的速度場(chǎng)，反演網(wǎng)格大小5 m×5 m，采用最小二乘共軛梯度(LSCG)反演算法共進(jìn)行了10次迭代處理。采用二維平滑濾波器對(duì)迭代的中間結(jié)果進(jìn)行平滑處理，以消除反演過程中的噪聲及野值，提高反演剖面的信噪比。圖6是兩種層析方法獲得的反演速度場(chǎng)(兩張剖面速度色標(biāo)一致)，圖中黑線為模型中速度異常體的邊界，從該圖可看到：兩種層析成像方法獲得圖像的背景速度大約為2 200m·s-1，與模型的背景速度一致，兩種方法4個(gè)高速異常體在成像速度平面上均有顯示；直射線層析剖面上高速異常體速度大約為2 600 m·s-1，較模型速度低了約400 m·s-1，異常體的邊界比較模糊，且較真實(shí)尺寸偏大，形態(tài)畸變嚴(yán)重，模型中間的長(zhǎng)條形低速異常體在反演剖面上并無明顯異常；在彎曲射線層析成像剖面上，模型頂部?jī)筛咚佼惓ｓw的顯示較直射線層析方法的分辨率高，邊界形態(tài)刻畫清晰，異常體速度大致為2 800 m·s-1，精度較直射線方法有所提高，異常體大小與其實(shí)際尺寸較接近，模型底部?jī)伤俣犬惓ｓw由于觀測(cè)方位的限制，畸變較嚴(yán)重，但仍較直射線成像方法更加收斂，模型中間的條形低速異常體成像也有所改善，異常體速度大約為1 900 m·s-1，較真實(shí)速度略高。
圖6 直射線(a)與彎曲射線(b)層析反演速度結(jié)構(gòu)Fig.6 Velocity structure by straight(a) and curve ray(b) tomography inversion
圖7為炮點(diǎn)位于(500 m，0 m)處直射線層析和彎曲射線層析速度模型初至波時(shí)間等值線圖，由該圖可見彎曲射線層成結(jié)果的等時(shí)線與理論等時(shí)線的吻合程度明顯優(yōu)于直射線法，由于模型底部沒有炮點(diǎn)和接收點(diǎn)，觀測(cè)方位受限，兩種反演方法獲得的速度均不夠準(zhǔn)確，初至?xí)r間的誤差較大。
圖7 直射線(a)與彎曲射線(b)層析等時(shí)線圖對(duì)比Fig.7 Comparison of isochron maps of straight (a) and curve (b) ray tomography
2.3 克?；舴蚱瞥上裥Ч麑?duì)比
分別以基于直射線和彎曲射線反演的速度作為偏移速度，對(duì)模擬單炮數(shù)據(jù)進(jìn)行時(shí)間偏移成像，成像網(wǎng)格為5 m×5 m。圖8是兩種不同偏移速度獲得的原始成像剖面，可以看到：兩種方法的原始偏移剖面上均存在較強(qiáng)的低頻背景干擾，剖面信噪比和分辨率較低，筆者借鑒逆時(shí)偏移的低頻噪聲壓制技術(shù)，應(yīng)用Laplace差分濾波器對(duì)原始成像剖面進(jìn)行了濾波處理，結(jié)果如圖9所示，濾波后圖像的信噪比和分辨率得到一定增強(qiáng)，速度異常體的邊界更加清楚。對(duì)比圖9中的兩張剖面可以看出：兩個(gè)偏移結(jié)果對(duì)所有速度異常體均有反映，基于直射線反演速度的成像剖面上4個(gè)高速異常體的形態(tài)發(fā)生了嚴(yán)重畸變，模型中間的長(zhǎng)條形低速異常體反射波聚焦較差；以彎曲線CT反演的速度為偏移速度的成像剖面上，4個(gè)高速異常體的形狀和大小均與原模型比較吻合，模型中間的長(zhǎng)條型低速異常體上半部分成像清晰，下半部由于觀測(cè)系統(tǒng)的原因，未能接收到有效反射波而不能成像。
a—直射線層析;b—彎曲射線層析a—imaging profile of straight ray result;b—imaging profile of curve ray result
a—直射線層析;b—彎曲射線層析a—imaging profile of straight ray result;b—imaging profile of curve ray result
3 結(jié)論
本文以層析反演速度作為偏移速度，應(yīng)用克?；舴蚍e分偏移方法對(duì)復(fù)雜二維模型進(jìn)行了全平面偏移成像，通過對(duì)層析反演結(jié)果及偏移成像剖面的對(duì)比分析得到如下認(rèn)識(shí)：
1) 對(duì)于復(fù)雜速度結(jié)構(gòu)，直射線層析反演方法獲得的速度異常的形態(tài)大小會(huì)發(fā)生嚴(yán)重畸變，且對(duì)小構(gòu)造刻畫不清，而彎曲射線層析算法可以較準(zhǔn)確地反映速度異常體的形態(tài)和位置，對(duì)小構(gòu)造的刻畫較直射線清晰。
2) 克?；舴蚍e分法偏移技術(shù)是一種觀測(cè)系統(tǒng)適應(yīng)性強(qiáng)的偏移成像方法，以彎線射線層析反演速度作為偏移速度，采用克?；舴虔B前偏移方法可以獲得較準(zhǔn)確的偏移成像剖面。
3) 應(yīng)用Laplace差分濾波器可以對(duì)原始偏移剖面上的低頻噪聲進(jìn)行有效壓制，提高成像剖面的信噪比和清晰度。
猜你喜歡
反演對(duì)稱變換在解決平面幾何問題中的應(yīng)用中等數(shù)學(xué)(2022年5期)2022-08-29 深層膏鹽體局部層析速度建模物探與化探(2022年4期)2022-08-26 地震折射層析法在紅壤關(guān)鍵帶地層劃分中的應(yīng)用研究*土壤學(xué)報(bào)(2022年1期)2022-03-08 一體化綠葉層析裝置發(fā)明與創(chuàng)新(2022年1期)2022-02-20 基于ADS-B的風(fēng)場(chǎng)反演與異常值影響研究成都信息工程大學(xué)學(xué)報(bào)(2021年5期)2021-12-30 Meteo-particle模型在ADS-B風(fēng)場(chǎng)反演中的性能研究成都信息工程大學(xué)學(xué)報(bào)(2021年5期)2021-12-30 利用錐模型反演CME三維參數(shù)空間科學(xué)學(xué)報(bào)(2021年4期)2021-08-30 多維空間及多維射線坐標(biāo)系設(shè)想學(xué)校教育研究(2021年21期)2021-05-24 包涵體蛋白的純化層析復(fù)性技術(shù)研究進(jìn)展商品與質(zhì)量(2018年38期)2018-04-16 話說線段、射線、直線初中生世界·七年級(jí)(2016年2期)2016-03-03
智敏,朱建剛
(中煤科工集團(tuán)西安研究院有限公司,陜西 西安 710077)
0 引言
克希霍夫積分偏移是一種基于波動(dòng)方程積分解的偏移成像方法，與其他偏移方法相比，該方法能夠適應(yīng)非規(guī)則觀測(cè)系統(tǒng)，具有較高的計(jì)算效率[1]。目前主流處理系統(tǒng)的克?；舴蚱颇K主要針對(duì)于地面地震數(shù)據(jù)成像，要求炮檢點(diǎn)均位于地表，波場(chǎng)信息來自地下半空間，屬于半空間成像方法。對(duì)于某些特殊地震觀測(cè)系統(tǒng)，如煤礦井下槽波觀測(cè)系統(tǒng)，地震信號(hào)可能來自以炮檢點(diǎn)為焦點(diǎn)的完整橢圓上，屬于全平面域偏移成像，主流處理軟件并不適用。王季、姬廣忠等基于直射線層析重構(gòu)的槽波速度場(chǎng)，采用直射線克?；舴蚍e分偏移方法對(duì)井下反射槽波數(shù)據(jù)進(jìn)行了偏移成像，取得了一定的效果[2-3]，該方法與地面地震直射線偏移成像一樣，將繞射點(diǎn)與炮檢點(diǎn)之間的射線路徑視為直線，沒有考慮速度場(chǎng)彎化引起的射線彎曲及旅行時(shí)變化，造成成像點(diǎn)與原繞射點(diǎn)不重合，降低了成像分辨率。
針對(duì)直射線層析成像速度的這一缺陷，本文采用彎曲射線層析成像方法反演成像區(qū)域速度。對(duì)于彎曲射線層析成像，射線追蹤的精度直接影響著成像精度。常用的射線追蹤方法有試射法、旅行時(shí)線性插值法(LTI)、快速步進(jìn)法(FMM)、有限差分求解程函方程法。試射法對(duì)于簡(jiǎn)單速度模型具有較高精度，但不能處理折射波問題，對(duì)于復(fù)雜構(gòu)造常存在追蹤盲區(qū)，并且計(jì)算效率低下[4]。LTI、FMM算法及有限差分法均首先計(jì)算從震源傳至各節(jié)點(diǎn)的最小旅行時(shí)，并基于旅行時(shí)間采用反向追蹤技術(shù)求取各檢波點(diǎn)到炮點(diǎn)的射線路徑。在計(jì)算最小旅行時(shí)時(shí)，LTI算法常采用按列(行)掃描的方式計(jì)算全局最小旅行時(shí)[5]，該方法未能考慮來自各個(gè)方向的射線，對(duì)于橫向速度變化較大的介質(zhì)，計(jì)算結(jié)果往往并非全局最小，葉佩等采用了全方位循環(huán)法求得了網(wǎng)格節(jié)點(diǎn)的全局最小旅行時(shí)[6]，張建中等采用Dijkstra算法逐點(diǎn)計(jì)算了各節(jié)點(diǎn)的最小旅行時(shí)[7]；FMM算法通過上風(fēng)差分近似和窄帶技術(shù)近似模擬地震波前的傳播，具有計(jì)算精度高，運(yùn)算速度快、無條件穩(wěn)定的優(yōu)點(diǎn)，可以計(jì)算網(wǎng)格節(jié)的全局最小旅行時(shí)，但該方法近源區(qū)域網(wǎng)格節(jié)點(diǎn)旅行計(jì)算誤差較大，并對(duì)后續(xù)節(jié)點(diǎn)的計(jì)算精度有影響，需要對(duì)近源區(qū)的網(wǎng)格進(jìn)行行加密處理，以提高近源區(qū)旅行時(shí)的計(jì)算精度[8]；有限差分求解程函方程法是一種計(jì)算精度高、運(yùn)算量小的旅行時(shí)計(jì)算方法[9]，對(duì)于速度變化大的介質(zhì)，采用傳統(tǒng)的波前矩形擴(kuò)展方可能因違返波傳播的因果關(guān)系導(dǎo)致計(jì)算不穩(wěn)定，潘紀(jì)順等通過在算法中額外考慮首波、體波和散射波對(duì)旅行時(shí)影響的方式有效地解決了這一問題[10]。本文采用有限差分法求解程函方程的方法求取成像平面內(nèi)各節(jié)點(diǎn)的初至波旅行時(shí)場(chǎng)，算法上借鑒了潘紀(jì)順的處理方法，波前擴(kuò)展方法采用Disjkstra算法逐點(diǎn)計(jì)算各節(jié)點(diǎn)的全局最小旅行，射線反向追蹤采用了Aldridge提出的最速下降法求取射線路徑的方法[11]，并基于全局最小旅行時(shí)對(duì)二維復(fù)雜速度模型進(jìn)行全平面域克?；舴蚩臻g域偏移成像，通過對(duì)偏移結(jié)果的對(duì)比分析，認(rèn)為該方法可以較真實(shí)地反映異常速度體的形態(tài)大小和位置，成像結(jié)果分辨率較高。
1 方法原理
與常規(guī)地面地震克?；舴蚍e分法疊前深度偏移類似，全平面域克?；舴蚍e分疊前偏移成像方法由偏移速度求取和克?；舴蚍e分偏移兩部分組成，兩者在實(shí)現(xiàn)過程中均涉及成像區(qū)域旅行時(shí)場(chǎng)計(jì)算，本文采用有限差分求解程函方程的方法計(jì)算節(jié)點(diǎn)旅行時(shí)。
1.1 有限差分程函方程初至旅行時(shí)計(jì)算
旅行時(shí)計(jì)算是積分法偏移成像的重要組成部分，其計(jì)算精度直接影響著成像效果，在高頻近似下，平面波波前的傳播符合程函方程式(1)[9]：
(1)
式中,t為初至?xí)r間；s(x,y)為平面模型的慢度場(chǎng)。
采用有限差分近似方法求解式(1)，將速度模型離散成規(guī)則的矩形網(wǎng)格(圖1所示)，X、Y方向網(wǎng)格邊長(zhǎng)分別為Δx、Δy，將每個(gè)網(wǎng)格單元內(nèi)的慢度視為常數(shù)，將方程中的兩個(gè)時(shí)間微分算子用有限差分算子近似，如式(2)～(3)所示[9,12]:
圖1 矩形網(wǎng)格剖分示意Fig.1 Diagram of rectangular mesh partition
(2)
(3)
將式(2)、(3)代入式(1)，經(jīng)推導(dǎo)可以得到式(4)[9]:
(4)
式中,A=Δx2-Δz2，B=2ΔxΔz，C=Δx2+Δz2，已知矩形網(wǎng)格3個(gè)節(jié)點(diǎn)旅行時(shí)t0、t1、t2可以利用式(4)得到旅行時(shí)t3，計(jì)算中常常會(huì)出現(xiàn)根號(hào)下為負(fù)值而出現(xiàn)計(jì)算不穩(wěn)定的情況，此時(shí)可將節(jié)點(diǎn)0、1、2視為散射源，按直射傳播方式，分別計(jì)算節(jié)點(diǎn)3的初至?xí)r間，結(jié)合節(jié)點(diǎn)3的已有旅行時(shí)間，取5個(gè)值里面的最小值作為節(jié)點(diǎn)3的全局最小旅行時(shí)。
在旅行時(shí)場(chǎng)外推過程中，筆者采用Dijkstra算法逐點(diǎn)計(jì)算所有節(jié)點(diǎn)的全局最小旅行時(shí)，將所有計(jì)算節(jié)點(diǎn)分成3個(gè)子集，分別用P、Q、R表示，其中P為尚未計(jì)算旅行時(shí)的節(jié)點(diǎn)集合(節(jié)點(diǎn)旅行為無窮大)，Q為已計(jì)算旅行時(shí)的節(jié)點(diǎn)集合(不確定是否為全局最小旅行時(shí))，R為已確定為全局最小旅行時(shí)節(jié)點(diǎn)的集合，采用如下方法對(duì)各節(jié)點(diǎn)的集合屬性進(jìn)行更新：
1)將所有節(jié)點(diǎn)置于P集合，且所有節(jié)點(diǎn)旅行時(shí)均置為無窮大，計(jì)算炮點(diǎn)所在網(wǎng)格4個(gè)節(jié)點(diǎn)的旅行時(shí)，將其置于Q集；
2)在Q集中選取旅行時(shí)最小的節(jié)點(diǎn)作為次級(jí)震源，并將其置于R集，計(jì)算該節(jié)點(diǎn)周圍8個(gè)節(jié)點(diǎn)的旅行時(shí)，若待計(jì)算節(jié)點(diǎn)屬于P集合，則將其置于Q集合，更新該節(jié)點(diǎn)旅行時(shí)；若節(jié)點(diǎn)屬于Q集合，則將計(jì)算結(jié)果與先前結(jié)果比較，取小值作為該節(jié)點(diǎn)旅行時(shí)；若節(jié)點(diǎn)屬于R集，則不改變?cè)担?/p>
3)重復(fù)步驟2)直至所有節(jié)點(diǎn)均屬于R集為止，即可求得成像區(qū)域內(nèi)所有節(jié)點(diǎn)的全局最小旅行時(shí)。
1.2 初至波層析成像反演速度場(chǎng)
偏移速度場(chǎng)的計(jì)算精度直接影響著偏移成像效果，是克?；舴蚍e分法偏移成像的關(guān)鍵，對(duì)于非地面地震觀測(cè)系統(tǒng)，采用地面地震觀測(cè)系統(tǒng)的偏移速度分析方法獲取偏移速度比較困難，通常采用初至波層析反演方法獲得偏移速度，其中直射線層析成像方法應(yīng)用比較普遍，如文獻(xiàn)[2-3]中的井下槽波繞射偏移成像，當(dāng)反演區(qū)域速度差異較小時(shí)，可近似認(rèn)為射線在介質(zhì)中沿直射傳播，如圖2所示，炮檢點(diǎn)間的射線路徑為直線，ai,j為第i條射線在第j個(gè)網(wǎng)格中的射線長(zhǎng)度，Sj為第j個(gè)網(wǎng)格的慢度，ti為第i條射線的初至?xí)r，利用直射線層析反演方法可以獲得較可靠的成像速度，當(dāng)成像區(qū)域速度變化較大時(shí)，由于實(shí)際射線路徑常繞過低速區(qū)，并向高速區(qū)集中，射線路徑發(fā)生較大彎曲，不滿足直射線假設(shè)條件，常導(dǎo)致高速異常體偏大，低速異常體偏小[13]。
圖2 離散網(wǎng)格中直射線層析反演示意Fig.2 Sketch of straight-ray tomography inversion in regular grids
為了克服直射線層析的這一缺陷，將成像區(qū)域進(jìn)行網(wǎng)格剖分后，采用有限差分法求得各成像網(wǎng)格節(jié)點(diǎn)的全局最小旅行時(shí)，應(yīng)用最速下降法反向求取炮檢點(diǎn)間的射線路徑[11]，在網(wǎng)格單元內(nèi)認(rèn)為射線路徑為直線段，這樣整體而言，炮檢點(diǎn)之間的射線路徑則是由一系列的直線段首尾依次連接形成的彎曲射線(圖3)，不妨設(shè)二維平面被剖分成N個(gè)網(wǎng)格，每個(gè)網(wǎng)格的慢度為sj，獲得了M條射線，假設(shè)第i條射線在第j個(gè)網(wǎng)格的射線長(zhǎng)度為ai,j，該射線的初至波旅行時(shí)為ti，則可以得到該射線的旅行時(shí)方程如式(5)所示:
圖3 離散網(wǎng)格中彎曲射線層析反演示意Fig.3 Sketch of curve-ray tomography inversion in regular grids
(5)
式中，ai,j為第i條射線在第j個(gè)網(wǎng)格中的射線長(zhǎng)度，Sj為第j個(gè)網(wǎng)格的慢度，ti為第i條射線的初至?xí)r，
將所有射線方程聯(lián)立可以得到線性方程組(6)：
A·S=T,
(6)
其中,A=(ai,j)M,N為射線路徑矩陣，S=(sj)N為慢度，列向量T=(ti)M為初至?xí)r間列向量，式中A通常為一大型稀疏病態(tài)矩陣，可采用ART、SVD、SIRT、LSQR、LSCG等算法對(duì)其進(jìn)行求解，本文采用LSCG算法對(duì)慢度場(chǎng)進(jìn)行更新，該算法是共軛梯度法的一種改進(jìn)算法，對(duì)于大型線性方程組Ax=b，可將其轉(zhuǎn)化為二次函數(shù)的極小點(diǎn)問題如式(7)所示[14]：
(7)
函數(shù)f(x)的梯度g(x)為：
g(x)=f(x)=ATAx-ATb。
(8)
對(duì)于任意非零向量p和實(shí)數(shù)t，有：
于是有：
(10)
考慮到p的隨機(jī)性，必須有g(shù)(u)=0成立，因而u為方程組的解，對(duì)于共軛向量序列{p(k)}，假定x(0)是任一給定的初始向量，對(duì)于k=0，1，2，3，…，以x(k)作為起點(diǎn)，沿方向向量p(k)求函數(shù)f(x)在直線x=x(k)+tp(k)上的極小點(diǎn)，可以得到：
x(k+1)=x(k)+αkp(k)，
(11)
s(k)=AT(b-Ax(k)) ，
(12)
(13)
以上3式中p(k)表示搜索方向，如果取p(0)=s(0)，則有：
p(k+1)=s(k+1)+βkp(k)，
(14)
(15)
式(11)～(15)構(gòu)成了LSCG算法的迭代格式，在迭代之前常給定初始解x(0)=0。從該迭代格式可知LSCG算法主要由大量Ax、ATx及矢量?jī)?nèi)積運(yùn)算構(gòu)成，為了提高計(jì)算速度降低內(nèi)存消耗，在程序設(shè)計(jì)中采用三元組存儲(chǔ)方式存儲(chǔ)距離陣A[15]。為了壓制迭代過程中出現(xiàn)的反演噪聲及野值，提高反演剖面的信噪比，采用均值濾波器對(duì)迭代的中間結(jié)果進(jìn)行了平滑處理[16]。
1.3 克?；舴蚍e分偏移
應(yīng)用有限差分法分別計(jì)算以炮點(diǎn)和檢波點(diǎn)作為激發(fā)點(diǎn)時(shí)各成像網(wǎng)格點(diǎn)的旅行時(shí)場(chǎng)后，可基于Kirchhof積分公式通過邊界積分法實(shí)現(xiàn)地震數(shù)據(jù)成像，對(duì)于二維全平面域偏移成像問題，其離散形式如式(16)所示[17]：
(16)
值得注意的是：由于在偏移過程中對(duì)數(shù)據(jù)已做擴(kuò)散補(bǔ)償，因此在預(yù)處理時(shí)不應(yīng)對(duì)數(shù)據(jù)再做擴(kuò)散補(bǔ)償處理。
2 數(shù)值模擬
2.1 模型建立
建立二維速度模型如圖4所示，該模型橫向長(zhǎng)1 100 m，深700 m，背景速度2 200 m·s-1，在模型上設(shè)置5個(gè)速度異常體，其中4個(gè)正方形速度異常體邊長(zhǎng)均為100 m，速度為3 000 m·s-1，模型中間長(zhǎng)條狀異常體長(zhǎng)500 m、寬10 m、速度1 700 m·s-1。在模型左側(cè)及頂部每隔20 m激發(fā)一單炮，地震子波為主頻60 Hz的Ricker子波，共激發(fā)了92張模擬記錄，并于模型頂部及右側(cè)每隔10 m安置檢波器接收地震信號(hào)，記錄長(zhǎng)度1 500 ms，采樣間隔0.5 ms；圖5為激發(fā)點(diǎn)位于(0 m，120 m)處的模擬單炮記錄。由該單炮記錄可以看出，對(duì)于圖4的簡(jiǎn)單模型，地震波場(chǎng)特征比較復(fù)雜，單炮記錄上發(fā)育有直達(dá)波、反射波、繞射波等。
圖4 模型及觀測(cè)系統(tǒng)Fig.4 Model and observation system
圖5 模擬單炮記錄Fig.5 Simulate single shot record
2.2 層析反演速度結(jié)構(gòu)
分別采用直射線和彎曲射線層析成像方法反演得到了成像區(qū)域的速度場(chǎng)，反演網(wǎng)格大小5 m×5 m，采用最小二乘共軛梯度(LSCG)反演算法共進(jìn)行了10次迭代處理。采用二維平滑濾波器對(duì)迭代的中間結(jié)果進(jìn)行平滑處理，以消除反演過程中的噪聲及野值，提高反演剖面的信噪比。圖6是兩種層析方法獲得的反演速度場(chǎng)(兩張剖面速度色標(biāo)一致)，圖中黑線為模型中速度異常體的邊界，從該圖可看到：兩種層析成像方法獲得圖像的背景速度大約為2 200m·s-1，與模型的背景速度一致，兩種方法4個(gè)高速異常體在成像速度平面上均有顯示；直射線層析剖面上高速異常體速度大約為2 600 m·s-1，較模型速度低了約400 m·s-1，異常體的邊界比較模糊，且較真實(shí)尺寸偏大，形態(tài)畸變嚴(yán)重，模型中間的長(zhǎng)條形低速異常體在反演剖面上并無明顯異常；在彎曲射線層析成像剖面上，模型頂部?jī)筛咚佼惓ｓw的顯示較直射線層析方法的分辨率高，邊界形態(tài)刻畫清晰，異常體速度大致為2 800 m·s-1，精度較直射線方法有所提高，異常體大小與其實(shí)際尺寸較接近，模型底部?jī)伤俣犬惓ｓw由于觀測(cè)方位的限制，畸變較嚴(yán)重，但仍較直射線成像方法更加收斂，模型中間的條形低速異常體成像也有所改善，異常體速度大約為1 900 m·s-1，較真實(shí)速度略高。
圖6 直射線(a)與彎曲射線(b)層析反演速度結(jié)構(gòu)Fig.6 Velocity structure by straight(a) and curve ray(b) tomography inversion
圖7為炮點(diǎn)位于(500 m，0 m)處直射線層析和彎曲射線層析速度模型初至波時(shí)間等值線圖，由該圖可見彎曲射線層成結(jié)果的等時(shí)線與理論等時(shí)線的吻合程度明顯優(yōu)于直射線法，由于模型底部沒有炮點(diǎn)和接收點(diǎn)，觀測(cè)方位受限，兩種反演方法獲得的速度均不夠準(zhǔn)確，初至?xí)r間的誤差較大。
圖7 直射線(a)與彎曲射線(b)層析等時(shí)線圖對(duì)比Fig.7 Comparison of isochron maps of straight (a) and curve (b) ray tomography
2.3 克?；舴蚱瞥上裥Ч麑?duì)比
分別以基于直射線和彎曲射線反演的速度作為偏移速度，對(duì)模擬單炮數(shù)據(jù)進(jìn)行時(shí)間偏移成像，成像網(wǎng)格為5 m×5 m。圖8是兩種不同偏移速度獲得的原始成像剖面，可以看到：兩種方法的原始偏移剖面上均存在較強(qiáng)的低頻背景干擾，剖面信噪比和分辨率較低，筆者借鑒逆時(shí)偏移的低頻噪聲壓制技術(shù)，應(yīng)用Laplace差分濾波器對(duì)原始成像剖面進(jìn)行了濾波處理，結(jié)果如圖9所示，濾波后圖像的信噪比和分辨率得到一定增強(qiáng)，速度異常體的邊界更加清楚。對(duì)比圖9中的兩張剖面可以看出：兩個(gè)偏移結(jié)果對(duì)所有速度異常體均有反映，基于直射線反演速度的成像剖面上4個(gè)高速異常體的形態(tài)發(fā)生了嚴(yán)重畸變，模型中間的長(zhǎng)條形低速異常體反射波聚焦較差；以彎曲線CT反演的速度為偏移速度的成像剖面上，4個(gè)高速異常體的形狀和大小均與原模型比較吻合，模型中間的長(zhǎng)條型低速異常體上半部分成像清晰，下半部由于觀測(cè)系統(tǒng)的原因，未能接收到有效反射波而不能成像。
a—直射線層析;b—彎曲射線層析a—imaging profile of straight ray result;b—imaging profile of curve ray result
a—直射線層析;b—彎曲射線層析a—imaging profile of straight ray result;b—imaging profile of curve ray result
3 結(jié)論
本文以層析反演速度作為偏移速度，應(yīng)用克?；舴蚍e分偏移方法對(duì)復(fù)雜二維模型進(jìn)行了全平面偏移成像，通過對(duì)層析反演結(jié)果及偏移成像剖面的對(duì)比分析得到如下認(rèn)識(shí)：
1) 對(duì)于復(fù)雜速度結(jié)構(gòu)，直射線層析反演方法獲得的速度異常的形態(tài)大小會(huì)發(fā)生嚴(yán)重畸變，且對(duì)小構(gòu)造刻畫不清，而彎曲射線層析算法可以較準(zhǔn)確地反映速度異常體的形態(tài)和位置，對(duì)小構(gòu)造的刻畫較直射線清晰。
2) 克?；舴蚍e分法偏移技術(shù)是一種觀測(cè)系統(tǒng)適應(yīng)性強(qiáng)的偏移成像方法，以彎線射線層析反演速度作為偏移速度，采用克?；舴虔B前偏移方法可以獲得較準(zhǔn)確的偏移成像剖面。
3) 應(yīng)用Laplace差分濾波器可以對(duì)原始偏移剖面上的低頻噪聲進(jìn)行有效壓制，提高成像剖面的信噪比和清晰度。
猜你喜歡
中等數(shù)學(xué)(2022年5期)2022-08-29
物探與化探(2022年4期)2022-08-26
土壤學(xué)報(bào)(2022年1期)2022-03-08
發(fā)明與創(chuàng)新(2022年1期)2022-02-20
成都信息工程大學(xué)學(xué)報(bào)(2021年5期)2021-12-30
成都信息工程大學(xué)學(xué)報(bào)(2021年5期)2021-12-30
空間科學(xué)學(xué)報(bào)(2021年4期)2021-08-30
學(xué)校教育研究(2021年21期)2021-05-24
商品與質(zhì)量(2018年38期)2018-04-16
初中生世界·七年級(jí)(2016年2期)2016-03-03
