LANDMARK系統巖性層位解釋質量控制及調整程序開發

王昭勤，鐘吉太

(大慶鉆探工程公司物探研究院 黑龍江 大慶 163357)

0 引 言

當前大慶油田原油要持續穩產，任務非常艱巨，油田公司和管理局無論是開發環節，還是勘探系統都極為重視。對于我們物探公司從事找油、找氣的排頭兵，也深知自己的重任，作為從事石油勘探一個重要環節的地震資料解釋工作，我們不僅要繼續發揚傳統高精度構造解釋優勢，而且面臨當前的新情況更要注重巖性解釋和地震反演工作。從葡萄花油層、扶楊油層尋找突破，在大慶油田尋找更多的可采儲量，因此，巖性層位解釋的正確與否直接關系到地震資料反演水平和儲層預測的精度，以及后續的油田開發環節。

LandMark系統是美國哈里伯頓公司的地震綜合解釋軟件，在全球廣泛應用。LandMark系統依靠其特色解釋技術(SeisWorks)、相干體技術、可視化技術、屬性分析、分頻處理等技術領先石油勘探開發行業，成為我們地震資料解釋的絕對主流軟件。從常規的構造解釋，到精細的巖性解釋都發揮著重要作用。隨著勘探開發周期縮短。為了加快解釋進度，提高工作效率，巖性層位解釋很多都是在構造解釋基礎之上進行層位漂移，以及層位間插值或者自動追蹤，再進行修改。由于巖性層位彼此之間非常薄，這樣就很容易造成巖性層位間存在竄層“打架”情況，如圖1所示，這就影響到后續的反演處理和巖性預測精度。如果進行巖性層位人工檢查，從平面圖上無法甄別，全工區剖面瀏覽，既費時又費力，難免還有遺漏情況。所以，為了精細解釋，高精度進行儲層預測和地震反演，特開發了巖性解釋層位質量控制及調整程序，目的在于發現巖性層位的質量問題，并進行人工或者自動修正，提高解釋的工作效率，減輕解釋人員的勞動強度。

圖1 巖性層位竄層“打架”情況圖

1 實現過程及程序編制

1.1 數據分析及原理

首先針對OpenWorks系統二維/三維巖性層位數據分別進行數據結構分析，通常的層位數據包括x、y、z即可滿足資料解釋需求，而本程序的要求，需要更多的地震數據信息參與運算，需要包括CDP信息，也就是OpenWorks系統的Single通用格式。程序在設計過程中既要考慮三維工區線、道數不確定性及工區不規則，也要周全二維工區測線長短不一亦即每條線的道數不確定。

其次對“竄層”情況的深入研究，巖性層位竄層通常就是地震數據的某一地震道上的巖性層位，下層的巖性層位漂移到上層巖性層位之上，從T0值上說，下層(H2)時間小于或者等于上層(H1)時間，如圖2所示。通過對下層T0值將超出上層(或重合)的部分自動向下漂移適當的時距(△t)，實現下層巖性層位自動修正[1]，如圖3所示。

圖2 巖性層位漂移圖

用上層加上時距取代原來下層不合理的部分，公式：H2(cdp)=H1(cdp)+△t，同時保持上層數據不變化，下層正常的部分不變化，并產生下層的質量控制記錄文件(QC3D_層名或者QC2D_層名)。質量控制記錄文件是用來檢查竄層的具體線、道位置以及上下層時間和修正后的時間。

圖3 巖性層位自動修正

1.2 程序設計

程序采用非常適合科學計算的Visual Fortran語言編寫，為了方便運行，利于擴展，既可以運行在Windows操作系統上，也可以運行在Linux操作系統上。程序主要包括三部分，第一、第二部分主要用來對巖性層位的上、下層數據進行整理處理，為后續工作做準備，第三部分程序的核心，又來查找竄層部分，并進行修正處理[2]。程序名稱Horizon_QC，流程如圖4所示。

圖4 程序流程圖

1.3 實現過程

1.3.1 二維工區(SORT2D子程序)

首先對上下層巖性層位數據進行分割處理(cut2d)，分割成以線名方式獨立存在臨時性文件，并記錄上下層的線名文件(moni111.dat、moni222.dat)。逐條線進行檢查(compare)，讀取相同線的巖性層位數據，分別賦予不同的一維數組(read_data)，逐道判別是否：VAL=T2(I)-T1(I) 小于零，如果是小于零：T2(I)=T1(I)+SF0 (其中T2下層巖性層位時間，T1上層巖性層位時間)。并將這條線此道相關信息記錄在QC2d質量控制文件[3-5]。

1.3.2 三維工區(SORT3D子程序)

通過給出的上下層巖性層位數據，確定巖性數據有效值的線道范圍(MAX_MIN2d)，用作二維數組下標值上限下限。二維數組用于臨時存儲上下層巖性層位信息，通過逐線(I)逐道(J)，對上下層位數值檢查，發現是否有竄層情況，即VAL=T2(I,J)-T1(I,J)小于零，如果是小于零T2(I,J)=T1(I,J)+SF0，實現對下層校正，并將此點記錄在QC3d質量控制文件，以供檢查[6-7]。

1.4 程序部分核心源代碼

Main horizon_qc()

PRINT*," This program is for searching Horizons reversal "

PRINT*,"writer:wangzhq Email:wangzhq@cnpc.com.cn"

PRINT*,"wwwwwwwwwwwwwwwwwwwwwwwwwww"

print*," ww"

print*," w 1. Processing Lmk2d horizon data w "

print*," ww"

print*," w 2. Processing Lmk3d horizon data w"

print*," ww"

print*," w 0. exit w"

print*," ww"

print*,"wwwwwwwwwwwwwwwwwwwwwwwwwwww"

print*," "

print*,"Pease input a number(0,1,2):"

read*,nu

end

SUBROUTINE SORT3D()

DIMENSION x1(5000,5000),y1(5000,5000),T1(5000,5000)

DIMENSION MCDP1(5000,5000),LINE1(5000,5000)

DIMENSION x2(5000,5000),y2(5000,5000),T2(5000,5000)

DIMENSION MCDP2(5000,5000),LINE2(5000,5000)

character file1*20,file2*20

print*,’ENTER First Horizon(up) FILE(LandMark Single):’

read*,file1

print*,’ENTER Second Horizon(down) FILE(LandMark Single):’

read*,file2

print*,’Enter A FB value(ms):’

read*,FB0

print*,’Enter A shift value(ms):’

read*,sf0

OPEN(UNIT=11,STATUS=’OLD’,FILE=file1)

READ(11,21)m,n

L0=m

L9=m

K0=n

K9=n

close(11)

CALL MAX_MIN2d(L0,L9,K0,K9,FILE1)

CALL MAX_MIN2d(L0,L9,K0,K9,FILE2)

OPEN(UNIT=11,STATUS=’OLD’,FILE=file1)

C read first horizon

print*," >1111111111111111111111111<"

99 READ(11,200,END=88)Lt,kt,Xt,Yt,Tt

LINE1(Lt,Kt)=Lt

MCDP1(Lt,Kt)=Kt

X1(Lt,Kt)=Xt

Y1(Lt,Kt)=Yt

T1(Lt,Kt)=Tt

GOTO 99

C read second horizon

88 print*," >2222222222222222222222222<"

OPEN(UNIT=12,STATUS=’OLD’,FILE=file2)

999 READ(12,200,END=77)Lt,kt,Xt,Yt,Tt

LINE2(Lt,Kt)=Lt

MCDP2(Lt,Kt)=Kt

X2(Lt,Kt)=Xt

Y2(Lt,Kt)=Yt

T2(Lt,Kt)=Tt

goto 999

77 PRINT*," >3333333333333333333333333<"

OPEN(UNIT=13,STATUS=’UNKNOWN’,FILE=’lmk_’//file2)

OPEN(UNIT=14,STATUS=’UNKNOWN’,FILE=’QC3d_’//file2)

WRITE(14,203)’LINE’,’CDP’,

2 ’FirstTime’,"SecondTime","TOP+",sf0

DO 10 I=L0,L9

DO 20 J=K0,K9

IF ((X1(I,J).NE.0).AND.(X2(I,J).NE.0)) THEN

VAL=T2(I,J)-T1(I,J)

IF (VAL.LT.FB0) THEN

A=T2(i,j)

T2(I,J)=T1(I,J)+sf0

WRITE(14,105)i,j,T1(I,J),A,T2(i,j)

ENDIF

fline=LINE2(I,J)*1.0

fmcdp=MCDP2(I,J)*1.0

WRITE(13,201)fline,fmcdp,X2(I,J),

1 Y2(I,J),T2(I,J)

ELSE

IF (X2(I,J).NE.0) THEN

fline=LINE2(I,J)*1.0

fmcdp=MCDP2(I,J)*1.0

WRITE(13,201)fline,fmcdp,X2(I,J),

1 Y2(I,J),T2(I,J)

ENDIF

ENDIF

20 CONTINUE

1.5 程序運行說明

運行平臺為Linux操作系統平臺，運行Horizon_QC，主菜單如圖5所示。

圖5 程序主菜單

在第一級菜單下輸入“1”,二維工區巖性層位處理；輸入“2”是三維工區巖性層位竄層處理。在二級菜單下按程序提示，需要輸入上層和下層巖性層位文件以及時移量(大于0，但不要太大)，程序運行后在當前目錄下產生新的下層巖性文件文件(lmk_原下層文件名)和質量控制文件(QC2d_原下層文件名或QC3d_原下層文件名)[8-9]。

2 應用實例

程序編制后在多個LandMark工區進行了應用，下面以肇源南三維工區為例說明本程序的應用效果。該工區共有7個巖性層位，其中葡萄花油層有p0、p1、p2，扶余油層f0、f1、f2、f3。在用戶完成階段對比后，應用Horizon_QC進行質量檢查。f0、f1是兩個連續巖性層位，運行Horizon_QC，按提示提交上、下層巖性文件名(f0.dat和f1.dat)以及時移量(2.0 ms)，程序自動生成f1層質量控制文件(QC3d_f1.dat)和新的下層巖性文件(lmk_f1.dat)。檢查QC3d_f1.dat文件，如圖6所示，發現其中line432、cdp564處有竄層現象，如圖7所示，直觀很難分辨，局部放大后才能看清f0、f1在cdp564有粘連，從控制文件獲悉下層f1比f0只高出0.01ms。圖8所示是利用本程序f1層自動修正后的情況。

圖6 f1層質量控制文件

圖7 f1在cdp564竄層情況

圖8 修正f1在cdp564竄層情況

3 建議及結論

通過程序的編制及實際解釋工區的應用，彌補了LandMark軟件的不足之處，解決了實際的生產問題，獲得了很好的效果，滿足了用戶解釋要求，提高了解釋人員的工作效率。對于巖性層位的質量監控起到了非常好的作用，具有很好的推廣應用前景。

1)在地震屬性參數提取前，應用本程序(Horizon_QC)進行巖性層位的質量檢查，通過查看質量控制記錄文件QC3D.dat(QC2D.dat)，無問題了再進行后續工作；

2)在為地震反演提供巖性約束層位時，應用本程序Horizon_QC，以新產生的層位文件取代原文件參與地震反演過程。