基于多維尺度變換和自適應增強回歸進行插值
——以三維地質屬性體為例

李楊， 鐘寶榮， 許曉宏*， 梁梓君

(1.長江大學地球科學學院, 武漢 430002； 2.長江大學計算機科學學院, 荊州 434100; 3. 中國石油新疆油田公司, 克拉瑪依 834000)

在油氣勘探開發過程中，三維地質建模能夠幫助地質人員了解地質發育特征，輔助油氣開采決策[1-2]?？臻g插值是三維地質建模的關鍵部分，提高插值的精度可以節省油氣勘探開發成本。

目前應用于地質建模的插值方法分為確定性方法和地統計方法。確定性方法通過數學函數，根據相似程度或平滑程度進行插值。地統計方法利用測量點的統計屬性進行插值。

在確定性插值方面，趙仕威等[3]利用反距離加權和多重二次曲面函數進行三維地層可視化，在稀疏樣本的條件下取得了良好的效果，但是并沒有考慮地質發育的各向異性特征。Hillier 等[4]以徑向基(radial basis function，RBF)為基礎，利用梯度約束的數據驅動特征分析各向異性，能夠應用于分散的多元地質結構。龐慶剛等[5]利用緊支撐徑向基(compactly supported radial basis functions,CSRBFs)進行三維地質建模，通過樣本的法向量來考慮各向異性，實現了在非均勻數據中快速插值。然而基于RBF的地質建模方法普遍過于平滑，在固體礦產方面有優勢，但是在油氣存儲區域不連續且封閉的地質狀況中適用性欠佳。另外，RBF中不同的基函數插值效果差異過大。黃牧等[6]利用Voronoi圖實現了快速自動化三維地質建模，提高了建模效率，然而Voronoi圖是以三角網為基礎，與常用的規則六面體相比，必然伴隨著數據結構復雜，兼容性差，轉換流程耗時等問題。

在地統計插值方面，插值方式主要以克里金(Kriging)為主。張小艷等[7]利用變異過程中修正變異方向的自適應差分進化算法來優化Kriging，進行煤層的三維可視化地質建模。楊陽等[8]通過利用半變異函數理論對樣本進行分析，采用普通克里金進行海洋地層厚度的地質建模。蔡文軍等[9]利用Kriging方法在含油的灰巖儲層中進行三維孔隙壓力建模。Li等[10]利用深度神經網絡模擬半變異函數進行Kriging插值，減少了理論Kriging方差。然而，Kriging方法具有一定的局限性：首先，Kriging方法應用的前提是樣本必須符合特定的隨機場(如普通克里金必須要正態分布)；其次，Kriging方法在插值過程需要經過繁瑣的人為分析流程，如半變異函數分析，塊金、基臺、變程等參數的確定。

地理學定律[11-13]認為地理空間分布具有空間相關性和異質性。從空間相關性可以得出地理分布與距離相關，從空間異質性可以得出地理分布具有各向異性特征。在地質統計學中通常利用傾斜角和方位角來計算各向異性特征。由此得出，影響插值的因素包括樣本分布的距離、傾斜角、方位角、樣本點屬性值。相對于其他降維方法，利用多維尺度變換[14-15](multidimemsional scaling，MDS)能夠最大限度地減少數據的損失的優勢，對影響插值的因素進行分析并降維，提取影響插值的特征，其結果作為訓練模型的特征集。相對于人工智能領域中其他的回歸算法，利用自適應增強回歸[16](adaptive boosting regressor，ABR)具有高精度、強學習器、防止過擬合等優勢，對已知樣本點的特征進行訓練，探索地理分布的空間相關性，建立插值模型。接著，把未知點與樣本點的空間相關性信息帶入到ABR模型中求得已知點的值。插值模型建立了多個ABR模型，每個ABR模型相對獨立，所以最后利用三維均值濾波(mean filtering，MF)進行平滑處理，以此來使結果符合地質屬性分布的連續性特征。

在地理學定律和地質統計學的前提下，針對上述文獻中的問題，以常見的規則六面體為基礎數據結構，提出利用人工智能算法(ABR算法)進行三維地質屬性體插值，提高了插值精度，實現自動化插值。本文的研究對于推進地質建模的智能化發展具有一定的意義。

1 技術路線

利用sklearn算法庫中的MDS和AdaBoostRegressor(ABR)進行插值。插值方法跟大多數插值方法類似，采用距離未知點最近的n個已知樣本點進行未知點的估算。三維空間插值方法概述流程如圖1所示，其中P0為未知點P1,P2,…,Pn為距離P0最近的n個已知樣本點，其順序以由近到遠排序。

圖1 利用MDS、ABR和MF插值流程概述Fig.1 Summary of MDS,ABR and MF interpolation process

2 MDS降維

空間插值中未知點的屬性值受其最近鄰的若干個已知樣本點的距離、已知點的屬性值、拓撲關系影響，然而，一些變量對插值結果影響很小。為了提出影響插值的主要因素，減少次要因素的影響，提高插值結果的穩定性，將這些因素作為MDS的降維對象，得到一維數據。

采用未知點(P0)最近鄰的n個樣本點(P1,P2,P3,…,Pn)來預測未知點的屬性值。假設每個已知樣本點為未知點，然后根據這些已知樣本點建立ABR模型。例如，假設P1為未知點，P2,P3,…,Pn為已知點(共n-1個)。把這些n-1個已知點根據到P1的距離由近到遠重新進行排序，排序后的點為Psort1,Psort2,…,Psort(n-1)。對n個已知樣本點(P1,P2,P3,…,Pn)依次假設每個樣本點為未知點。每次假設就能得到一個訓練樣本。共產生n個訓練樣本(每個訓練樣本的特征有n-1個)，把這些訓練樣本帶入到ABR模型中。當假設P1為未知點時，利用MDS降維的描述如圖2所示。

(1)Z表示屬性值。Zsort1表示Psort1的屬性值，Psort1為距離P1最近的點；Zsort2表示Psort2的屬性值，Psort2為距離P1第2近的點；以此類推，Zsort(n-1)表示Psort(n-1)的屬性值，Psort(n-1)為距離P1第n-1近的點(距離P1最遠的點)。

(2)h表示兩個點之間的距離。h1(sort1)表示P1～Psort1的距離；h1(sort2)表示P1～Psort2的距離；以此類推，h1[sort(n-1)]表示P1～Psort(n-1)的距離。

(3) 由于地理現象中的屬性特征在不同的方向上呈現出不同的變化速度，所以本文用方位角ɑ表示各向異性特征。α1(sort1)表示P1與Psort1的方位角；α1(sort2)表示P1與Psort2的方位角；以此類推，α1[sort(n-1)]表示P1與Psort(n-1)的方位角。

(4)X1(P1Psort1)表示P1與Psort1的h、Z、α進行降維后的數據。

(5) ABR的訓練集中包含n個訓練樣本，每個訓練樣本包含n-1個特征。特征根據已知樣本點與未知樣本點的距離遠近獲取。特征1表示距離假設的未知點P1最近的樣本點Psort1對假設未知點P1的影響，用X1(P1Psort1)表示。特征2表示距離假設的未知點P1第二近的樣本點Psort2對假設未知點P1的影響，用X1(P1Psort2)表示。以此類推，共產生n-1個特征。

在圖2中需要注意的是，在對所有樣本點假設為未知點時，每次假設，在計算特征1時都會產生一個Psort1，這些Psort1并不是代表同一個樣本點，而是對假設未知點后，其余的已知點到未知點的距離排序由近到遠而產生的。例如，假設P1為未知點時，Psort1表示其余已知點中距離P1最近的點。假設P2為未知點時，Psort1表示其余已知點中距離P2最近的點。

圖2 n個樣本點MDS降維處理Fig.2 MDS dimension reduction for n sample points

3 ABR空間插值

AdaBoost(adaptive boosting)[16]算法是在多個弱學習模型的基礎上建立起來的，把樣本帶入弱學習模型中進行訓練及驗證，根據誤差對樣本權重進行更新，然后帶入到新的弱學習模型中訓練，最后對這些弱學習模型的結果進行加權求和(或者投票)的方式組合，求得模型的結果。AdaBoost在每次訓練弱學習模型時，都對樣本進行了驗證以及誤差計算，所以樣本進行了多次驗證及誤差分析。常用的分類回歸算法一般都需要用到訓練集及驗證集，以此來更新模型權重獲取最優模型。AdaBoost算法在計算過程中已經把訓練樣本當作驗證集多次進行驗證來獲取最優模型，所以額外的驗證集并不是必須的。

AdaBoost可以運用于分類(adaptive boosting classifier)和回歸(adaptive boosting regressor)，本文中主要運用其回歸特性，具體推導過程及應用見文獻[16]。

在利用ABR建立模型時，由于樣本量較少，很容易產生較好的訓練模型，選用的弱學習模型為CART決策樹[17]，損失模型為線性模型，最大迭代次數為50，學習率為0.1。

3.1 ABR模型訓練

采用ABR方法建立插值模型，一個未知點建立一個ABR模型。若采用n鄰域樣本點進行插值則訓練ABR模型所需的數據集個數為n，訓練集的輸入中每個樣本的特征為n-1個。ABR模型訓練時需要確定訓練集的輸入和輸出(圖3)。

(1) ABR訓練集輸入：ABR訓練集的輸入為經過MDS降維到一維的數據集(圖2中降維后的特征，也就是圖3中的輸入)。

(2) ABR訓練集輸出：ABR訓練集的輸出為圖2中假設的未知點的真實值(圖3中的輸出)。

圖3 訓練ABR插值模型Fig.3 Training ABR interpolation model

所有的訓練集的特征之間可能有一定的關聯，因為它們都是從P1,P2,…,Pn獲取的。然而每個訓練樣本所表達的意義以及獲取的計算方式都是獨立的(圖2中的特征)。因此，盡管它們具有一定的關聯，但是仍然可以用于建立ABR模型。

3.2 未知點預測

通過上述訓練得到ABR模型，然后獲取未知點P0與其最近的n-1個點之間的相對距離h，方位角α，樣本點的屬性值Z。P1,P2,…,Pn-1是距離P0最近的n-1個樣本點，并且這些樣本點是根據到P0的距離，由近到遠排序的。如圖4所示。

圖4 未知點值的預測Fig.4 Prediction of unknown point

(1)Z1表示已知樣本點P1的真實值。

(2)h01為未知點P0與已知樣本點P1的距離。

(3)α01為未知點P0與已知樣本點P1的方位角,用來表示各向異性特征。

(4)X0(P0P1)為未知點P0與已知樣本點P1的Z1、h01、α01進行MDS降維后的結果，同理得到X0(P0P2),…,X0(P0Pn-1)。

(5) 把X0(P0P1),X0(P0P2),…,X0(P0Pn-1)代入ABR模型中預測出待估計點P0的屬性值Z0。

(6) 特征1, 特征2,…, 特征n-1表示已知樣本點距未知點P0由近到遠的影響特征。

4 三維MF平滑處理

均值濾波(MF)是一種鄰域運算，鄰域運算是通過對某一像素的鄰域進行加權運算求得新的值代替該像素，加權運算所用到的矩陣叫作局部算子。本文的三維MF采用3×3×3的局部算子，具體計算流程如下。

(1) 設定局部算子為3×3×3，同時在X、Y、Z方向分別擴展一行(列)，新增行的像素值賦值為已知的最近像素值。

(2) 以待處理像素為中心，選取3×3×3的矩陣。

(3) 計算M矩陣中3×3×3個矩陣的均值，并替換掉待處理像素。

(4) 依次遍歷其他像素，重復步驟(2)和步驟(3)。

圖5中是ABR進行插值的結果(MF處理前)與MF處理后效果對比，可以看出經過MF處理后可以減少多個ABR模型計算的插值結果產生的不連續性，使插值結果符合地質屬性分布具有連續性的特征。

圖5 均值濾波前后對比Fig.5 Comparison result before and after MF

5 插值結果分析

孔隙度是指巖石孔隙體積與巖石表面體積之比，是評價油氣分布的重要指標。以遼河油田某區塊孔隙度數據為例，進行了三維空間插值。由于設備或其他原因，某些孔隙度數據可能無效。圖6顯示了11口井的分布，其中驗證樣品(3口井，共76個樣品)用紅色矩形標記。其他9口井用于插值(總共203個樣本)。由于鉆井位置與含油量、孔隙度、滲透率、飽和度等因素有關，無法保證井的均勻分布。在圖6中，左側的孔隙度較大，這有利于儲油、高含油量和更多的井。其他地區孔隙度小，含油量低，井數較少。插值結果的步長在X方向為5 m，Y方向為5 m，Z方向為1 m(每個最小單位網格的插值結果為5×5×1=25 m3)，總共有27 000個網格。工區面積為150×150 m2，地層厚度為30 m(海拔為-1 640～-1 670 m)。

圖6 測試樣本分布Fig.6 Test sample points distribution

由于反距離加權插值沒有考慮地質發育的各向異性特征；Voronoi圖插值結果數據結構復雜，可重用性較弱；RBF插值結果過于平滑，只適用于固態礦產的地質體建模，而地質統計學中的克里金方法在含油氣層的地質體建模中最為常用。所以將本文方法與地統計學插值方法中的普通克里金(ordinary Kriging，OK)的插值結果進行比較，分析插值效果。OK使用Python中的pykrige包。

當存在多個樣本點時，并非所有的樣本點都將用于插值。通常選取未知點附近的n個已知樣本點建立插值模型，用于計算未知點的值。在本文中，搜索距離未知孔隙度點(圖1～圖4中的P0)最近的n個已知孔隙度樣本點(圖1～圖4中的P1,P2,…,Pn)。圖7顯示了MDS-ABR-MF(由MDS、ABR和MF處理的插值結果)和OK插值結果之間的比較，其中n=5、n=10和n=15的插值效果。

圖7 不同方法插值效果對比Fig.7 Comparison of interpolation results by different methods

用3口測試井中的76個樣本點進行誤差分析。通過均方根誤差(root mean squared error, RMSE)和平均絕對誤差(mean absolute error, MAE)來測試本文方法的插值精度。由表1可以看出MDS-ABR-MF運用本文中的測試數據進行插值比OK插值結果在精度有所提升。

表1 MDS-ABR-MF和OK插值結果與3口測試井之間的RMSE和MAETable 1 The RMSE and MAE between the MDS-ABR-MF and OK interpolation results and 3 test wells

OK與MDS-ABR-MF插值結果的誤差產生的原因可能如下。

(1) OK方法假設在標準正態分布的情況下進行插值。然而，無論原始數據如何處理，這些數據都只是近似正態分布，很難達到標準正態分布。

(2) 在各向異性分析和半變異函數分析中，OK插值結果受參數選擇的影響。例如，孔隙度范圍為0～1，通過測試OK方法，不合理的人為調整可能會產生負值，因為OK插值中計算的權重可能是負值。然而，在ABR插值中，由于弱學習模型為決策樹，算法的基本原理決定了未知點的值是已知點的加權計算，所以當ABR算法中決策樹為弱學習模型時，決定了權重不能為負，所以最終結果不可能為負值。

(3)圖7中，OK插值中的不連續性是因為選取的n個已知樣本點的孔隙度值可能會相同。此時，在OK算法中，計算的協方差矩陣無法獲得逆矩陣，因此本文取平均孔隙度值，這也是OK算法的一個弊端。

6 結論

利用MDS、ABR和MF進行探索三維空間樣本點之間的相關性來實現三維地質屬性體插值。擺脫了反距離加權方法中不考慮各向異性特征的弊端。通過對比傳統的普通克里金(OK)插值，得出以下結論。

(1) OK插值需要特定的隨機場(正態分布)，在非正態分布區域執行OK插值時，插值結果精度可能不高。MDS-ABR-MF方法不受特定隨機場的約束。

(2) OK插值需要人工半變異函數分析、各向異性分析、塊金、基臺和變程等繁瑣流程分析，不適當的參數調整也可能產生異常。MDS-ABR-MF插值可以縮短繁瑣的參數調整過程，在保證精度的同時實現自動插值。