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 引言

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

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

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

1 網(wǎng)格粗化算法及數(shù)據(jù)

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

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

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

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

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

2 網(wǎng)格粗化結(jié)果

2.1 一維空間單一屬性參數(shù)的網(wǎng)格粗化

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

圖1 利用核函數(shù)法中的自適應(yīng)帶寬對一維空間兩個屬性參數(shù)同時(shí)粗化流程圖

圖2 孔隙壓力測井?dāng)?shù)據(jù)（a）及電阻率測井?dāng)?shù)據(jù)（b）

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

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

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

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

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

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

2.2 兩個屬性參數(shù)的同時(shí)粗化

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

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

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

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

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

綜上所述，采用最小帶寬法和最大帶寬法對兩個儲集層屬性參數(shù)同時(shí)粗化得到的模型粗化網(wǎng)格單元數(shù)分別為848和144個?；诤撕瘮?shù)的帶寬，應(yīng)用這兩種不同方法建立1個具有相同網(wǎng)格單元數(shù)量和相同網(wǎng)格單元位置的粗化模型。其中，最大帶寬法中的網(wǎng)格單元尺寸遠(yuǎn)遠(yuǎn)大于最小帶寬法。

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

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

3 結(jié)論

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

符號注釋：