史敬華,，李少華*，束青林，張以根，于金彪，盧昌盛

1長(zhǎng)江大學(xué)地球科學(xué)學(xué)院，湖北 武漢

2中石化勝利油田分公司，山東 東營

1. 引言

儲(chǔ)層隨機(jī)建模技術(shù)在現(xiàn)代油藏描述中得到越來越廣泛的應(yīng)用[1] [2]，該項(xiàng)技術(shù)能夠更好地刻畫儲(chǔ)層的非均質(zhì)性和定量描述儲(chǔ)層中的不確定性[3] [4] [5]。儲(chǔ)層隨機(jī)建模方法按照模擬對(duì)象的不同，可以分為基于像元(Pixel)、目標(biāo)(Object)和樣式(Pattern)三類[6] [7] [8]。基于像元的建模方法在模擬的過程中是一個(gè)網(wǎng)格一個(gè)網(wǎng)格逐個(gè)的模擬，模擬的順序是隨機(jī)的，因此建立的模型通常會(huì)含有一部分“噪點(diǎn)”，也就是分布很離散、缺乏地質(zhì)意義的一些相對(duì)孤立的網(wǎng)格點(diǎn)或少量網(wǎng)格點(diǎn)組合。針對(duì)這一問題，加拿大學(xué)者Clayton提出了MAPS 平滑方法[9]，該方法通過一個(gè)加權(quán)平均窗口對(duì)初始模擬結(jié)果進(jìn)行平滑，能夠過濾掉認(rèn)為不符合地質(zhì)認(rèn)識(shí)的孤立的網(wǎng)格點(diǎn)。國內(nèi)學(xué)者尹艷樹在此基礎(chǔ)上提出了一種基于信息度的平滑方法[10]，在權(quán)值確定中考慮了條件數(shù)據(jù)的影響，進(jìn)一步完善了該平滑方法。在商業(yè)化軟件中也增加了平滑處理的模塊。在多數(shù)情況下，平滑方法能夠取得不錯(cuò)的效果。當(dāng)模型中含有薄的呈線性(二維)或片狀(三維)分布的地質(zhì)體時(shí)，例如薄的泥巖夾層，這些地質(zhì)體對(duì)滲流又有重要的影響[11] [12] [13] [14]，需要在模型中進(jìn)行表征。MAPS 方法會(huì)把這種地質(zhì)體完全平滑掉，而不僅僅是把噪點(diǎn)去掉。針對(duì)這一問題，對(duì)平滑方法MAPS 進(jìn)行了改進(jìn)，在存在薄的地質(zhì)體時(shí)，針對(duì)性地設(shè)計(jì)平滑的窗口和權(quán)值，取得了較好的效果。同時(shí)對(duì)于邊界上的噪點(diǎn)也進(jìn)行了處理，原始的MAPS 方法沒有考慮這一問題。

2. 問題的提出

序貫指示模擬(SISIM)、多點(diǎn)地質(zhì)統(tǒng)計(jì)學(xué)SNESIM 等基于像元的建模算法，在模擬離散變量(如巖性、巖相等)過程中，采用的是根據(jù)隨機(jī)路徑確定的網(wǎng)格點(diǎn)逐點(diǎn)進(jìn)行模擬。因?yàn)槭菃蝹€(gè)網(wǎng)格單元的模擬，容易導(dǎo)致模擬結(jié)果出現(xiàn)相對(duì)離散的點(diǎn)，這些離散的點(diǎn)不符合地質(zhì)認(rèn)識(shí)，是所謂的“噪點(diǎn)”。Clayton 教授提出了一種模型平滑方法MAPS 可以在一定程度上解決該問題。在MAPS 方法中，主要采用了一種平滑窗口，如二維情況可以采用3 × 3、5 × 5 窗口對(duì)模型進(jìn)行平滑處理。圖1 為平滑前后的對(duì)比，一般情況下MAPS可以取得較好的效果。但是在某些特殊情況下，現(xiàn)有的平滑方法不適用。如圖2 所示，該圖為阿爾伯塔大學(xué)地質(zhì)統(tǒng)計(jì)學(xué)授課中的平滑原理示意圖，采用3 × 3 的模板，權(quán)值分布采用的是各向同性的方法，即中心點(diǎn)權(quán)值最大，距離中心點(diǎn)距離相同的網(wǎng)格的權(quán)值一樣。當(dāng)平滑窗口處于如圖中紅色虛線框所示位置時(shí)，中心點(diǎn)所屬相類型的概率采用加權(quán)平均的方法。例如平滑窗口共有9 個(gè)網(wǎng)格，權(quán)值之和為15，而屬于代碼1 的有兩個(gè)網(wǎng)格，其權(quán)值分別為1 和3，因此中心點(diǎn)屬于代碼1 的概率為4/15 = 0.27。同理，屬于代碼2 的概率為(2 + 1 + 2 + 2)/15 = 0.47，屬于代碼3 的概率為(1 + 2 + 1)/15 = 0.27。因此，平滑后紅色虛線框內(nèi)的中心網(wǎng)格的屬性將用概率最大的代碼2 所替代，這樣就去除了離散的點(diǎn)。在研究中發(fā)現(xiàn)，當(dāng)平滑窗口位于左下角時(shí)，如圖中藍(lán)色虛線框所示，則中心點(diǎn)屬于代碼2 的概率為8/15 = 0.53，而屬于代碼3 的概率為7/15 = 0.47，因此中心點(diǎn)會(huì)被代碼3 所替代。繼續(xù)平移平滑窗口(向右平移一個(gè)網(wǎng)格)，下一個(gè)平滑窗口中心點(diǎn)的綠色(代碼3)也將被白色(代碼2)所替代。如此重復(fù)下去，整個(gè)綠色條帶都將被白色所替代。通過圖2 可以直觀的看出，綠色條帶并不是孤立的，連續(xù)性很好，只是因?yàn)槌示€性分布，所以這種平滑方法針對(duì)線性分布的地質(zhì)體是不適用的。

Figure 1. Comparison before and after smoothing 圖1. 平滑前后對(duì)比

3. 方法的改進(jìn)

改進(jìn)處理方法主要針對(duì)MAPS 方法中平滑窗口進(jìn)行重新設(shè)計(jì)，針對(duì)特殊的地質(zhì)情況設(shè)計(jì)合理的平滑窗口和權(quán)值。原始MAPS 算法中平滑窗口為各向同性，而泥巖夾層的分布主要呈線性分布(二維)或薄片狀分布(三維)，不同方向的連續(xù)性差別較大。采用各向同性的平滑窗口不僅把孤立的點(diǎn)平滑掉，而且也會(huì)把線性分布的泥巖平滑掉。以薄層泥巖的平滑處理為例說明，將各向同性平滑窗口改進(jìn)為長(zhǎng)方形的各向異性平滑窗口，如圖3 所示，窗口不再是規(guī)則的正方形網(wǎng)格，而是根據(jù)地質(zhì)體的空間分布特征來設(shè)置，如本例中水平方向設(shè)置為5 個(gè)網(wǎng)格，而在垂直方向上設(shè)置為3 個(gè)網(wǎng)格，并且權(quán)值不再是各向同性分布，例如本例中中間網(wǎng)格的權(quán)值明顯高于上下層網(wǎng)格的權(quán)值，這種窗口適合于圖2 中呈線性分布地質(zhì)體的平滑。在實(shí)際應(yīng)用中需要根據(jù)地質(zhì)體的特點(diǎn)來設(shè)計(jì)合理的平滑窗口。

Figure 3. 5 × 3 smooth window 圖3. 5 × 3 平滑窗口

使用改進(jìn)后的平滑窗口對(duì)圖2 進(jìn)行平滑，當(dāng)平滑窗口位于圖中黃色虛線時(shí)。平滑窗口共有15 個(gè)網(wǎng)格，權(quán)值之和為59。屬于代碼1 概率為1/59 = 0.0169，屬于代碼2 的概率為9/59 = 0.1525，屬于代碼3 的概率為49/59 = 0.8305。平滑后圖中綠色相帶就保留下來了，而圖中紅色虛線內(nèi)的離散點(diǎn)按照此平滑窗口也會(huì)正確的被白色所替代。提高平滑窗口中心點(diǎn)及其左右兩邊網(wǎng)格的權(quán)重可以滿足了線性分布巖相類型的平滑處理。

下面采用通用模板進(jìn)一步說明改進(jìn)平滑窗口的工作原理以及需要注意的問題：假設(shè)巖相模型中單層泥巖分布代碼如下，0 代表背景或者其他巖相類型，1 代表泥巖夾層，紅色表示泥巖分布的兩端，1 的數(shù)量為L(zhǎng)，表示泥巖長(zhǎng)度。

Figure 4. Distribution of single mudstone layer 圖4. 單層泥巖分布

Figure 5. General 3*m smooth window 圖5. 通用3*m 平滑窗口

使用圖5 所示3*m 的平滑窗口去平滑圖4 所示泥巖夾層。當(dāng)窗口中心點(diǎn)K 移動(dòng)到圖4 中左邊紅色1處時(shí)，屬于泥巖的權(quán)值為K + K1 + … + Kn，屬于其他相的權(quán)值為K1 + … + Kn + sum(‘1’)，sum(‘1’)表示平滑窗口中所有權(quán)值為1 的和。為保證泥巖夾層兩端(圖4 中紅色1)不會(huì)被其他相類型(圖4 代碼0)替代，屬于其他相的權(quán)值應(yīng)該小于屬于泥巖的權(quán)值，所以平滑窗口應(yīng)該滿足K > sum(‘1’)，這樣平滑過程中泥巖夾層兩端不會(huì)被平滑掉，泥巖夾層中間也不會(huì)被平滑掉，這樣整個(gè)泥巖夾層就完整保留下來。圖5 所示平滑窗口長(zhǎng)度為m，圖4 所示泥巖夾層長(zhǎng)度L，為保證窗口中心點(diǎn)K 移動(dòng)到圖4 中左邊紅色1 處時(shí)，窗口邊緣點(diǎn)Kn 在圖4 右邊紅色1 的左邊，需要滿足n ≥ L，其中m = 2*n + 1，這樣上述權(quán)值才能正確計(jì)算，泥巖夾層才能保留下來。所以若想在平滑過程中不保留長(zhǎng)度為L(zhǎng) 以下的泥巖夾層，則讓n < L，并且Kn > K-sum(‘1’)，這樣在移動(dòng)窗口的過程中所有泥巖夾層的權(quán)值都小于背景相權(quán)值而被平滑掉。

原始MAPS 方法在平滑過程中當(dāng)模板中心點(diǎn)K 位于模型的左上角時(shí)，平滑窗口的左上部分在泥巖夾層模型外部，只有中心點(diǎn)K 右面及下面平滑窗口權(quán)值參與計(jì)算。這樣平滑窗口中的權(quán)值沒有被完全使用，會(huì)損失部分信息。在改進(jìn)方法中對(duì)原始泥巖夾層模型進(jìn)行向外擴(kuò)充，比如使用3 × 5 平滑窗口對(duì)網(wǎng)格數(shù)量為200 × 200 的泥巖夾層模型進(jìn)行平滑，事先在模型上面和下面各增加一層背景相，在模型左面及右面各增加兩層背景相，整個(gè)模型的網(wǎng)格數(shù)量為202 × 204。這樣平滑窗口在202 × 204 模型網(wǎng)格內(nèi)部平滑就不會(huì)出現(xiàn)平滑窗口溢出到模型外部的情況，使得平滑窗口中每個(gè)權(quán)值都得到計(jì)算，避免了模型邊界處的平滑誤差。

4. 改進(jìn)前后的對(duì)比

在老油田SISIM 方法經(jīng)常被用于泥巖夾層模型的建立，一般來說由于網(wǎng)格大小的限制以及泥巖夾層厚度分布的特點(diǎn)，泥巖夾層一般呈條帶狀分布。如圖6 所示例子，網(wǎng)格大小為10 × 1 米，泥巖夾層的厚度為1～2 米，平均長(zhǎng)度為200 米，泥巖比例為10%。采用SISIM 模擬，變差函數(shù)設(shè)置為水平方向變程為200 米，垂直方向2 米，理論模型為球狀模型。圖6 左邊為模擬的結(jié)果，可以看出夾層呈水平條帶分布，局部存在一些離散的點(diǎn)或是橫向延伸很短的線。右邊為泥巖夾層的長(zhǎng)度分布直方圖。在統(tǒng)計(jì)泥巖延伸長(zhǎng)度的過程中，由于相鄰單層網(wǎng)格的泥巖會(huì)出現(xiàn)在垂向上疊加而形成較長(zhǎng)的泥巖夾層的情況，單獨(dú)統(tǒng)計(jì)某一層網(wǎng)格上的泥巖長(zhǎng)度與實(shí)際泥巖長(zhǎng)度不符，所以本文設(shè)置統(tǒng)計(jì)泥巖長(zhǎng)度規(guī)則如下：所有能夠與某一泥巖網(wǎng)格在上下左右或?qū)巧舷噜彽哪鄮r網(wǎng)格都屬于同一個(gè)泥巖夾層，網(wǎng)格在統(tǒng)計(jì)過程中不會(huì)被重復(fù)計(jì)算。統(tǒng)計(jì)夾層長(zhǎng)度算法先遍歷整個(gè)模型，提取出所有獨(dú)立的相互不連通的泥巖夾層，再計(jì)算每個(gè)泥巖夾層的最大長(zhǎng)度。計(jì)算的原始模型中泥巖長(zhǎng)度平均值為140.4 米，比輸入的期望值200 米(水平方向變程)小了近1/3，主要原因是因?yàn)榇嬖谝恍╇x散的噪點(diǎn)。采用MAPS 方法平滑，平滑樣板3 × 3，平滑結(jié)果如圖7 所示，泥巖長(zhǎng)度平均值為92.6 米，小于平滑前的泥巖長(zhǎng)度平均值。對(duì)比圖6 和圖7 可以看出，采用MAPS平滑泥巖夾層模型有如下特點(diǎn)，對(duì)于單層的泥巖夾層，MAPS 方法會(huì)將其完全平滑為背景相，對(duì)于離散噪點(diǎn)MAPS 方法去除不完全，處理后不僅僅導(dǎo)致泥巖長(zhǎng)度變短，而且很多泥巖夾層都被平滑掉了，也大大改變了原始的泥巖百分比。

針對(duì)MAPS 方法平滑存在的不足，采用上節(jié)提出的改進(jìn)方法，按照改進(jìn)平滑窗口設(shè)置規(guī)則，如果需要平滑掉長(zhǎng)度為30 米及30 米以下的泥巖夾層，設(shè)置3 × 9 大小的模板[1 1 1 1 1 1 1 1 1; 4 5 6 7 20 7 6 5 4; 1 1 1 1 1 1 1 1 1] 對(duì)泥巖夾層模型(圖6)進(jìn)行平滑，迭代兩次。平滑后的結(jié)果如圖8 所示?？梢钥闯銎交竽鄮r夾層最短長(zhǎng)度為30 米，所有30 米及以下的泥巖夾層噪點(diǎn)都被平滑掉了。新平滑方法保留了原始泥巖夾層的分布特點(diǎn)，針對(duì)噪點(diǎn)及認(rèn)為不合理的泥巖短線進(jìn)行了剔除。統(tǒng)計(jì)直方圖顯示平均泥巖長(zhǎng)度為164 米，相對(duì)于MAPS 處理后模型的140 米，更接近期望值200 米。

Figure 6. Interlayer distribution model simulated by SISIM 圖6. SISIM 模擬夾層分布模型

Figure 7. Smooth effect of MAPS method 圖7. MAPS 方法平滑效果

Figure 8. Smoothing results with mudstone length greater than 30 m 圖8. 泥巖長(zhǎng)度大于30 米的平滑結(jié)果

Figure 9. Smoothing results with mudstone length greater than 40 m 圖9. 泥巖長(zhǎng)度大于40 米的平滑結(jié)果

該算法設(shè)計(jì)靈活，可以根據(jù)實(shí)際資料的情況設(shè)置不同的下限值，例如設(shè)置3 × 11 大小的模板[1 1 1 1 1 1 1 1 1 1 1; 4 5 6 7 8 24 8 7 6 5 4; 1 1 1 1 1 1 1 1 1 1 1]對(duì)泥巖夾層模型(圖6)進(jìn)行平滑，可以平滑掉長(zhǎng)度40米及40 米以下的泥巖夾層，平滑結(jié)果如圖9 所示。對(duì)比圖8、圖9 和圖6 可以看出，新的平滑方法可以消除給定長(zhǎng)度的噪點(diǎn)，平滑處理后的泥巖夾層的平均長(zhǎng)度更接近模擬前預(yù)設(shè)的期望值，也就是說新的平滑方法效果明顯。通過改變平滑窗口的長(zhǎng)度能夠保留長(zhǎng)度在指定值以上的泥巖夾層，改進(jìn)算法更具有針對(duì)性。

5. 結(jié)論

基于像元的建模方法如SISIM、SNESIM 等經(jīng)常被用于建立離散變量的模型，如微相、巖相、巖性模型等。由于算法本身的特點(diǎn)，會(huì)導(dǎo)致模擬結(jié)果中會(huì)出現(xiàn)一些不符合地質(zhì)認(rèn)識(shí)的離散的“噪點(diǎn)”，為了使模型更合理，通常需要做平滑后處理。Clayton 教授提出的MAPS 平滑方法在多數(shù)情況下均能取得較好的效果，但是當(dāng)存在薄的地質(zhì)體時(shí)，例如薄的泥巖夾層，該平滑方法會(huì)把線性或片狀分布的夾層完全平滑掉。針對(duì)這種情況，本文提出了改進(jìn)方法，根據(jù)地質(zhì)體的特點(diǎn)設(shè)計(jì)平滑窗口，并根據(jù)需要平滑掉的最短泥巖長(zhǎng)度設(shè)計(jì)相應(yīng)的權(quán)值，改進(jìn)后的平滑方法不僅能夠去除“噪點(diǎn)”，而且能夠克服原始MAPS 方法把薄的地質(zhì)體平滑掉的不足。通過對(duì)比分析發(fā)現(xiàn)，改進(jìn)后方法得到的地質(zhì)模型更接近地質(zhì)真實(shí)情況。本文給出的算例是二維模型，提出的算法可以很容易擴(kuò)展到三維空間。

致 謝

本文是國家科技重大專項(xiàng)任務(wù)(辮狀河儲(chǔ)層構(gòu)型建模算法模塊編寫及知識(shí)庫應(yīng)用系統(tǒng)研制，編號(hào)2016ZX05011-001-004)的部分研究成果。