劉旭躍

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

2.1 線性插值

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

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

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

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

(1)

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

(2)

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

y=(1-α)y0+αy1

(3)

或者，

y=y0+α(y1-y0)

(4)

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

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

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

(5)

(6)

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

2.3 克里金插值方法

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

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

(7)

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

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

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

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

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

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

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

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

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

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

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

沿層速度模型加載時(shí)，首先從工區(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)行對(duì)應(yīng)繪制(圖6)。

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

當(dāng)鼠標(biāo)在繪圖區(qū)域滑動(dòng)時(shí)，在信息顯示區(qū)會(huì)根據(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)系，計(jì)算出區(qū)域內(nèi)的速度值，并進(jìn)行分組顯示。色標(biāo)采用繪圖色標(biāo)，保持速度分組的顏色與原始速度模型的顏色一致，如圖8所示。

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

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

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

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

4 結(jié)論

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