半變異函數及取樣間距對克里金法在海洋地層分析中的影響研究*

楊 陽 李 颯 何福耀 施展玲

(①天津大學建筑工程學院 天津 300350)(②中海油上海分公司 上海 200032)

0 引 言

在海洋工程建設前，對施工場地進行工程地質勘察，了解海底地形、地貌以及海底地層分布特征尤為重要。工程鉆探是一種常用的勘察技術手段，取樣獲得的鉆孔數據可以很好地反映地層垂向分布特征(Lemon et al., 2003)。陸地施工方便，獲得鉆孔數據較容易且成本低。而對于海洋工程，由于勘察條件的限制及海底地理環境的復雜，獲取鉆孔數據極其困難，且成本較高。相比陸地，海上鉆孔取樣密度大幅度降低，僅依靠少量的鉆孔數據顯然不能完整地獲取地層厚度分布特征(李一保等， 2007； 李安龍等， 2016)。淺剖數據在測線方向上具有良好的連續性，若將極有限的鉆孔數據和連續的淺剖數據有機結合起來，將有利于提高地層厚度判別精度(周曉， 2014)。

地層厚度在空間分布上具有一定的隨機性和結構性。因此，可將地層厚度視為區域化變量，采用克里金法進行空間插值估計地層厚度(Cressie， 1991； Chilès et al.，1999； 周小文等， 2005； 許國等， 2013； 劉曉磊等， 2018)。在利用克里金法估計地層厚度方面，國內外學者展開了廣泛的研究。李曉軍等(2014)采用IK-OK方法較準確地估計了崇明越江隧道工程段地層厚度； 杜文鳳等(2010)將地層厚度和地震振幅作為源數據，利用協同克里金法預測煤層厚度； Hengl et al. (2004)則采樣回歸克里金法預測西伯利亞某土層厚度； Jang et al. (2016)采用指示克里金多閾值組合的方法在估計沖積扇地區低滲性土材料厚度中取得良好的效果。但大多數學者在利用克里金法估計地層厚度方面，研究區域多局限在陸地，針對海底地層厚度估計的研究較少。正確且準確地識別海底地層厚度分布對海洋工程非常重要。因此，將克里金法應用在海底地層厚度識別中具有重要的意義。

對地層厚度進行克里金插值前，首先需要對已有數據進行變異結構分析，通過變異結構分析構建一個半變異函數模型； 再利用該模型計算未知區域地層厚度，從而得到地層厚度分布特征(張征等， 1995； 崔洪慶等， 2003)。目前國內外針對地層厚度插值的研究多在分析比較不同克里金法插值的優劣性(Penizek et al.，2006； Kuriakose et al.，2009； 易湘生等， 2012； 張琳娜等， 2016)，如：Penizek et al. (2006)采用了不同的克里金法對土壤厚度進行插值，發現協同克里金優于普通克里金、線性克里金和回歸克里金； Kuriakose et al. (2009)則認為回歸克里金優于普通克里金。對半變異函數模型選取的研究較少，而半變異函數模型選取適當與否將直接影響插值精度。因此，有必要對最優半變異函數模型作進一步探討(Kambhammettu et al.，2011)。同時，對于半變異函數模型最優性的判別，通常采取決定系數越大越佳，殘差越小越佳的原則確定(劉愛利等， 2012)。但是，很多時候會出現各模型的決定系數和殘差的差異很小而無法擇優(施小清等， 2009； 張濤等， 2016)，甚至出現決定系數最大的模型，其插值結果卻不是誤差均值最小的情況(巫振富等， 2012)。若僅依據上述原則選擇模型可能無法保證所選模型為最優，造成預測結果的偏差，所以對于最優的模型需進一步驗證其可靠性。

此外，插值精度不僅受半變異函數模型的影響，也受樣本數量、取樣間距、原始數據等的影響(朱會義等， 2004)。通過目前的研究資料發現，對插值精度的研究主要涉及取樣空間尺度(Robeson， 1993)、樣本容量(Journel et al., 1978； Bourennane et al.， 2000)等的影響。因此，了解取樣間距對插值精度的影響尤為重要。

為了研究適合海洋地層厚度的最佳半變異函數模型，本文基于GS+8.0軟件，以某井場3個反射層組的淺剖數據為源數據，鉆孔數據為實測驗證數據，采用克里金法，利用不同的半變異函數模型分別對地層厚度變異函數進行擬合對比分析。一方面探討更適合該井場地層厚度插值的模型，同時研究取樣間距對預測誤差的影響，為該井場地層厚度估計提供依據； 另一方面給出海底地層厚度估計中半變異函數模型選擇的一般思路，為其他海域地層厚度插值提供參考。

1 地統計學方法

克里金插值法屬于地質統計法的一種，它是以無偏為約束，尋求估計方差最小的一種估計方法(孫洪泉， 1990； 王鵬等， 2013)。根據地質統計理論，半變異函數既能考慮空間統計特征，又能同時描述區域化變量在空間分布的隨機性和結構性，完整反映區域化變量的空間變異規律。由于地層厚度符合區域化變量的特征，因此可以用地層厚度取樣值建立半變異理論函數，利用半變異理論函數求出研究區域內任意一點處的地層厚度。

1.1 變異結構分析

半變異函數是計算研究區域內地層厚度變化的關鍵。關于地層厚度的理論半變異函數定義為地層厚度Z(x)在點x和x+h處數值之差的方差之半，其表達式為：

(1)

在實際計算γ(x，h)時，采用地層厚度實驗半變異函數計算。測線方向上相距h的xi和xi+h的取樣點地層厚度值為Z(xi)和Z(xi+h)，測線方向上相距h的取樣點的數據對數目為N。對于任意的h，實驗半變異函數γ*(x，h)可以表示為：

(2)

假定任意一點x0處的地層厚度待估值Z*(x0)是周圍實測值Z(xi)的線性組合，即：

(3)

式中，λi為權重系數；Z*(x0)為點x0處的地層待估值；Z(xi)為已知點i=1，2，3，…，n處的地層厚度實測值。根據無偏約束條件可推導出：

(4)

根據最小方差條件，可得：

(5)

式中，γ(xi，xj)為地層厚度采樣點xi和xj的理論變異函數。

(6)

1.2 半變異函數理論模型

半變異函數是半變異值γ(x，h)與距離之間的函數。其中半變異值表示空間相關性，等于所有取樣點中距離為h的厚度值之差的期望。實驗半變異函數通常由3個參數表征：(1)塊金值c0，主要是由于空間變異和測量誤差導致； (2)基臺值c0+c，一般來說表征數據對間的最大方差； (3)變程a，代表土性參數存在相關性的最遠距離。實驗變異函數是根據采集的數據獲得的非連續分布函數，需采用理論變異函數對實驗變異函數進行擬合以獲得連續的變異性評價，巖土工程中常用的半變異函數的理論模型見表 1。

表 1 巖土工程常用變異函數模型Table 1 Typical theoretical semi-variogram models in geotechnical engineering

2 地層的厚度

地層厚度是識別地層結構和構造的重要依據之一。對于海洋地質勘查，通常通過淺地層剖面儀所觀測到的海底地層聲波接收反射特征得到地層分布特征，反映地層厚度變化規律。

淺地層剖面儀是一種基于水聲學的走航式水下淺部地層結構和構造的地球物理探測方法，其工作原理是利用聲波在不同介質中傳播速度的差異性。淺地層剖面儀的發射基陣(震源)按一定時間間隔向下發射聲波，聲波穿透海水抵達海底后，一部分聲波在分界面處反射并返回，另一部分則被接受基陣(水聽器)接收。

產生反射的條件是分界面兩側存在波阻抗差(反射系數R)。在淺地層剖面調查中，通常認為聲波是垂直入射的，此時反射系數R有：

(7)

圖 1 井場過井剖面Fig. 1 Shallow profile section of over well

式中，ρ1V1、ρ2V2分別是上下兩層介質的聲阻力率，ρ和V分別表示介質的密度和聲速。當相鄰地層介質存在較大的聲阻力差，在界面就會有較強的聲強，接收到較強的反射信號，在淺地層剖面儀終端顯示器上反映灰度較強的分界線。通過觀測分析海底地層的聲波反射時間可計算出地層厚度。

3 地層數據分析

為盡可能滿足采樣點數據對的要求，選定采樣點的距離間隔為50im，每組地層厚度數據點數目為78個，符合數據配對數目應大于30的原則(Robeson， 1993)。

3.1 半變異函數擬合

GS+是目前專業的常用地統計學軟件之一，在對區域化變量的變異結構分析方面具有較明顯的優勢。因此，本文利用GS+8.0軟件，選取高斯函數、球狀函數和指數函數進行地層厚度變異函數擬合。有效滯后距為1950im。擬合結果如圖 2所示，模型具體參數見表 2。

圖 2 擬合半變異函數圖Fig. 2 Fitted semivariograma. 第Ⅰ層； b. 第Ⅱ層； c. 第Ⅲ層

表中，r2為決定系數，RSS為殘差。塊金值c0與基臺值c0+c的比值稱為空間相關度，可以反映影響因素中結構性因素和隨機性因素的作用。根據Cambardella et al.(1994)的劃分原則，空間相關度小于25%，表明變量空間相關性強； 大于75%，表明空間相關性較弱； 處于兩者中間，說明空間相關性中等。從3種變異函數模型參數可以看出，3組地層厚度在3種模型下的塊金值c0與基臺值c0+c的比值均小于10%，說明研究區域地層厚度分布以結構性變異為主，具有強烈的自相關性。

表 2 變異函數參數統計Table 2 Variance function parameter statistics

海洋地層在沉積過程中大致受巖土體形成及后期改造過程中各種宏觀規律控制，使得不同點的地層厚度之間存在極強的空間結構性。當c0和c0+c發生變化時，變程a的變化幅度較大，其中變程大小依次為指數模型、球狀模型和高斯模型，但均未超過研究區域范圍。通過決定系數越大，殘差越小的原則確定模型最優性(劉愛利等， 2012)。對比模型擬合參數，同組地層各模型的決定系數較為接近，高斯模型略高于球狀模型和指數模型。

表 3 交叉驗證結果分析Table 3 Analysis of cross validation results

對比殘差，指數模型最大，球狀模型和高斯模型次之。僅從擬合效果角度分析，3種函數模型中最優的是高斯模型，其次是球狀模型。

3.2 交叉驗證

在對3組地層厚度進行克里金插值前，應對實測采樣數據模型進行整體的精度評價。交叉驗證是一種常用的精度驗證方法，即每次移除一個或多個實測采樣數據，然后使用其他位置的數據來預測與其相關的數據。重復操作，最后將每個數據點的真實值和估計值進行統計回歸分析，進而評價模型的精度。表 3為交叉驗證誤差統計結果。

圖 3 第Ⅰ層插值結果Fig. 3 Ordinary Kriging of layer Ⅰa. 指數模型； b. 球狀模型； c. 高斯模型

圖 4 第Ⅱ層插值結果Fig. 4 Ordinary Kriging of layer Ⅱa. 指數模型； b. 球狀模型； c. 高斯模型

表3中，SE是回歸系數的標準誤差，其值越小說明精度越高；r2是相關系數的平方，其值越大說明回歸效果越好；SE預測值是標準誤差預測值，其值越大說明誤差越大。通過這3個指標可以綜合評價不同半變異理論模型的預測精度。從表3中可以看到，在第Ⅰ、Ⅱ層中，3種模型的r2接近，較難看出3種模型回歸效果的優劣，但從SE預測值中可以看出球狀模型的誤差明顯小于其余兩種模型。在模型數據交叉驗證環節，球狀模型表現出的整體精度最高，高斯模型次之，指數函數的誤差最大。可見，結構分析最佳的模型不一定是誤差最小的模型。在出現結構分析和誤差分析的擇優結果不一致的情況下，很難說明球狀模型和高斯模型哪一個是絕對最優。因此，需對插值效果進一步對比分析以選擇最適合的模型。

圖 5 第Ⅲ層插值結果Fig. 5 Ordinary Kriging of layer Ⅲa. 指數模型； b. 球狀模型； c. 高斯模型

3.3 地層厚度插值效果分析

在對數據完成結構分析的基礎上，經交叉驗證的方法對不同理論模型精度評定后，采用普通克里金法對3組地層厚度(沿測線方向)進行插值。插值結果(圖 3～圖 5)。對比3組地層厚度插值云圖，可以看到，球狀模型和高斯模型在第Ⅰ反射層組0、690im、1300im和1600im處，第Ⅱ反射層組910im、1560im、1950im和2340im處，第Ⅲ反射層組300im、520im、1410im和3380im處表現出較高的相似性，與結構分析和交叉驗證分析結果相吻合。

表 4 K1插值結果對比Table 4 The actual thickness of K1 and prediction

利用研究區域內的淺鉆孔K1(非采樣點)的實測厚度數據和利用普通克里金法插值得出的地層厚度進一步檢驗(表 4)。從圖中可發現，不同理論模型插值的結果在淺鉆實測點均存在不同程度的誤差。其中，球狀模型的最小，在3組地層厚度誤差分別是1.9%、4.4%和2.3%。在地層厚度變化幅度較大的第Ⅰ、Ⅱ層中尤為明顯。由以上分析可知，球狀模型在該井場地層厚度離散數據插值中具有較高的實用性。

從K1插值結果數據對比可以看到，結構分析最佳的模型，即決定系數最大，殘差最小的模型，不一定是誤差最小的模型。因此，在進行半變異函數的選擇時，對于其可靠性驗證需要通過以下步驟：首先進行結構分析，其次通過交叉驗證評定模型整體精度，最后進一步驗證克里金插值效果。通過上述方法綜合分析來選擇最為適宜的半變異函數模型。僅根據結構分析，就有可能造成判斷結果的偏差。圖 6為最優半變異函數模型選擇流程圖。

圖 6 最優半變異函數模型選擇流程圖Fig. 6 Flow diagram of selection of optimal Semi-variogram function

4 取樣間距的影響

圖 7 第Ⅲ層交叉驗證Fig. 7 Cross validation of layer Ⅲa. 50 m； b. 75im； c. 100im； d. 125im； e. 150im

求解半變異函數值時，若取樣間距增大，取樣點空間分布離散稀疏，無法反應小范圍內的空間變異性，即無法準確反映地層厚度變異規律(Cambardella et al.， 1994)。本文通過設置不同的取樣間距來研究其對插值結果的影響。為滿足數據對數目大于30的要求，故將取樣間距最大值控制在150im以內，分別為50im、75im、100im、125im、150im。半變異理論模型選取上文中得到的最優模型球狀模型。將不同取樣間距獲得的地層厚度數據進行交叉驗證。圖 7分別為第Ⅲ組地層不同取樣間距時計算出的預測值和實際值的回歸曲線，其中橫坐標為預測值，縱坐標為實際值，虛線為理想回歸曲線。從3組地層厚度數據交叉驗證誤差統計結果(圖 8)可以看出，取樣間距越大，SE值越大，擬合精度越低。第Ⅰ、Ⅱ層的相關系數r2受采樣間距影響較小。第Ⅲ層的相關系數有一定程度的降低，當取樣間距為150im時，相關系數為0.801。通常認為r2>0.8，即可認為具有較好的回歸效果，因此該條件下地層厚度仍具有較強的自相關性。同時也說明對于球狀模型，增大取樣間距對地層厚度變化較大的地層回歸效果影響較小，而對于地層厚度變化較平緩的地層回歸效果影響較大。

圖 8 交叉驗證Fig. 8 Cross validationa. SE； b. 相關系數r2； c. SE預測值

SE預測值隨著取樣間距的增大而增大，說明預測誤差逐漸增大。定義：

SE預測值變化率=

(8)

第Ⅰ層SE預測值變化率分別為77.38%、54.7%、45.77%和24.48%； 第Ⅱ層SE預測值變化率分別為119.11%、25.45%、16.25%和12.56%； 第Ⅲ層SE預測值變化率分別為25.32%、20.01%、7.33%和7.23%。同一取樣間距下，地層厚度變化較劇烈的地層SE預測值變化率大于地層厚度變化不大的地層。對于地層厚度變化劇烈的地層，減小取樣間距可以有效的減少插值誤差，而對于地層厚度變化不大的地層，減小取樣間距對插值精度提高的意義不大。不同取樣間距的插值結果表明，取樣間距在一定程度上會影響地層厚度的插值精度。通過交叉驗證數據比較分析，可以看到，合理選擇取樣間距有利于提高預測結果的可靠性。取樣間距的選取取決于使用者對精度的要求。對于本算例，當選擇r2不小于0.8且SE預測值不大于1作為精度要求，對于土層厚度變化較劇烈的土層(第Ⅰ、Ⅱ層)，取樣間距控制在100im以內； 對于土層變化不大的土層(第Ⅲ層)，本文的最大取樣間距150im仍可以滿足要求。

5 結 論

本文基于GS+8.0軟件，以淺剖數據為源數據，鉆孔實測數據為驗證數據，采用克里金法對海洋地層厚度進行空間插值。通過比較不同變異函數模型和不同采樣間距的結果，得到以下結論：

(1)海洋地層厚度在空間分布上具有一定結構性與局部的變異性。因此，可將淺剖數據得到的地層厚度作為區域化變量，采用普通克里金進行空間插值。插值結果表明，普通克里金是一種有效的海底地層厚度預測方法。

(2)利用不同變異函數理論模型對實際半變異函數擬合，在結構變異分析的基礎上，進行交叉驗證獲得模型精度評價，并將鉆孔K1實測數據作為驗證。結果表明，決定系數最大，殘差最小的模型不一定是誤差最小的模型，應對不同模型下的插值結果進行綜合分析來選擇最合適的模型。在本文的井場地層厚度估計中，球狀模型為最優模型，高斯模型次之。

(3)取樣間距對插值精度有一定的影響。隨著取樣間距增大，插值擬合精度越低，預測誤差越大。對于球狀模型，增大取樣間距對地層厚度變化劇烈的地層回歸效果影響較小，對于地層厚度變化不大的地層回歸效果影響較大。同時，減小取樣間距可有效減小地層厚度變化劇烈的地層的插值誤差，而對地層厚度變化不大的地層插值精度的提高意義不大。