劉旭躍

(中國石化 石油物探技術(shù)研究院,南京 211103)

0 引言

隨著地震勘探技術(shù)快速發(fā)展，勘探區(qū)域逐漸擴(kuò)大，地表條件和地質(zhì)構(gòu)造越來越復(fù)雜。通過收集分析地震波，能夠反演地下構(gòu)造的分布情況，有助于針對性地進(jìn)行油氣勘探，提高開采效率。因此精確地計算出地震波的傳播速度是地震資料處理和解釋的關(guān)鍵問題。地震勘探的采集精度一直在提高，分析面元不斷減小，常規(guī)的速度分析方法不能滿足精細(xì)構(gòu)造成像的要求[1]，利用有效的地震速度分析方法，建立精確的速度模型成為地震勘探的核心問題之一，關(guān)系著整個地震成像的質(zhì)量和最終解釋結(jié)果[2]。地震速度建模方法的相關(guān)研究非?；钴S[3]，就當(dāng)前的實際應(yīng)用而言，基于層析反演理論的速度建模方法仍是主流應(yīng)用技術(shù)[4-5]。

在常規(guī)疊加處理時，都會進(jìn)行速度分析，速度分析時，不論是自動或手動拾取速度點(diǎn)，拾取的速度譜以及沿層速度內(nèi)插后，會存在速度誤差[6]，為減少速度異常值，通過分析速度拾取層的譜信息，對速度值進(jìn)行修整，進(jìn)一步優(yōu)化速度場，更新速度模型，以達(dá)到提高速度精度的目標(biāo)。筆者研究開發(fā)了層速度模型編輯系統(tǒng)，采用Qt面向?qū)ο缶幊陶Z言，讀取繪制層速度模型圖像，通過鼠標(biāo)滑動，實時列出CDP和line對應(yīng)的速度值，圖像的不同顏色代表不同的速度值范圍，直觀地顯示速度變化。在異常值附近選擇指定區(qū)域，就會出現(xiàn)速度矩形列表分析圖，選定某列矩形，就會選定速度區(qū)域，清除后，可采用線性插值或優(yōu)化后的克里金插值法進(jìn)行速度插值，實現(xiàn)精準(zhǔn)去除速度異常值，提高整體速度精度。與國外主流的商業(yè)軟件相比，本系統(tǒng)采用文件式I/O磁盤讀取和自定義數(shù)據(jù)結(jié)構(gòu)體管理，針對自主研發(fā)的復(fù)雜區(qū)域速度建模方法研究而開發(fā)，操作簡單，采用Qt語言開發(fā)，支持跨平臺應(yīng)用，易于移植、擴(kuò)展和后期維護(hù)。

圖1 系統(tǒng)架構(gòu)圖Fig.1 System architecture diagram

1 系統(tǒng)設(shè)計

根據(jù)沿層速度模型編輯的需求，系統(tǒng)開發(fā)的功能有：速度模型加載和輸出模塊、層速度顯示模塊，速度直方圖顯示模塊、層速度編輯模塊、色標(biāo)顯示模塊(圖1)。

圖1中，數(shù)據(jù)I/O從磁盤中根據(jù)層速度SEGY文件的結(jié)構(gòu)和層位文件格式讀取數(shù)據(jù)或輸出數(shù)據(jù)；層速度模型繪制：加載的層速度文件經(jīng)過坐標(biāo)轉(zhuǎn)化，將速度值對應(yīng)于屏幕坐標(biāo)值，映射為像素顏色值RGB(紅、綠、藍(lán))，轉(zhuǎn)化為二進(jìn)制圖像繪制出來?？梢赃x擇層位文件，從文件中讀取該層的速度值，并在原來圖像上繪制層位圖；層速度實時讀?。和ㄟ^圖像交互功能，鼠標(biāo)在圖像上移動時，把屏幕坐標(biāo)轉(zhuǎn)為圖像坐標(biāo)，搜索對應(yīng)的速度值，在顯示區(qū)域把速度值顯示出來；速度直方圖顯示：根據(jù)存儲在二維數(shù)組中的速度值和一維數(shù)組中的色標(biāo)顏色值。用鼠標(biāo)在速度圖上拖動一個矩形框，在主界面的右下角區(qū)域繪制速度直方圖。直方圖由多個小矩形組成，每個矩形代表一個速度區(qū)域，鼠標(biāo)選中矩形，在速度圖上突出顯示相應(yīng)的速度；層速度編輯：對層速度圖像中速度值編輯。

層速度編輯模塊操作流程如下：

圖2 速度編輯流程Fig.2 Speed editing process

圖3 線性插值圖Fig.3 Linear interpolation graph

層速度編輯分為兩個步驟：①清空;②插值。可以選擇3種方式進(jìn)行插值:①線性插值;②權(quán)重線性插值;③克里金插值。色標(biāo)顯示模塊可以選擇色標(biāo)面板中的任意一組，色標(biāo)值變化，直方圖也會改變。清空時，會在插差值顯示區(qū)把清空區(qū)域的色標(biāo)顯示出來，便于對比效果。

2 關(guān)鍵技術(shù)

2.1 線性插值

線性插值方法在數(shù)學(xué)和計算圖形學(xué)等學(xué)科領(lǐng)域應(yīng)用非常廣泛,它簡便易用，連續(xù)性好[7]。

如果有兩個點(diǎn)坐標(biāo)(x0,y0)與(x1,y1),要得到[x0,x1]區(qū)間范圍內(nèi)的x在兩點(diǎn)直線上的值，如圖3所示，計算得到(y-y0)(x-x0)/(y1-y0)(x1-x0)

如果方程兩邊的值是α，則從x0到x的距離與從x0到x1距離的比值就是求取得到的插值系數(shù)。因為x值是之前已經(jīng)知道的，所以從公式(1)計算得到α的值。

α=(x-x0)/(x1-x0)

(1)

同理,α=(y-y0)/(y1-y0)

(2)

那么，用數(shù)學(xué)方程式考慮就可以表示成為：

y=(1-α)y0+αy1

(3)

或者，

y=y0+α(y1-y0)

(4)

采用這種計算方法，求取α值就可以直接知道y。事實上，假如x不是[x0,x1]區(qū)間范圍內(nèi)的值，并且α也不在[0,1]范圍內(nèi)，上述的數(shù)學(xué)公式也是成立的。

2.2 距離反比權(quán)重插值

距離反比權(quán)重插值(Inverse Distance Weighted，IDW)基于樣點(diǎn)相近相似的原理，也是一種廣泛使用的簡單空間插值方法，權(quán)重采用插值點(diǎn)與樣本間的距離得到，然后進(jìn)行加權(quán)平均，由此可見，距離插值點(diǎn)越近的樣本點(diǎn)，對它賦予的權(quán)重值就越大,該算法簡單且時間、空間復(fù)雜度相對都很小[8]。假如二維平面上有一系列的無規(guī)則的離散點(diǎn)，給出它們坐標(biāo)和值為Xj、Yj、Zj(j=1,2,3,…,n)。z點(diǎn)值通過距離加權(quán)值求取得到。Z值的計算公式如下：

(5)

(6)

它是一種比較簡單又有效的數(shù)據(jù)內(nèi)插方法，而且運(yùn)算速度也算較快。距離反比權(quán)重插值的重要影響因子除了權(quán)重距離外，查找半徑和冪次也是它非常重要的影響因子。定長查找是在指定半徑范圍內(nèi)的所有采樣點(diǎn)都要運(yùn)用到柵格單元的插值運(yùn)算中。假設(shè)在之前假定的半徑范圍內(nèi)參與內(nèi)插計算的采樣點(diǎn)個數(shù)比指定的最小數(shù)目小，那么將查找半徑繼續(xù)擴(kuò)大，使得它能夠包含更多的采樣點(diǎn)，從而確保參加計算的采樣點(diǎn)個數(shù)滿足之前指定的最小數(shù)目。

2.3 克里金插值方法

克里金方法是根據(jù)協(xié)方差函數(shù)對隨機(jī)過程或隨機(jī)場進(jìn)行空間建模和預(yù)測的回歸算法[9-10]。它是從地統(tǒng)計學(xué)里面逐漸發(fā)展改進(jìn)而形成的，能夠?qū)值的離散點(diǎn)通過數(shù)學(xué)公式進(jìn)行分析運(yùn)算，它可以被稱為一個線性的數(shù)學(xué)估計系統(tǒng)。只要固有平穩(wěn)隨機(jī)場能用各向同性假設(shè)滿足的話，都可以采用克里金算法計算。它包含多種算法，也可以叫做空間最優(yōu)無篇估計器，經(jīng)過多年發(fā)展，衍生出很多改進(jìn)算法，比如協(xié)同克里金、泛克里金等，近年來隨著人工智能的發(fā)展，它逐漸與其他算法結(jié)合起來，形成新型算法，比如神經(jīng)網(wǎng)絡(luò)克里金、回歸克里金、貝葉斯克里金等。無論怎么改進(jìn)，它的算法是為了精確地產(chǎn)生預(yù)測表面，不斷提高度量預(yù)測的準(zhǔn)確性以及確定性?？死锝鸱椒ㄕf明表面變化的空間相關(guān)性通過假設(shè)采樣點(diǎn)之間的距離或者方向。它確定每個位置的輸出值，采用數(shù)學(xué)擬合過程，將指定數(shù)量的值或者是一定范圍內(nèi)的全部值通過數(shù)學(xué)函數(shù)來計算。而且它不是簡單的計算，要經(jīng)過一個以上的操作才能產(chǎn)生。比如說，它研究方差表面、變異函數(shù)等，克里金方法比較適合數(shù)據(jù)中存在空間相關(guān)距離或者方向偏差，通常應(yīng)用的領(lǐng)域有土壤科學(xué)和地質(zhì)研究中。

克里金方法與反距離權(quán)重法有點(diǎn)相似，都可以對附近的測量值進(jìn)行加權(quán)來得出未測量位置的預(yù)測，通常將加權(quán)總和算法運(yùn)用到這兩種插值方法里面：

(7)

式中：Z(si)為第i個位置處的測量值;λi為第i個位置處的測量值的未知權(quán)重;s0為預(yù)測位置;N為測量值數(shù)。

從上數(shù)公式可知，權(quán)重的值是由計算點(diǎn)在整個隨機(jī)場里的位置關(guān)系而確定的，不單單由計算點(diǎn)自身所處的空間位置和計算點(diǎn)間的距離得到的，它的固有平穩(wěn)隨機(jī)場的數(shù)學(xué)期望都是一樣的，本系統(tǒng)可以選定需要計算點(diǎn)的范圍，把數(shù)值點(diǎn)的 CDP和Line保存起來，得到計算點(diǎn)的空間位置。在各向同性假設(shè)下，計算點(diǎn)的數(shù)學(xué)期望在固有平穩(wěn)隨機(jī)場中，與它自身的空間方位沒有任何關(guān)系。一般采用變異函數(shù)當(dāng)做它的近似值，在運(yùn)算中還可以采用拉格朗日乘數(shù)法來進(jìn)行計算，得到克里金算法的方程組。

采用克里金插值方法處理數(shù)學(xué)計算任務(wù)時，一般都是要完成兩個目標(biāo)，首先是要通過分析數(shù)據(jù)，找到數(shù)據(jù)相互間的聯(lián)系和規(guī)律。然后就是開始對待求解的數(shù)值點(diǎn)通過擬合等過程，進(jìn)行預(yù)估，使之不斷逼近精確值。在算法實現(xiàn)方面，最開始要設(shè)置一系列運(yùn)算中需要采用的變量及函數(shù)值，比如說變異函數(shù)以及協(xié)方差函數(shù)，擬合模型的統(tǒng)計相關(guān)性有關(guān)的一些值可以用創(chuàng)建的函數(shù)來預(yù)估。接著就是通過公式，對需要求解的數(shù)值進(jìn)行預(yù)測。在實現(xiàn)的過程中，算法采用的數(shù)據(jù)是不一樣的，最初是估算要求解值的空間自相關(guān)，得到該值后，隨后再對要求解值預(yù)測。

擬合模型或空間建?？梢苑Q為變異分析或者是結(jié)構(gòu)分析[11-12]，需要求解的數(shù)值，必須放在它自身所處的空間模型中來分析，圍繞指定區(qū)域內(nèi)的所有空間點(diǎn)來計算，利用經(jīng)驗半變異函數(shù)的圖像進(jìn)行處理。

圖4 經(jīng)驗半變異函數(shù)示例Fig.4 Examples of empirical semi-variograms

一般情況下，每個位置的數(shù)值點(diǎn)之間的空間位置都是固定的，每個點(diǎn)是不一樣的，這樣的話，它們的空間點(diǎn)對就有很多種表達(dá)方式。在很短時間內(nèi)要找出空間內(nèi)全部數(shù)值點(diǎn)對，并且在圖像中顯示出來，這將是件困難的事情。比方說，給定一個范圍內(nèi)的數(shù)值點(diǎn)，要找出這個范圍內(nèi)全部數(shù)值點(diǎn)對的平均方差，采用經(jīng)驗半變異函數(shù)來表示的話，如圖4所示，橫向數(shù)值是距離或者表示步長，縱向數(shù)值代表平均半變異函數(shù)值。

空間自相關(guān)是表示變量數(shù)值在一個指定區(qū)域內(nèi)與其他數(shù)據(jù)之間存在的相互依賴關(guān)系。在進(jìn)行空間自相關(guān)量化時候，可以遵循空間位置的邏輯規(guī)律作為評判準(zhǔn)則。認(rèn)為離待求的數(shù)值越近，關(guān)聯(lián)性越高，離待求的數(shù)值越遠(yuǎn)的話，關(guān)聯(lián)性越小。在圖4中表示，在坐標(biāo)軸x軸的最左邊，空間數(shù)值對的距離相對來說會小，那么它們的關(guān)系就更加密切，縱向看的話，就是在y軸的靠下方。在坐標(biāo)軸x軸的最右邊，空間數(shù)值對的距離相對來說會大，那么它們之間的關(guān)系就疏遠(yuǎn)，縱向看的話，就在y軸的靠上方。

將經(jīng)驗半變異函數(shù)生成的數(shù)值進(jìn)行擬合，得到數(shù)值模型。這個過程中，若想精確的對數(shù)值空間進(jìn)行估算和預(yù)測，就要注重半變異函數(shù)建模的環(huán)節(jié)。通過對區(qū)域內(nèi)數(shù)值點(diǎn)的空間屬性估算，利用經(jīng)驗半變異函數(shù)來得到待求數(shù)值的屬性以及它們之間的空間自相關(guān)信息。采用計算連續(xù)函數(shù)來得到經(jīng)驗半變異函數(shù)擬合模型，能夠讓克里金法計算中的克里金方差大于零值，即便不知道所有數(shù)值點(diǎn)的空間屬性的情況下，也能擬合出連續(xù)的函數(shù)曲線，以方便分析和預(yù)測。

如果在計算時發(fā)現(xiàn)經(jīng)驗半變異函數(shù)產(chǎn)生的數(shù)值與數(shù)值模型之間有一些不同，也可認(rèn)為是存在一些誤差，擬合模型的曲線中，有部分?jǐn)?shù)值大于擬合曲線的值，有部分?jǐn)?shù)值小于擬合曲線的值?？梢越o出一個參考的空間距離值，待求數(shù)值就會全部大于擬合曲線的值，或者給出一個參考的空間距離值，待求數(shù)值就會全部小于擬合曲線的值，滿足這種情況時，這個半變異函數(shù)模型就達(dá)到要求。

3 系統(tǒng)實現(xiàn)

本系統(tǒng)基于Linux平臺開發(fā)，采用文件存儲方式。由于系統(tǒng)需要圖形顯示和交互，鑒于Qt跨平臺的C++圖形用戶界面應(yīng)用程序開發(fā)框架，適用于GUI程序開發(fā)，故采用Qt作為開發(fā)環(huán)境。

沿層速度模型加載時，首先從工區(qū)里篩選出符合沿層速度模型文件后綴的數(shù)據(jù)。

1)從中選出沿層速度模型文件。

2)根據(jù)該文件結(jié)構(gòu)體，找出與速度模型文件相關(guān)的層位文件。

3)選出層位文件，讀出層位的速度值。

它們之間的關(guān)系采用結(jié)構(gòu)體鏈表表示(圖5)。

圖5 速度模型結(jié)構(gòu)體關(guān)系Fig.5 Velocity model structural body relation

在繪圖區(qū)域里，用不同顏色代表不同速度值進(jìn)行對應(yīng)繪制(圖6)。

圖6 層速度圖Fig.6 Interval velocity map

當(dāng)鼠標(biāo)在繪圖區(qū)域滑動時，在信息顯示區(qū)會根據(jù)CDP和Line值顯示速度值(圖7)。

圖7 信息顯示圖Fig.7 Information Display Chart

直方圖顯示是根據(jù)速度范圍進(jìn)行繪制。操作步驟如下：

表1 插值結(jié)果

插值完成后，保存速度模型文件

1)速度模型圖繪制完成后， 點(diǎn)擊直方圖按鈕。

2)用鼠標(biāo)在圖形中選擇需要瀏覽的區(qū)域。

3)利用窗口坐標(biāo)與速度值之間的關(guān)系，計算出區(qū)域內(nèi)的速度值，并進(jìn)行分組顯示。色標(biāo)采用繪圖色標(biāo)，保持速度分組的顏色與原始速度模型的顏色一致，如圖8所示。

速度編輯時，按如下步驟進(jìn)行操作：

1)首先選擇速度區(qū)域，右擊圖像區(qū)域，選擇清空菜單。

2)同時會在旁邊空白繪圖區(qū)域里把清空的部分顯示出來。如圖9所示。

清空之后，右擊圖像區(qū)域，在彈出的菜單中選擇插值方法，進(jìn)行插值。插值方法分別是線性插值、反比權(quán)重插值、克里金插值。三種插值方法效果如表1所示。線性插值的光滑度不夠，速度值變化趨勢不明顯。距離反比權(quán)重插值的冪次越高，結(jié)果精確性越高，但計算量越大?？死锝鸩逯挡粌H考慮樣點(diǎn)數(shù)據(jù)值還考慮周邊鄰近點(diǎn)的位置以及它們之間的關(guān)系，結(jié)果更為精確。

4 結(jié)論

根據(jù)速度分析建立精細(xì)、準(zhǔn)確速度場的需求，研究了對沿層速度模型速度值進(jìn)行編輯的插值算法和應(yīng)用工具，開發(fā)了具有加載速度模型數(shù)據(jù)、繪制層速度模型圖像、瀏覽速度值、速度范圍直方圖、速度編輯等功能的軟件系統(tǒng)，實現(xiàn)了有針對性的去除速度異常值，提高速度精度，為后續(xù)地震勘探處理和解釋工作取得有效結(jié)果提供借鑒。