地震層速度模型編輯系統應用研究

劉旭躍

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

0 引言

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

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

圖1 系統架構圖Fig.1 System architecture diagram

1 系統設計

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

圖1中，數據I/O從磁盤中根據層速度SEGY文件的結構和層位文件格式讀取數據或輸出數據；層速度模型繪制：加載的層速度文件經過坐標轉化，將速度值對應于屏幕坐標值，映射為像素顏色值RGB(紅、綠、藍)，轉化為二進制圖像繪制出來。可以選擇層位文件，從文件中讀取該層的速度值，并在原來圖像上繪制層位圖；層速度實時讀取：通過圖像交互功能，鼠標在圖像上移動時，把屏幕坐標轉為圖像坐標，搜索對應的速度值，在顯示區域把速度值顯示出來；速度直方圖顯示：根據存儲在二維數組中的速度值和一維數組中的色標顏色值。用鼠標在速度圖上拖動一個矩形框，在主界面的右下角區域繪制速度直方圖。直方圖由多個小矩形組成，每個矩形代表一個速度區域，鼠標選中矩形，在速度圖上突出顯示相應的速度；層速度編輯：對層速度圖像中速度值編輯。

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

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

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

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

2 關鍵技術

2.1 線性插值

線性插值方法在數學和計算圖形學等學科領域應用非常廣泛,它簡便易用，連續性好[7]。

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

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

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

(1)

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

(2)

那么，用數學方程式考慮就可以表示成為：

y=(1-α)y0+αy1

(3)

或者，

y=y0+α(y1-y0)

(4)

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

2.2 距離反比權重插值

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

(5)

(6)

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

2.3 克里金插值方法

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

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

(7)

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

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

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

擬合模型或空間建模可以稱為變異分析或者是結構分析[11-12]，需要求解的數值，必須放在它自身所處的空間模型中來分析，圍繞指定區域內的所有空間點來計算，利用經驗半變異函數的圖像進行處理。

圖4 經驗半變異函數示例Fig.4 Examples of empirical semi-variograms

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

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

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

如果在計算時發現經驗半變異函數產生的數值與數值模型之間有一些不同，也可認為是存在一些誤差，擬合模型的曲線中，有部分數值大于擬合曲線的值，有部分數值小于擬合曲線的值。可以給出一個參考的空間距離值，待求數值就會全部大于擬合曲線的值，或者給出一個參考的空間距離值，待求數值就會全部小于擬合曲線的值，滿足這種情況時，這個半變異函數模型就達到要求。

3 系統實現

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

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

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

2)根據該文件結構體，找出與速度模型文件相關的層位文件。

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

它們之間的關系采用結構體鏈表表示(圖5)。

圖5 速度模型結構體關系Fig.5 Velocity model structural body relation

在繪圖區域里，用不同顏色代表不同速度值進行對應繪制(圖6)。

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

當鼠標在繪圖區域滑動時，在信息顯示區會根據CDP和Line值顯示速度值(圖7)。

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

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

表1 插值結果

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

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

2)用鼠標在圖形中選擇需要瀏覽的區域。

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

速度編輯時，按如下步驟進行操作：

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

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

清空之后，右擊圖像區域，在彈出的菜單中選擇插值方法，進行插值。插值方法分別是線性插值、反比權重插值、克里金插值。三種插值方法效果如表1所示。線性插值的光滑度不夠，速度值變化趨勢不明顯。距離反比權重插值的冪次越高，結果精確性越高，但計算量越大。克里金插值不僅考慮樣點數據值還考慮周邊鄰近點的位置以及它們之間的關系，結果更為精確。

4 結論

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