基于坐標(biāo)變換的三維旋轉(zhuǎn)曲率屬性計算方法

張曉琦, 文曉濤, 張超銘, 何易龍, 王錦濤

(成都理工大學(xué) a.地球物理學(xué)院,b.油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都 610059)

0 引言

裂縫是重要的油氣儲集空間與運(yùn)移通道，因此有效地識別裂縫已成為裂縫型油氣藏研究的重點與難點。一般而言，裂縫與構(gòu)造有密切的關(guān)系，在構(gòu)造高部位及斷層附近易發(fā)育裂縫，而在這些部位，同相軸的曲率與構(gòu)造其他部位也有差異，因此曲率屬性分析是一種較有效的裂縫識別方法[1-2]。Robert[3]總結(jié)了各種曲率屬性的定義和計算方法，分析了曲率在三維解釋中的應(yīng)用，發(fā)展了一代傳統(tǒng)的二維曲面曲率計算方法；Al-Dossary等[4]發(fā)展了第二代的體曲率計算方法，該方法直接利用地震振幅計算體曲率，驗證了曲率屬性對地下裂縫以及裂縫型碳酸鹽巖儲層有良好的反映能力；Marfurt等[5]將相干屬性和曲率屬性進(jìn)行對比融合顯示，能較好地突出地質(zhì)構(gòu)造特征，取得了很好的效果；陳學(xué)華[6]改進(jìn)了算法，將多尺度自適應(yīng)微分算子引入體曲率計算中，充分利用不同尺度的曲率來突出地質(zhì)異常，因此此方法可以用來描述地震在時間和空間上的多尺度特征；劉松鳴[7]利用推導(dǎo)的方位濾波算子對體曲率屬性進(jìn)行方位加強(qiáng)處理，其主要特點是突出了特定方向的地質(zhì)體異常信息，同時也抑制了地震背景干擾。

一般來說，常規(guī)曲率算法在描述平緩地層的幾何形態(tài)上還是比較可靠的，但如果地層陡峭傾斜甚至于反轉(zhuǎn)，就會有很多算法高估曲率屬性，而這常常是因為傳統(tǒng)算法對相關(guān)的方位角及傾角做出了不準(zhǔn)確的估計，進(jìn)而影響了區(qū)域曲率精度。研究者們對此也提出了諸多辦法來解決這一方面的問題，如進(jìn)行大傾角掃描[8]、離心窗掃描[9]，或提出為消除傾角影響采用振幅曲率加權(quán)算法[10]等各種優(yōu)化改進(jìn)方法。雖然上述方法在裂縫識別中可以取得一定的效果，但總體來看，常規(guī)的曲率計算在處理地震數(shù)據(jù)時均在地球直角坐標(biāo)系中進(jìn)行處理，計算方法復(fù)雜，且方位角傾角的計算與使用存在較大的誤差，不利于裂縫斷層的準(zhǔn)確描述。為了彌補(bǔ)這些方法的不足，增強(qiáng)曲率屬性的提取精度，筆者提出了基于坐標(biāo)變換的三維旋轉(zhuǎn)坐標(biāo)變換曲率計算方法，將構(gòu)造面進(jìn)行三維旋轉(zhuǎn)坐標(biāo)變換，借助構(gòu)造傾角將局部構(gòu)造曲面轉(zhuǎn)換為平面，間接減少角度影響。這樣既簡化了計算，減少了誤差，又提高了曲率屬性的分辨率，即增強(qiáng)了曲率屬性的提取精度。

1 基本原理

體曲率屬性[11]作為直接描述斷層和預(yù)測裂縫走向及分布的一個幾何屬性，在地震勘探中占有重要地位。三維地震數(shù)據(jù)的采集和顯示通常是在直角坐標(biāo)系中，如圖1所示，z軸為深度方向，θ角為真傾角，φ角為從x軸出發(fā)的方位角。

圖1 傾角及方位角關(guān)系圖

根據(jù)Marfurt等理論，三維空間中的特定某點處的曲率值[12-13]，是利用這一點及其相鄰的采樣點的視傾角的擬合方程進(jìn)行計算。然后通過使用最小二乘法近似的獲得二次曲面方程：

z(x,y)=ax2+by2+cxy+dx+ey+f

(1)

通過計算上述微分方程，趨勢面系數(shù)計算如下：

(2)

所以傾角θ和方位角φ分別表示為：

(3)

(4)

這里是在基于地球直角坐標(biāo)系旋轉(zhuǎn)[14]后的新坐標(biāo)系中對每個樣本的二階導(dǎo)數(shù)，采用復(fù)地震道分析方法進(jìn)行曲率屬性計算，在旋轉(zhuǎn)后的新系統(tǒng)中，其中X1軸沿傾斜方位角方向，Z1軸與過原點的表面垂直，x0、y0、z0表示原始坐標(biāo)系中坐標(biāo)參數(shù)，x1、y1、z1表示坐標(biāo)系更新后的坐標(biāo)參數(shù)。具體操作步驟為：將Z0軸作為旋轉(zhuǎn)軸，方位角φ為旋轉(zhuǎn)角度，對原始的坐標(biāo)系統(tǒng)X0-Y0-Z0沿逆時針(右側(cè))進(jìn)行旋轉(zhuǎn)，得到新的X1-Y1-Z0系統(tǒng)。 然后將Y1作為旋轉(zhuǎn)軸，傾角θ為旋轉(zhuǎn)角度對X1-Y1-Z0系統(tǒng)沿逆時針(右側(cè))進(jìn)行旋轉(zhuǎn)，經(jīng)過兩次旋轉(zhuǎn)得到新的X1-Y1-Z1系統(tǒng)，各步驟如圖2所示。按照右手坐標(biāo)系約定，逆向旋轉(zhuǎn)為正。原始坐標(biāo)系統(tǒng)X0-Y0-Z0經(jīng)過圖2的兩個步驟旋轉(zhuǎn)至新的X1-Y1-Z1坐標(biāo)系統(tǒng)，沿方位角旋轉(zhuǎn)后得到的X1軸與過傾斜面的Z1軸垂直。在數(shù)學(xué)中，上述兩步旋轉(zhuǎn)可以用兩個旋轉(zhuǎn)矩陣相乘的結(jié)果來表示：

圖2 坐標(biāo)轉(zhuǎn)換示意圖

(5)

其中旋轉(zhuǎn)系數(shù)分別為：

(6)

(7)

(8)

注意矩陣的右乘。

相應(yīng)的坐標(biāo)逆向旋轉(zhuǎn)為：

(9)

(10)

(11)

(12)

其中:θy和φz分別表示沿y1軸和z0軸逆時針旋轉(zhuǎn)的角度。

(13)

(14)

(15)

中間涉及的相關(guān)參數(shù)如下：

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

利用一階公式并借助鏈?zhǔn)椒▌t可用二次曲面擬合系數(shù)表示以上公式中的二階導(dǎo)數(shù)。

根據(jù)Di[16]提出的理論，我們將式(13)、式(14)、式(15)代入最大正曲率kpos式(24)和最小負(fù)曲率kneg式(25)中。

(24)

(25)

常規(guī)曲率計算是基于最小二乘法近似獲得二次曲面方程，并計算趨勢面系數(shù)，帶入常規(guī)曲率計算方程計算而出。而本次使用的三維坐標(biāo)變換方法所用的處理流程是在計算完趨勢面系數(shù)后，先計算出傾角，方位角，再依次以Z0為旋轉(zhuǎn)軸逆時針旋轉(zhuǎn)該點的方位角角度，以Y0為旋轉(zhuǎn)軸逆時針旋轉(zhuǎn)該點傾角角度，這樣就可以計算新曲面的一階導(dǎo)數(shù)，將一階導(dǎo)數(shù)為零，至此即為新坐標(biāo)系原點，完成三維旋轉(zhuǎn)，最后計算新曲面的二階導(dǎo)數(shù)，帶入曲率計算方程。在相同的情況下，本文方法計算的曲率屬性分辨率更高，效果更好，斷層、裂縫帶顯示更清晰，細(xì)節(jié)更為豐富。

2 模型方法測試

筆者借助三維French模型來驗證本文基于坐標(biāo)變換的旋轉(zhuǎn)體曲率計算方法，在提高曲率屬性分辨率上的有效性(圖3)。

圖3 三維French模型及目的層反射切片

圖4是常規(guī)體曲率和本文,基于坐標(biāo)旋轉(zhuǎn)后體曲率算法所分別得出的最小負(fù)曲率沿層切片與最大正曲率沿層切片。圖4(a)～圖4(b)表示常規(guī)體曲率算法的沿層切片計算結(jié)果，圖4(c)～圖4(d)表示本文基于坐標(biāo)旋轉(zhuǎn)的體曲率算法沿層切片模型計算結(jié)果。如圖4(a)所示，經(jīng)過本文方法計算后所獲得的曲率結(jié)果，將河道顯示地更為清晰，檢測精度也更高，沒出現(xiàn)常規(guī)算法中出現(xiàn)的部分河道缺失等問題。對比圖4(a)和圖4(b)可以發(fā)現(xiàn)，經(jīng)過文中方法計算后，圖4(d)中展示的細(xì)節(jié)更豐富，各細(xì)小構(gòu)造得到了加強(qiáng)，可以更加突出地顯示細(xì)小構(gòu)造的特征，對于我們精確識別微小斷層及裂縫帶很有幫助。

圖4 不同算法的曲率屬性沿層切片

因為實際地震數(shù)據(jù)會包含許多干擾信息(噪聲)，而噪聲對曲率計算是有一定影響的，當(dāng)噪聲過大，還會出現(xiàn)錯誤結(jié)果，因而誤導(dǎo)我們對地下真實構(gòu)造的判斷，所以需要檢驗此方法是否在現(xiàn)實生產(chǎn)中仍能保持相對穩(wěn)定。關(guān)于抗噪性分析(圖5、圖6)，圖5為模型加不同信噪比后的沿層切片圖，本實驗對模型添加隨機(jī)噪聲，添加的噪聲分別為信噪比5.0、10.0、20.0，圖6為當(dāng)模型加不同信噪比后，用本方法處理所得到的最大正曲率切片。對比圖6(a)、圖6(b)、圖6(c)可以發(fā)現(xiàn),模型隨噪聲的加劇逐漸變得很難直觀判斷出構(gòu)造分布，而經(jīng)過處理后的曲率屬性還可以在一定程度上識別出模型真實情況。雖然隨著噪聲的加大，細(xì)節(jié)開始逐漸變得模糊不可分辨，曲率效果也隨之越來越差，但是走向與大致輪廓依舊清楚可現(xiàn)，因此總體來說，經(jīng)過信噪比測試，可以肯定本方法相對穩(wěn)定，具有一定抗噪性。

圖5 模型加噪后沿層切片圖

圖6 加噪后經(jīng)曲率處理的最大正曲率沿層切片圖

3 實際地震資料測試

研究區(qū)位于川東北盆地，處在我國西南地區(qū)四川盆地的東北部，其北邊有秦嶺造山帶，東部為龍門山?jīng)_斷帶，西部為大巴山褶皺帶。研究區(qū)以位于四川A工區(qū)一定范圍的實際地震數(shù)據(jù)為實驗樣本進(jìn)行計算研究。目前已知A強(qiáng)形變地區(qū)分布須二段高產(chǎn)富集帶，富集最有利區(qū)為(海陸相烴源巖+烴源斷裂、陸相斷裂+砂體發(fā)育區(qū))?；煸闯渥ⅲ攲?dǎo)體系與裂縫發(fā)育帶控制“甜點”分布。此處M3井鉆遇良好油氣顯示，海相斷裂供烴，陸相斷裂控縫，海陸斷裂聯(lián)合控制天然氣富集帶分布，井中巖石采樣可見高角度裂縫(圖7)。

圖7 M3高角度裂縫

研究區(qū)須二段屬于裂縫發(fā)育帶控制著“甜點”的分布，該區(qū)域特點是“大面積含氣、局部富集高產(chǎn)”。須家河組二段作為晚三疊世出現(xiàn)的第一個穩(wěn)定層位，其累產(chǎn)氣量近90×108m3，探明產(chǎn)量可占須家河組的70%，此段已測試出了多口高產(chǎn)穩(wěn)定井，獲得工業(yè)氣流，由此可見須二段是整個地區(qū)須家河組天然氣勘探生產(chǎn)的主要產(chǎn)層段。

研究區(qū)裂縫帶走向呈北-東向分布，該方向裂縫形成時期較晚，與現(xiàn)在最大應(yīng)力方向一致，對于裂縫開啟是非常有利的。因為該區(qū)域須二段為大型巖性氣藏，裂縫控制油氣的運(yùn)移與儲集，所以對研究區(qū)裂縫帶精確分析有助于后期的井位部署。結(jié)合以上背景原因，對須二段的部分工區(qū)地震資料進(jìn)行基于坐標(biāo)變換的三維體曲率屬性提取，以獲得更為精確的曲率結(jié)果與裂縫分布顯示。圖7引自中國石化勘探分公司。

圖8(a)～圖8(d)分別為研究區(qū)域的地層T0圖、該區(qū)域的沿層振幅切片圖、最大正曲率沿層切片圖及最小負(fù)曲率沿層切片圖，從中我們可以看出該區(qū)域斷層、裂縫帶宏觀分布，但是很難發(fā)現(xiàn)一些小的斷層和裂縫帶，即常規(guī)方法提取的曲率屬性對細(xì)節(jié)的反映不夠。

圖8 研究區(qū)域地層T0圖、沿層振幅切片圖、最大正曲率沿層切片圖及最小負(fù)曲率沿層切片圖

圖9是使用坐標(biāo)旋轉(zhuǎn)體曲率方法對該區(qū)域所做的最小負(fù)曲率沿層切片和最大正曲率沿層切片。從圖9中可以明顯發(fā)現(xiàn)，最小負(fù)曲率與最大正曲率比常規(guī)方法分辨率更高，能夠顯示更為豐富的細(xì)節(jié)，斷層、裂縫帶分布情況更明顯，解釋結(jié)果與實際工區(qū)區(qū)域地質(zhì)構(gòu)造保持一致，與此地區(qū)裂縫體展布相同。而且通過測井?dāng)?shù)據(jù)可知本實驗地區(qū)M3井裂縫發(fā)育，測井裂縫走向與處理所得曲率屬性切片中的裂縫發(fā)育方向一致，曲率屬性計算預(yù)測結(jié)果與成像測井吻合(圖10)。說明了本方法所計算的曲率屬性較為準(zhǔn)確，證明了本文方法預(yù)測結(jié)果的可靠性。相關(guān)的地震剖面與本工區(qū)測井圖，如圖11、圖12所示。圖11中，原始地震剖面顯示為彩色，曲率屬性顯示為黑白色，從圖11(a)～圖11(d)對比可知，傳統(tǒng)的計算方法在剖面上的顯示較為粗略，相對來說不夠細(xì)致，而本方法計算的曲率屬性在剖面上有更多細(xì)節(jié)體現(xiàn)。

圖9 坐標(biāo)旋轉(zhuǎn)曲率屬性方法沿層切片

圖10 M3二段成像測井裂縫走向與裂縫傾角圖

圖11 曲率屬性與原始地震數(shù)據(jù)疊合的地震剖面局部圖

圖12 M3過井剖面圖

4 結(jié)論

通過模型試算和實際地震數(shù)據(jù)的分析可以得到以下結(jié)論：

1)本方法將三維坐標(biāo)旋轉(zhuǎn)應(yīng)用到曲率屬性計算中，借助構(gòu)造傾角將局部構(gòu)造曲面轉(zhuǎn)換為平面，使目的層在新的旋轉(zhuǎn)坐標(biāo)系中處于水平狀態(tài)，將計算中的一階導(dǎo)數(shù)變?yōu)榱?，既符合工區(qū)實際數(shù)據(jù)的應(yīng)用，又提高了計算效率，使精準(zhǔn)度有所提高。

2)噪聲對本方法也有影響，但在模型加噪后，經(jīng)本方法處理所得的曲率屬性切片其宏觀構(gòu)造形態(tài)依然清晰可見，具有一定程度的抗噪性，證明方法具有實用性與可靠性。

3)通過模型試算與實際地震數(shù)據(jù)分析可以看出，相對于傳統(tǒng)體曲率計算方法，本方法可提取出更多細(xì)節(jié)，更能詳細(xì)地描述地下構(gòu)造，描述裂縫和斷層等的區(qū)域分布情況。