利用核函數法中的自適應帶寬同時對一維空間兩個儲集層屬性參數進行粗化

MOHAMMAD Reza Azad,ABOLGHASEM Kamkar Rouhani,BE HZAD Tokhmechi,MOHAMMAD Arashi,EHSAN Baratnezhad

（1.Shahrood University of Technology,Shahrood 3619995161,Iran; 2.Tarbiat Modarres,Tehran 114-14115,Iran）

0 引言

多孔介質中流體流動性研究對于解決環(huán)境和能源問題非常重要[1-2]。但是從計算工作量的角度，無法求解精細尺度地質模型的流動方程。為了保留儲集層屬性等重要信息并建立流動方程，需要將精細尺度模型轉化為大尺度模型，該過程稱為網格粗化[3-4]。網格粗化可以減少計算工作量并縮短計算時間[5-10]。

許多學者對網格粗化方法進行了改進并取得了非常好的效果，但都具有一定的局限性，并且在屬性粗化過程中存在一些問題[11-18]?；谄骄捣ㄊ莻鹘y(tǒng)的粗化方法[19-21]，該方法不適用于非均質性強的情況。拉普拉斯算子-表皮因子法對三維導流系數粗化結果進行隨機分析，并利用多重邊界法對三維裂隙多孔介質的滲透率進行粗化[17-18]。基于小波的網格粗化法對閾值的選取要求很高[22-24]，在傳統(tǒng)的閾值函數（分為硬閾值函數和軟閾值函數）概念的基礎上[25]，Ge等[26-27]在2015年提出了混合閾值函數的概念，避免了單一使用硬閾值函數或軟閾值函數的缺點。目前，大多數網格粗化方法只能對一維或多維空間1個儲集層屬性進行粗化，無法同時對2個或多個儲集層屬性進行粗化。本文利用核函數法中的自適應帶寬技術，在一維空間同時對2種儲集層屬性進行粗化。在初始模型中，粗化與網格單元的可變性有關，因此需要采用可變帶寬方法，并將核函數法中的帶寬視為系統(tǒng)的變化程度。用于計算可變帶寬的方法主要有：①在動態(tài)或自適應帶寬中，通過設置門檻值，將每個網格單元的粗化表示為網格單元的可變性；②通過誤差值將門檻值限定在所需頻帶或網格單元的數量上。在前述小波變換的網格粗化法中，以二進制格式進行網格粗化，但在此過程中不知道粗化網格單元的數量。前人研究均無法做到從粗化過程的初始狀態(tài)就確定網格單元的數量。本文提出的基于核函數的粗化方法可以從粗化開始就確定模型的粗化網格單元的數量。

本文首先從理論角度介紹核函數法中基于自適應帶寬的網格粗化算法，然后分別給出每個屬性參數的粗化結果。最后，采用兩種方法（兩種方法具有相同的粗化網格單元數量）建立儲集層兩個屬性參數同時粗化的最終粗化模型。

1 網格粗化算法及數據

核函數方法的主要思路是利用鄰域觀測來估計某一點（x點）的分布函數。核函數法是一種靈活估算特定數據分布密度的方法。核函數在結構上與柱狀圖相似，但與柱狀圖的主要區(qū)別在于核函數的密度計算基于x點周圍及具有不同的權重的指定范圍，而在柱狀圖中，密度計算范圍包括x點及具有相同權重的x中心位置，且取決于每個領域。假設X1，X2，…，Xn是種群n的隨機樣本，這個樣本的結果用x1，x2，…，xn表示，則核函數計算的分布函數可用以下方程表示[28]：

核函數通常是概率密度函數，即積分為1，且K域中的x滿足的一個直接結果是核函數密度估算量也是一個概率密度函數。核函數的類型和帶寬是兩個主要參數。在利用核函數方法進行密度函數估計時，核函數本身的作用很小，但是正確選擇帶寬h至關重要[29]?？紤]到系統(tǒng)粗化與網格單元的可變性有關，可以將核函數密度估計法中的帶寬視為系統(tǒng)可變程度的函數。在含有大量信息且具很強可變性的區(qū)域，考慮到任何一種粗化方法，無論其精度和操作方式如何，都應盡量保持系統(tǒng)的主要特性屬性參數，短程帶寬可以達到這個粗化的目的，且該區(qū)域的粗化程度最小。因此，這些區(qū)域將保持小尺度。反之，在較小可變性或平滑變化的區(qū)域，帶寬選擇足夠大的情況下，被合并的原始網格單元最多，這時則由小尺度變?yōu)榇蟪叨取?/p>

可變帶寬可以直接指示系統(tǒng)的可變程度。用動態(tài)方法來確定最佳帶寬，假設隨機樣本為x1，x2，…，xn，使用可變帶寬的核函數方法估算密度的函數公式如下：

此時定義的門檻值或帶寬是可控的，可以從數據中獲得變化量及其最大變化量?；诤撕瘮捣ǖ淖赃m應帶寬對兩個屬性參數的粗化步驟如圖1所示。與傳統(tǒng)的粗化方法相比，該方法的網格粗化是嚴格基于屬性的固有可變性。通過確定可變帶寬，在可變性較大的區(qū)域，選擇較小帶寬，網格單元應盡可能??；相反，在儲集層特征變化較小的地區(qū)，選擇較大帶寬，粗化的網格單元較多。這與核函數帶寬方法的多分辨率特性有很大關系。

本文方法適用于任何類型的測井數據，放射性測井、電測井、聲波測井或機械式測井等數據類型均可。本文以孔隙壓力和電阻率這兩種測井資料為例進行粗化，這兩個屬性參數均為井內測量的兩套不同的數據集。沿450 m井段測量了2 413個點的兩類測井數據（見圖2），根據最大值和最小值之差得到最大門檻值和門檻值的變化。

2 網格粗化結果

2.1 一維空間單一屬性參數的網格粗化

在基于核函數法的自適應帶寬進行網格粗化時，如果兩個相鄰網格單元的差異性小于門檻值或帶寬，則可將這兩個小網格單元合并為一個大網格單元。網格粗化與門檻值關系密切，如果選擇較低的門檻值，則被粗化的網格單元較少；如果選擇較高的門檻值，被粗化的網格單元的數量將會增加，最終網格單元的數量將會大大減少。

圖1 利用核函數法中的自適應帶寬對一維空間兩個屬性參數同時粗化流程圖

圖2 孔隙壓力測井數據（a）及電阻率測井數據（b）

在基于核函數帶寬的網格粗化中，通過定義門檻值可控制計算效率。在本文方法中，門檻值可以用所需網格單元數量的函數進行計算，而網格單元的數量由時間函數表示。通過更改門檻值，誤差值和網格單元的數量將發(fā)生變化。網格單元越多，計算時間越長，計算精度越高。隨著網格單元的變大，計算時間會減少，但計算精度也會降低。使用核函數帶寬的主要思路是在合理的時間以最高的精度進行模擬。提高精度需要增加網格單元的數量，因此在保持網格單元較小的情況下，網格單元過多會大大增大計算量和降低計算效率。計算效率完全取決于網格單元的數量。如果網格單元的數量是固定，例如，如果需要n個網格單元，并且計算時間為t，那么在模擬過程中，當基于可變性進行多分辨率粗化時，就變成n個網格單元，然后運行該模型，以便用網格單元的數量控制原始模型的最大方差并保證最終模型的精度在可接受范圍內。為此，通過考慮核函數評價標準，如均方誤差（MSE）、累積或平方誤差之和（SSE）及定義的門檻值，可以控制粗化網格單元的數量。在核函數粗化方法中，儲集層的可變性可認為是門檻值的函數，通過控制該門檻值，確定模型的最終網格單元數量。這是核變換函數粗化方法和小波變換粗化方法的主要區(qū)別。

在本文算法中，可變帶寬值從0到最大值不等，可等于數據的最小值和最大值之差。對于每個任意帶寬，假設數據x1，x2，…，xn位于帶寬范圍內。如果這些數據的平均值等于則SSE可由以下方程式得出：

因此，對于每個帶寬或門檻值，設一個SSE值和有限粗化帶數量。在對數據的第1個屬性參數進行網格粗化時，門檻值變化范圍為0～1 400，步長為2.5。第2個屬性參數網格粗化的門檻值變化范圍為0～1 000，步長為3。如果選擇較小的門檻值，模型會保留更多的儲集層信息，粗化過程將從更多的非均質區(qū)域開始。零門檻值將使整個網格單元保持細粒度且不會進行粗?；Ｗ畲箝T檻值可將網格單元整體收縮到一個均勻區(qū)域，整個模型被轉換成一個網格單元。網格粗化過程中的另一個關鍵點是設置最優(yōu)門檻值或適當的帶寬，對于提高計算效率是非常有效的。為此，提出了3種不同的方法：SSE變化-帶寬關系法、連續(xù)SSE變化法、SSE可變性與帶寬關系曲線和網格單元數量-變化性與帶寬關系曲線相交法。利用SSE差分法確定最優(yōu)門檻值或最優(yōu)帶寬。圖3顯示了第1個屬性參數（孔隙壓力）的粗化結果，當門檻值為50時，SSE值為86，網格單元數量為940；當門檻值為80時，SSE值為176，網格單元數量為682；當門檻值為100時，SSE值為331，網格單元數量為554。結果表明，最優(yōu)門檻值為80。圖4顯示了第2個屬性參數（電阻率）的粗化結果，當門檻值為15時，SSE值為14，網格單元數量為800；當門檻值為30時，SSE值為48，網格單元數量為590；當門檻值為80時，SSE值為192，網格單元數量為391。結果表明，最優(yōu)門檻值為30。

圖3 第1個屬性參數（孔隙壓力）在不同門檻值下的粗化結果

圖4 第2個屬性參數（電阻率）在不同門檻值下的粗化結果

基于核函數法中的自適應帶寬進行網格粗化取決于計算時間和模型精度。例如，如果粗化的目的僅僅是降低計算成本，可以選擇較大的帶寬；反之，如果粗化的目的是提高模型的精度，則應選擇較小的帶寬。

2.2 兩個屬性參數的同時粗化

同時對兩個或多個儲集層屬性參數進行網格粗化的結果更加可靠，這種方法的目的是采用兩種不同方法獲得兩個相同數量的網格單元，然后根據這兩個屬性參數建立流體流動方程。本文首先計算儲集層兩個屬性參數的最優(yōu)帶寬并得到粗化結果，根據最小帶寬法和最大帶寬法對兩個屬性參數同時粗化獲得最終的粗化模型，同時保證這兩種方法的網格單元數和網格位置都相同。

如前所述，第1個屬性參數的最優(yōu)門檻值是80，粗化網格單元的數量為682；第2個屬性參數的最優(yōu)門檻值是30，粗化網格單元的數量為590。在同時對這兩個屬性參數進行粗化時，需要保證按比例粗化的網格單元數量相等且網格單元的位置也相同。在最小帶寬法中，以兩個屬性參數中帶寬較小的一個作為粗化標準。通過比較兩個屬性參數的粗化網格單元，將其中較小的一個作為最優(yōu)網格進行計算。顯然，兩個屬性參數最終粗化模型的網格單元數大于單一屬性參數的數量。圖5為采用最小帶寬法對兩個屬性參數同時粗化的結果，最終的粗化網格單元數為848，原始模型的網格單元有2 413個。應用最小帶寬法時，第1個屬性參數的粗化誤差為390，第2個屬性參數的粗化誤差為187。在兩個屬性參數的粗化模型具有相同的網格單元數量和相同的網格位置。

圖5 最小帶寬法第1個屬性參數（a）和第2個屬性參數（b）粗化結果

采用最大帶寬法時，使用兩種屬性參數中帶寬較大的一個作為粗化標準。通過比較兩種屬性的粗化網格單元，將較大的一個作為兩個屬性參數同時粗化的最優(yōu)網格進行計算。顯然，兩個屬性參數的最終粗化模型的網格單元數量小于單一屬性參數的數量。圖6為采用最大帶寬法對兩個屬性參數同時粗化的結果，最終的粗化網格單元數量為144。應用最大帶寬法時，第1個屬性參數的粗化誤差為1 000，第2個屬性參數的粗化誤差為1 532。

圖6 最大帶寬法第1個屬性參數（a）和第2個屬性參數（b）粗化結果

綜上所述，采用最小帶寬法和最大帶寬法對兩個儲集層屬性參數同時粗化得到的模型粗化網格單元數分別為848和144個。基于核函數的帶寬，應用這兩種不同方法建立1個具有相同網格單元數量和相同網格單元位置的粗化模型。其中，最大帶寬法中的網格單元尺寸遠遠大于最小帶寬法。

為了驗證方法的可靠性，將粗化結果與測井和試井解釋結果進行對比（見表1），在某些點上，估算壓力和粗化壓力相等。此外，實際壓力值與基于粗化模型得到的壓力值之間的差異不大，可以接受。

表1 不同方法得到的孔隙壓力值對比

3 結論

本文介紹了一種利用核函數法中的自適應帶寬同時對一維空間兩個儲集層屬性參數進行粗化的方法。該方法既減少了網格數量，同時也保留了原始精細模型的主要非均質性特征。本文方法采用最小帶寬法和最大帶寬法對兩個屬性參數同時進行粗化，保留了屬性參數高可變性的精細尺度，并對低可變性區(qū)域進行粗化，由此建立兩種屬性參數的最終粗化模型。最小帶寬法和最大帶寬法的粗化誤差率由最終的網格單元數決定。從計算的角度來看，帶寬粗化法得到了具有多個表示對數網格單元的非均勻對數，且最小帶寬法的粗化誤差小于最大帶寬法。

符號注釋：