地質統計學變程優化研究

郝紅瑞潘懋劉振曹凱

（北京大學地球與空間科學學院北京100871）

地質統計學變程優化研究

郝紅瑞潘懋劉振曹凱

（北京大學地球與空間科學學院北京100871）

三維地質屬性建模需要采用各種地質屬性的樣本數據，經過變異函數分析以及合適的插值算法，對未知區域進行相應插值，以獲取研究區內的屬性分布。為了構建更為精確的屬性模型通常在建模過程中需要一定的人機交互。但在以往屬性建模中，變程的有效選取一般根據屬性建模的結果予以不斷調整來獲取，而沒有一個參考的標準，這就給屬性建模的效率和精度造成很大影響。針對這一問題，嘗試進行地質統計學變程優化研究，尋求有效變程同相關影響因素之間的理論模型，從而提高屬性建模的效率和精度。

屬性建模；變異函數；有效變程；屬性變化；地質體規模；序貫高斯

Class NumberP628

1 引言

三維地質屬性建模的目的是獲得能夠表達建模區域內任意一點屬性值的模型，原始數據來源于采樣點樣品數據的屬性值，因此屬性建模的核心內容是從已有樣品點的屬性值得到建模區域內任意一點屬性值所采用的方法［1］。即在原有樣品點的基礎上根據一定的插值方法獲得整個研究區域的屬性特征。

最適用于三維地質屬性建模的方法為地質統計學方法。地統計學方法中引入變異函數的概念來量化地質體內部的空間相關性，考慮研究區域的空間變異性，從而得到更為精準的屬性模型。一般來說，變異函數可以由兩個方面來獲取：可以通過地質學專家對研究區域的專業認識來給予指導；也可以根據獲取到的原始數據進行分析得到實驗變異函數，再利用最小二乘法［2～4］、線性規劃［5］、加權線性規劃［6～7］等方法進行擬合得到表征變異性的各個參數。

但是，為了得到最佳屬性模型需要不斷調整變異函數中的各個參數，其中最為重要的是表征空間變異性最大距離的參數，即變程，參數不斷調整的過程就會影響屬性建模整體的速度，調整不當亦會影響屬性建模的精度。因此，本文旨在探索地質統計學中的變程優化研究，并且給定變程與相關影響因素之間的理論模型，以期對于快速精確地構建屬性建模提供一定程度的指導。

2 理論背景

地質統計學由法國統計學家G.Matheron創立，是一門以區域化變量為基礎，借助變異函數來研究空間變異性和空間結構性的學科，最早應用于地質學領域。在運用地質統計學方法時需要用到一個重要工具——變異函數。變異函數代表兩點之間屬性值方差的二分之一，用來定量表征兩點之間的相互關系，一般用塊金值、基臺值和變程三個參數來表達。塊金值指距離無限接近0時的變異函數值；基臺值是變異函數在距離大于變程時的變異函數值，用來刻畫區域變量的總體變異性［8］。其中最為重要的參數為變程a。

變程是指變異函數在達到基臺值時所對應的距離，即自相關距離。當距離h≤a時，兩點間的屬性值存在相關性，并且該相關性隨著h的增大而呈現逐漸減小的趨勢；而當h＞a時，兩點不再具有相關性。即變程a是區域變量從空間相關狀態轉向不相關狀態的轉折點。因此變程是反映區域變量變化程度的重要參數。

變程的獲取通常需要三步。首先，對采樣數據進行統計分析處理；然后應用常用的變異函數理論模型（球狀模型、指數模型、高斯模型、空穴模型等）進行擬合，從而得到研究區域整體的變量相關關系；最后，從理論模型中得到需要的變程。但該變程只是一個理論變程，在實際應用中存在一定偏差。因此，對于構建合理的屬性模型，選取一個最優變程非常關鍵，該變程稱為有效變程。如何高效獲取有效變程對于快速構建合理的屬性模型起著至關重要的作用，這也是本文的研究重點。

3 有效變程獲取

在研究空間變異性時通常先需要進行變異函數的擬合，然后對重要參數進行相應調整，以便更準確地刻畫變量的空間變異情況。本文的有效變程獲取分為擬合變異函數和校正有效變程兩部分。具體過程如下。

3.1 擬合變異函數

在進行插值計算前，需要一個確定的理論變異函數，用來量化數據間的相關關系。一般采取的方法是：首先設置合適的單位滯后距，利用實驗變異函數的計算公式計算對應于每個滯后距的變異函數值，即得到一系列散點圖；然后，任意選取一種理論模型，其中本文中以球狀模型為例說明，采用線性規劃法求解變異函數模型中的各個參數，從而得到理論變異函數。具體過程如下：

1）設定合適的單位滯后距和點對數量，每個滯后距對應一個實驗變異函數值，從而得到一組滯后距h和實驗變異函數γ*() h的對應散點圖。

2）選用一種理論模型——球狀模型進行擬合，數學表達式［9～10］為

擬合變異函數即根據實驗數據點對進行計算進而求取未知參數C0，C和變程a。將球狀模型簡化為線性函數。即令：

則球狀模型簡化為：y=k1+k2·h+k3·(-h3)。

將求取到的實驗數據點對(hi，γ(hi))帶入上述函數中，并且令

3.2 校正有效變程

由于探索空間變異程度時，首先要利用原始數據計算實驗變異函數，然后進行相關的分析和擬合，進而推斷出整個研究區域的空間變異性。其中，實驗變異函數計算公式如下：

表示所有位置之差為向量h的兩個樣品點屬性值之差的平方的均值的一半。由式（2）可以看出，對應不同的屬性參數（如：孔隙度、滲透率、含油飽和度等）屬性變化不同，求取到的變異函數亦會不同。因此，有效變程與屬性的變化范圍有關。

另外，除了屬性變化，控制實驗變異函數求取的還有一個重要參數：滯后距h，因此h的選取亦會影響到實驗變異函數的準確性。通常h的選取原則為：滯后距的數量nh與滯后距h的乘積nh·h大約為被考慮方向上地質體規模的一半。因為滯后距大于地質體規模的一半時不能保證地質體中央的數據能被應用到，進而變異函數變得不穩定并且不能代表整個地質體。而且，這樣長距離的變異函數在后續的地質統計學建模中也通常不需要。因此，變異函數隨地質體規模的變化而變化。換句話說，有效變程與所研究地質體的規模，特別是地質體沿水平方向延伸的距離有關。為方便起見，下文均采用地質體規模代替該距離。

為尋求有效變程同屬性變化和地質體規模兩個參數之間的關系，得到各參數之間的理論模型，本文采用實驗分析的方法進行探索，應用大量實驗數據并且在不同的研究區域內進行模擬，并對結果進行相關回歸分析得到最終理論模型。

3.2.1 建立有效變程理論模型

首先，選取一組測試數據進行相應試驗，實驗證明得到的有效變程與擬合變異函數得到的理論變程存在一定的差距，因此求取有效變程的理論模型是非常必要的。然后，對多組不同數據進行實驗，并且運用回歸擬合的方法得到理論模型。最后，對得到的理論模型進行驗證，證明該方法的可行性。下面分別從研究思路和技術路線兩方面進行闡述。

1）研究思路

本部分的測試是在長度為3.8km的區域進行，包括原始樣品點97個。首先對該研究區進行三維地質構造建模，然后進行網格剖分，每個網格屬性均用中心位置的屬性值代替，對所有網格應用序貫高斯模擬的方法進行插值得到整個區域的屬性分布情況。對測試數據進行變異函數分析擬合，得到理論變程a=1282m，并對該理論變程進行適當放大和縮小，這里以0.1的倍數逐次遞增和遞減，計算得到的插值結果數據與原始數據進行整體相對誤差分析，以誤差最小原則來獲取最優變程。不同變程對應的相對誤差如圖1。

圖1屬性相對誤差圖

圖1 中，橫坐標代表屬性建模時應用的變程大小，其中a表示理論變程，0.1a、0.2a、…、2a表示對變程適當的進行縮小和放大；縱坐標表示針對相應變程的插值結果同原始數據的整體相對誤差。分析圖中結果，可以得出：當變程為理論變程的0.3倍時為最低點，即插值結果與原始數據的相對誤差最小；而當變程大于理論變程時會呈現誤差不斷增大的趨勢，并且逐漸趨向平緩。由此可見，有效變程并不完全等于擬合計算出來的理論變程，而是存在一定偏差。因此研究怎樣通過有限的已知變量來獲取到更為合理的有效變程，對于構建精確的三維地質屬性模型至關重要。

2）技術路線

由上節可知，有效變程與屬性變化范圍和地質體規模兩個參數相關，且屬性變化與地質體規模不具有相關性、即為相對獨立參數，所以可以對兩個參數進行線性組合。有效變程與屬性變化范圍和地質體規模的關系式求解方法為：首先，保持地質體規模不變，對不同類型的屬性分別求解有效變程a，得到a與屬性變化的相關關系；其次，保持同一種屬性不變，對不同大小的研究區域分別求解有效變程a，得到a與地質體規模的相關關系；再次，將屬性變化和地質體規模同時進行變化，求解方程組，得到所有權重系數。流程圖如圖2所示，具體步驟如下：

圖2 理論模型求取流程圖

（1）保持地質體規模不變，對8組不同的屬性數據分別進行三維地質屬性建模，以整體相對誤差最小為標準來獲取每種屬性數據對應的有效變程。方法同上部分研究思路，每組數據均通過調整理論變程大小構建20個模型來尋找有效變程，共計構建屬性模型160個。建模結果如圖3（部分）。

這八組實驗的原始數據屬性變化范圍與其對應的有效變程如表1所示，其中有效變程代表建模結果最逼近于原始數據的變程。

對這八組數據進行散點圖的繪制，觀察地質體大小不變時，有效變程同屬性變化范圍之間的關系，如圖4所示。并對其進行分析擬合，得到兩者之間的數學關系。

由圖4可以看出，有效變程隨著屬性變化范圍的增大而呈現逐漸減小的趨勢，并且減小的幅度越來越小，形狀接近于拋物線形式，即有效變程同參數屬性變化大致具有二次關系，經過分析擬合處理，得到有效變程與屬性變化的相關關系式為

圖4 有效變程同屬性變化相關關系

（2）對于同一種屬性，在八個不同的研究區域進行模擬，即保持屬性不變，改變地質體的規模。方法同上，每組建立模型20個、共計構建屬性模型160個。這八個區域的三維地質屬性建模效果圖如圖5（部分）所示。

圖5 不同區塊三維地質屬性建模效果圖

測試數據選取了地質體規模不等的八個研究區域，對每個區域進行多次插值，對比選取使得插值結果相對原始數據整體相對誤差最小的合理變程，即有效變程，不同的地質體規模與有效變程的對應如表2。

表2 地質體規模與有效變程對應表

對這八組數據進行散點圖的繪制，觀察同一種屬性、即屬性變化不變的情況下，有效變程同地質體規模之間的關系，如圖6所示。

圖6 有效變程與地質體規模相關關系

由圖6可以看出，有效變程隨地質體規模的變化呈逐漸增大的趨勢，但是達到一定程度后，又逐漸減小。整體近似于正弦曲線的分布。經過分析處理以及正弦函數擬合，得到最優變程同地質體規模的相關關系為

（3）根據上述兩次實驗，可以總結出有效變程和屬性變化呈二次關系、和地質體規模呈近似正弦關系。又由于這兩個參數共同作用的效果影響最終的有效變程獲取，且各部分不相關。所以將上面的兩個方程進行線性組合，假定理論模型為

公式中有四個未知系數需要求取，為了估算每個參數對應的權重系數，需要求解四元一次方程組。本文通過選取不同屬性值、不同研究區域的四組測試數據進行研究，為求取每組參數對應的有效變程，每組構建20個屬性模型，來擇優選取有效變程。共計構建屬性模型80個。四組數據的屬性變化、地質體規模和有效變程對應如表3。

表3 屬性變化、地質體規模與有效變程對應表

將前四組數據分別代入公式中，通過矩陣計算的方式求解線性方程組得到每個變量的權重系數分別為：k1=0.0057，k2=-4.99，k3=-85，k4= 1415.64。綜上所述，有效變程同屬性變化和地質體規模兩個影響因素的最終理論模型為

3.2.2 理論模型驗證

為了驗證上節得到的理論模型、即式（6）的正確性，需要對其進行結果的驗證。驗證過程為：選取兩個不同的研究區域，并且針對兩種不同的屬性進行模擬，首先對其進行實驗人工調整來獲取有效變程，然后根據本文提出的理論模型進行計算得到有效變程，對比這兩種方法得到的有效變程。來證明模型的準確性。

本節利用表4中的2組數據進行結果驗證，分別將兩組數據的屬性變化和地質體規模帶入到式（6）中，得到理論模型求取出的有效變程分別為521.1m和404.1m。而多次實驗得到的結果分別為513.6m和378.7m。兩組數據對應的柱狀圖如圖7。可以看出計算出的結果同實驗結果非常接近，證明該公式是正確可行的。

表4 驗證數據

驗證結果證明有效變程理論模型是可行的，并且模型的提出克服了每次建模需要不斷調整參數來獲取最優插值效果的問題，從而大大提高了建模的效率，也為快速地構建合理的屬性模型提供一定的理論指導。

圖7 公式求取結果驗證圖

4 實例應用

為確保結果的實用性，本文選取了安塞區王窯油田的一個區塊作為實驗數據予以應用。該區構造變化較為簡單，構造背景為平緩的西傾單斜［11］。地層傾角為0.5°左右，埋藏深度約為1005～1060m。地質體平均孔隙度為13.7%，平均滲透率為2.27×10-3μm2，油藏平均厚度約為12m，屬于低孔低滲低產油藏［12～13］。進行測試的屬性數據為該區域的孔隙度參數。建模過程分為數據統計分析及預處理、擬合變異函數、鄰域搜索和插值計算四部分［1，14］進行。

1）數據統計分析及預處理

在進行插值前，首先需要對原始數據進行統計分析及轉換，提供平均數、最大值、最小值、峰度、偏態等多種統計量，并可進行數據轉換；此外，還要對原始數據進行預處理，如異常值的剔除或樣品屬性值過于密集時進行適當抽稀等。計算得到該屬性的屬性變化范圍約為0.08。

2）擬合變異函數

應用原始數據進行實驗變異函數計算，并選取球狀模型應用線性規劃法進行變異函數的擬合。由于該區域地質體長度為6.4km，根據3.2.1節所述參數選取原則，設置滯后距為128m。

擬合出理論實驗變異函數，得到理論變程為1177m，然后進行有效變程的校正，將屬性變化0.08和地質體規模6.4km分別帶入上節得到的公式（6）中，得到有效變程為1355m。

3）鄰域搜索

對于每一個待插點，給定搜索范圍，即鄰域半徑，并設置最大搜索點數，最大搜索點數決定最終求解的方程組個數。通過鄰域搜索獲取未知區域周圍的原始數據，剔除與待插點距離較遠相關性較弱的數據，本文選取計算出的有效變程1355m作為鄰域搜索的搜索半徑。對每個待估點進行鄰域搜索，得到與未知點相關的樣品值，用于后續待估值點的插值計算。

4）插值計算

選取合適的插值方法，如確定性建模的克里金插值法或者隨機建模的序貫高斯模擬方法，得到待估點周圍相關樣品數據的權重值，并進行加權求和得到所有網格的屬性值，最終形成連續、完整的全區三維地質屬性模型。本文采取序貫高斯模擬方法對每個待插點進行隨機模擬，序貫高斯模擬［15～18］是一種基于高斯場的序貫模擬方法，多用于符合高斯分布的連續型變量（如孔隙度、滲透率、飽和度等），能夠給予不確定性評價，并且相較克里金插值［19～20］更能真實地反映地質情況，因此應用非常廣泛。插值完成最終得到的屬性建模效果圖如圖8所示。

圖8 王窯區屬性建模效果圖

5 結語

1）通過實驗數據回歸分析、擬合的方法得到有效變程同屬性變化以及地質體規模兩個參數之間的理論數學模型，從而更加快速地獲得有效變程，為更加高效地構建合理、精確的三維地質屬性模型提供理論指導。

2）將該方法應用到油田的實際工作中，獲得更為合理的建模結果。能更好地為油藏數值模擬提供數據準備，也為致密砂巖地質體的評價與預測、水平井井段優化設計等方面提供可靠的依據，有一定的應用價值。

3）研究得到的有效變程理論模型針對所有的地質統計學方法均適用，因此可以應用于礦產資源預測和評價、礦山地質與礦產經濟、石油地質與煤田地質勘探與開發、水文工程地質等多個領域，具有很強的普適性。

［1］吳耕宇.基于規則網格的三維地質屬性建模關鍵技術研究［D］.北京：北京大學，2015.

WU Gengyu.Research on key technology of 3D geological attribute modeling based on regular grid［D］.Beijing：Peking University，2015.

［2］Cressie N.Fitting variogram models by weighted least squares［J］.Jounal Of The International Association For Mathematical Geology，1985，17（5）：563-586.

［3］黃滄鈿.變異函數的自動擬合方法［J］.新疆石油地質，2001，22（2）：155-157.

HUANG Cangdian.Automatic fitting method of variation function［J］.Xinjiang Petroleum Geology，2001，22（2）：155-157.

［4］孫洪泉.地質統計學及其應用［M］.徐州：中國礦業大學出版社，1990：20-45.

SUN Hongquan.Geological statistics and its application［M］.Xuzhou：China Mining University Press，1990：20-45.

［5］胡小榮，俞茂宏.理論變異函數球狀模型的加權線性規劃法擬合［J］.地質與勘探，2001（5）：45-48.

HU Xiaorong，YU Maohong.Weighted linear programming method for the spherical model of the theoretical variation function［J］.Geology And Exploration，2001（5）：45-48.

［6］程勖，楊毅恒.變差函數中加權系數擬合方法的改善［J］.地球物理學進展，2009，24（1）：362-366.

CHENG Xu，YANG Yiheng.Improvement of fitting method of weighted coefficient in variation function［J］.Progress In Geophysics，2009，24（1）：362-366.

［7］滕召良，呂文生，楊鵬，等.基于加權線性規劃法的變異函數球狀模型自動擬合研究［J］.西部探礦工程，2014，26（9）：73-76.

TENG Zhaoliang，LV Wensheng，YANG Peng，et al.Automatic fitting of the spherical model of the variation function based on the weighted linear programming method［J］.West-China Exploration Engineering，2014，26（9）：73-76.

［8］張仁鐸.空間變異理論及應用［M］.北京：科學出版社，2005：13-18.

ZHANG Renduo.Spatial variability theory and its application［M］.Beijing：Science Press，2005：13-18.

［9］PYRCZ M J，DEUTSCH C V.Geostatistical reservoir modeling［M］.New York：Oxford University Press，2002：40-70.

［10］DEUTSCH C V，JOURNEL A G.GSLIB Geostatistical software library and user's guide［M］.New York：Oxford University Press，Second edition，1998：11-14.

［11］田豐.安塞油田王窯區長6地質體微觀特征及氣驅提高采收率研究［D］.西安：西安石油大學，2013.

TIAN Feng.EOR research on microscopic characteristics and 6 gas geological body of Ansai oilfield Wangyao district［D］.Xi'an：Xi'an Petroleum University，2013.

［12］趙向原，曾聯波，靳寶光，等.裂縫性低滲透砂巖油藏合理注水壓力——以鄂爾多斯盆地安塞油田王窯區為例［J］.石油與天然氣地質，2015（5）：855-861.

ZHAO Xiangyuan，ZENG Lianbo，JIN Baoguang，et al. Rational water injection pressure in fractured low permeability sandstone reservoir——A case study of Ansai oil field in Ordos Basin［J］.Petroleum And Natural Gas Geology，2015（5）：855-861.

［13］許永濤，朱玉雙，張洪軍，等.安塞油田王窯區長6儲層特征及孔滲特性控制因素［J］.石油地質與工程，2011，25（4）：25-28.

XU Yongtao，ZHU Yushuang，ZHANG Honghui，et al. The characteristics of controlling factors of Ansai oilfield infiltration Wangyao 6 reservoir characteristics and hole［J］.Petroleum Geology and Engineering，2011，25（4）：25-28.

［14］姚凌青.基于地質統計學的三維屬性建模方法研究［D］.北京：北京大學，2008.

YAO Lingqing.Research on the method of 3D attribute modeling based on geo statistics［D］.Beijing：Peking University，2008.

［15］趙彥鋒，孫志英，陳杰.Kriging插值和序貫高斯條件模擬算法的對比分析［J］.地球信息科學學報，2010，12（6）：767-776.

ZHAO Yanfeng，SUN Zhiying，CHEN Jie.Comparison and analysis of Kriging interpolation and sequential Gauss conditional simulation algorithm［J］.Journal of Earth Information Science，2010，12（6）：767-776.

［16］BOISVERT J B，DEUTSCH C V.Programs for kriging and sequential Gaussian simulation with locally varying anisotropy using non-Euclidean distances［J］.Computers And Geosciences，2011，37（4）：495-510.

［17］胡先莉.序貫條件模擬方法研究及應用［D］.成都：成都理工大學，2007.

HU Xianli.Research and application of sequential conditional simulation method［D］.Chengdu：Chengdu University of Technology，2007.

［18］姚興苗.各向異性序貫高斯隨機模擬研究［D］.成都：電子科技大學，2014.

YAO Xingmiao.Anisotropic sequential Gauss simulation［D］.Chengdu：University of Electronic Science and Technology，2014.

［19］張靖.基于克里金算法的點云數據插值研究［D］.西安：長安大學，2014.

ZHANG Jing.Study on Kriging algorithm of point cloud data based on interpolation［D］.Xi'an：Chang'an University，2014.

［20］王亭.顧及各向異性的三維Kriging空間插值方法研究［D］.南京：南京師范大學，2013.

WANG Ting.Research on 3D Kriging spatial interpolation method considering anisotropy［D］.Nanjing：Nanjing Normal University，2013.

Variable Range Optimization in Geostatistics

HAO HongruiPAN MaoLIU ZhenCAO Kai
（School of Earth and Space Sciences，Peking University，Beijing100871）

It is necessary to adopt sample data of various geological attributes for 3D geological attribute modeling，the interpolation of unknown regions is performed by variogram analysis and proper interpolation algorithm to obtain the distribution of attributes in the study area.In order to build more accurate attribute models，a certain amount of human-computer interaction is usually needed during the modeling process.But in the past property modeling，the effective selection of variable range is usually adjusted according to the result of attribute modeling，without a reference standard，and this has great influence on the efficiency and precision of attribute modeling.Based on the problems above，this paper attempts to study the optimization of geostatistical variable range，and seek the theoretical model between effective variable range and related influencing factors，so as to improve modeling efficiency and accuracy.

attribute modeling，variation function，effective variable range，attribute change，geological scale，sequential Gaussian

P628

10.3969/j.issn.1672-9722.2017.06.035

2016年12月3日，

2017年1月30日

國家自然科學基金項目“基于語義的多分辨率儲層數據組織與管理”（編號：41472113）資助。

郝紅瑞，女，碩士研究生，研究方向：三維地質屬性建模及插值算法。潘懋，男，博士，教授，研究方向：石油地質、環境地質和信息地質。劉振，男，博士研究生，研究方向：地應力建模及巖石力學參數建模。曹凱，男，博士研究生，研究方向：三維地質構造建模和相建模。