楊 陽 李 颯 何福耀 施展玲

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

0 引 言

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

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

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

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

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

1 地統(tǒng)計學(xué)方法

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

1.1 變異結(jié)構(gòu)分析

半變異函數(shù)是計算研究區(qū)域內(nèi)地層厚度變化的關(guān)鍵。關(guān)于地層厚度的理論半變異函數(shù)定義為地層厚度Z(x)在點(diǎn)x和x+h處數(shù)值之差的方差之半，其表達(dá)式為：

(1)

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

(2)

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

(3)

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

(4)

根據(jù)最小方差條件，可得：

(5)

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

(6)

1.2 半變異函數(shù)理論模型

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

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

2 地層的厚度

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

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

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

(7)

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

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

3 地層數(shù)據(jù)分析

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

3.1 半變異函數(shù)擬合

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

圖 2 擬合半變異函數(shù)圖Fig. 2 Fitted semivariograma. 第Ⅰ層； b. 第Ⅱ?qū)樱?c. 第Ⅲ層

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

表 2 變異函數(shù)參數(shù)統(tǒng)計Table 2 Variance function parameter statistics

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

表 3 交叉驗(yàn)證結(jié)果分析Table 3 Analysis of cross validation results

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

3.2 交叉驗(yàn)證

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

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

圖 4 第Ⅱ?qū)硬逯到Y(jié)果Fig. 4 Ordinary Kriging of layer Ⅱa. 指數(shù)模型； b. 球狀模型； c. 高斯模型

表3中，SE是回歸系數(shù)的標(biāo)準(zhǔn)誤差，其值越小說明精度越高；r2是相關(guān)系數(shù)的平方，其值越大說明回歸效果越好；SE預(yù)測值是標(biāo)準(zhǔn)誤差預(yù)測值，其值越大說明誤差越大。通過這3個指標(biāo)可以綜合評價不同半變異理論模型的預(yù)測精度。從表3中可以看到，在第Ⅰ、Ⅱ?qū)又校?種模型的r2接近，較難看出3種模型回歸效果的優(yōu)劣，但從SE預(yù)測值中可以看出球狀模型的誤差明顯小于其余兩種模型。在模型數(shù)據(jù)交叉驗(yàn)證環(huán)節(jié)，球狀模型表現(xiàn)出的整體精度最高，高斯模型次之，指數(shù)函數(shù)的誤差最大?？梢?，結(jié)構(gòu)分析最佳的模型不一定是誤差最小的模型。在出現(xiàn)結(jié)構(gòu)分析和誤差分析的擇優(yōu)結(jié)果不一致的情況下，很難說明球狀模型和高斯模型哪一個是絕對最優(yōu)。因此，需對插值效果進(jìn)一步對比分析以選擇最適合的模型。

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

3.3 地層厚度插值效果分析

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

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

利用研究區(qū)域內(nèi)的淺鉆孔K1(非采樣點(diǎn))的實(shí)測厚度數(shù)據(jù)和利用普通克里金法插值得出的地層厚度進(jìn)一步檢驗(yàn)(表 4)。從圖中可發(fā)現(xiàn)，不同理論模型插值的結(jié)果在淺鉆實(shí)測點(diǎn)均存在不同程度的誤差。其中，球狀模型的最小，在3組地層厚度誤差分別是1.9%、4.4%和2.3%。在地層厚度變化幅度較大的第Ⅰ、Ⅱ?qū)又杏葹槊黠@。由以上分析可知，球狀模型在該井場地層厚度離散數(shù)據(jù)插值中具有較高的實(shí)用性。

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

圖 6 最優(yōu)半變異函數(shù)模型選擇流程圖Fig. 6 Flow diagram of selection of optimal Semi-variogram function

4 取樣間距的影響

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

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

圖 8 交叉驗(yàn)證Fig. 8 Cross validationa. SE； b. 相關(guān)系數(shù)r2； c. SE預(yù)測值

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

SE預(yù)測值變化率=

(8)

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

5 結(jié) 論

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

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

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

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