張龍飛

摘 要：土層剪切波速是重要的工程地質(zhì)參數(shù)，對場地類別劃分及工程地質(zhì)條件分析具有重要意義。以大同盆地為研究對象，對盆地進行了地質(zhì)單元劃分及地層巖性特征分析，依據(jù)盆地內(nèi)已實施的大量地質(zhì)勘察鉆孔的剪切波速資料，采用地統(tǒng)計學(xué)方法中的空間分析、探索性統(tǒng)計分析、Kriging插值、確定性插值方法，給出了大同盆地范圍內(nèi)任意場點的Vse20及建筑場地類別預(yù)測分布。分析結(jié)果表明大同盆地內(nèi)區(qū)域范圍內(nèi)任意場點的等效剪切波速及建筑場地類別數(shù)據(jù)與地質(zhì)單元類型具有明顯的相關(guān)性，對于盆地內(nèi)不具備實施鉆孔及無法進行原位測試試驗的場點來說，本文提供了一種很好的預(yù)測分析方法，研究成果可為當(dāng)?shù)氐膰临Y源規(guī)劃及巖土勘察設(shè)計提供基礎(chǔ)資料，同時可為大同盆地區(qū)域性地質(zhì)條件分析、三維地質(zhì)模型及大震情景模型構(gòu)建提供一定的借鑒意義。

關(guān)鍵詞：大同盆地;等效剪切波速Vse20;建筑場地類別;地統(tǒng)計;插值

Prediction of? Vse20 and Building Site Classification in Datong Basin based

on Geostatistical Analysis

ZHANG Longfei

（Shanxi Institute of Seismic Engineering Investigation， Taiyuan 030002）

Abstract： Shear wave velocity of soil layer is an important engineering geological parameter， which is of great significance for site classification and engineering geological condition analysis. Taking Datong Basin as the research object， this paper divides the geological units and analyzes the stratigraphic and lithologic characteristics of the basin. Based on the data of shear wave velocity of the boreholes that have been implemented in the basin， through the spatial analysis， exploratory statistical analysis， Kriging interpolation and deterministic interpolation methods in geostatistics， the Vse20 of any field point and the prediction of building site types in Datong Basin are given. The results show that the data of equivalent shear wave velocity and construction site type of any field point in Datong Basin have obvious correlation with the type of geological unit. For the field points without drilling and in-situ test， this paper provides a good prediction and analysis method. The research results can be used for local land and resources planning and rock resources planning. The basic data of soil investigation and design can be used for reference in regional geological condition analysis， three-dimensional geological model and large earthquake scenario model construction of Datong Basin.

Keywords： Datong Basin; Equivalent shear wave velocity Vse20; Construction site category; Geostatistics; Interpolation

0 引言

土層剪切波速是重要的場地工程地質(zhì)條件參數(shù)，其能夠表征場地地層的軟硬狀態(tài)及相對密實程度，進行土層剪切波速研究對場地類別劃分、地質(zhì)條件分析及地質(zhì)模型建立都具有重要意義（彭艷菊等，2009）。通常土層剪切波速數(shù)據(jù)是通過鉆孔及孔內(nèi)原位試驗來實現(xiàn)的，然而對于大范圍的區(qū)域性場地來說，這種方法顯然是不現(xiàn)實的，為此我們需要尋求一種基于現(xiàn)有場地剪切波速及場地類別數(shù)據(jù)來預(yù)測估算未知場點數(shù)據(jù)的方法。而地統(tǒng)計學(xué)及其插值方法給我們提供了一種能夠?qū)崿F(xiàn)上述目的的途徑。

地統(tǒng)計學(xué)（又稱地質(zhì)統(tǒng)計學(xué)）是在傳統(tǒng)的統(tǒng)計學(xué)基礎(chǔ)上，以區(qū)域化變量為基礎(chǔ)，采用變異函數(shù)為主要工具，研究在空間分布上具有自相關(guān)性的自然現(xiàn)象的科學(xué)（呂連宏等，2006）。地統(tǒng)計是統(tǒng)計的一類，用于分析和預(yù)測與空間或時空現(xiàn)象相關(guān)的數(shù)值。它將數(shù)據(jù)的空間坐標(biāo)納入分析中。其不僅能夠提供插值，還可以衡量所插入的值的不確定性。地統(tǒng)計分析已從一元演化為多元，并提供了可融入用于補充主要目標(biāo)變量的輔助數(shù)據(jù)集的機制，從而可以構(gòu)建更準確的插值和不確定性模型。目前地統(tǒng)計學(xué)已在采礦、環(huán)境、土壤、氣象、地理及地質(zhì)等許多科學(xué)和領(lǐng)域中得到了廣泛應(yīng)用（梅志雄等，2008;陳勇等，2011;劉穎等，2011;晉銳等，2012;楊永川等，2012;熊林華等，2015;徐占軍等，2018;郎藝超等，2018）。

大同盆地近年來已經(jīng)實施了大量的地質(zhì)勘察鉆孔，積累了豐富的鉆孔勘察及測試數(shù)據(jù)。本文以大同盆地為研究對象，對盆地進行了區(qū)域地質(zhì)單元劃分及地層特征分析，對地統(tǒng)計學(xué)中的空間分析、探索性統(tǒng)計分析、Kriging插值、確定性插值方法進行了簡要介紹，并以大同盆地已有鉆孔數(shù)據(jù)為基礎(chǔ)，基于地統(tǒng)計學(xué)方法預(yù)測分析了整個大同盆地范圍內(nèi)任意場點的Vse20及建筑場地類別分布，預(yù)測結(jié)果表明，盆地內(nèi)的地質(zhì)單元類型與Vse20及建筑場地類別具有很好的相關(guān)性，在無任何測試數(shù)據(jù)時，可基于本文研究成果對場等效剪切波速及場地類別指數(shù)進行預(yù)測。研究成果可以為大同盆地三維地質(zhì)模型及大震情景構(gòu)建提供基礎(chǔ)數(shù)據(jù)，并可作為大同地區(qū)建筑場地類別劃分及城市規(guī)劃的重要參考依據(jù)。

1 大同盆地區(qū)域地質(zhì)背景

（1）新構(gòu)造單元劃分

大同盆地雛形形成于中新世時期，上新世以來受本區(qū)地殼整體抬升及盆地主控邊界斷裂影響形成不對稱箕狀斷陷盆地。盆地總體走向NE-NNE，長約225km，寬約60km，受盆地主控邊界斷裂及盆地隱伏斷裂控制，盆地內(nèi)發(fā)育了5個次級凹陷構(gòu)造單元（圖1），其中馬營莊凹陷第四系最大沉積厚度為900m，是大同斷陷盆地的沉降中心，是大同盆地斷陷最深的構(gòu)造單元;朔州斷階第四系最大沉積厚度為100m，是大同盆地凹陷區(qū)斷陷最淺的構(gòu)造單元。

（2）第四系特征

大同盆地內(nèi)廣泛堆積覆蓋著第四紀沉積地層，其第四系沉積物的厚度、地層結(jié)構(gòu)及分布規(guī)律受盆地基底形態(tài)、盆地斷裂的共同控制。從盆地周邊的丘陵邊山區(qū)、洪積扇及傾斜平原區(qū)向盆地中心過渡為沖積平原及湖積平原區(qū)，表現(xiàn)為明顯的同時異相特征，總體上盆地中心粉土、粉質(zhì)黏土分布較多，由盆地中心向盆地邊緣過渡砂礫石粒徑逐漸增大，砂礫石層逐漸增厚，直至山前洪積扇地帶轉(zhuǎn)變?yōu)槁训[石層。根據(jù)盆地地形地貌及地層巖性特征，可將大同盆地劃分為Qhal-l、Qp3al-l、Qhal+pl、Qp3al、Qp3pl五類地質(zhì)單元。

（3）勘察鉆孔

通常場地剪切波速與場地類別數(shù)據(jù)是通過鉆孔及孔內(nèi)原位測試的技術(shù)方法來獲取的，作者基于地震安全性評價、地震小區(qū)劃、工程地質(zhì)勘察等項目，收集了近十年來于大同盆地范圍內(nèi)已經(jīng)實施的地質(zhì)鉆孔，共計1488個，由圖1可以看出大多數(shù)鉆孔位于城區(qū)及其周邊，城區(qū)以外的地質(zhì)鉆孔密度較小，分布較為稀疏，若要獲取區(qū)域范圍內(nèi)大面積的場地剪切波速與建筑場地類別數(shù)據(jù)，通過鉆孔來一一實現(xiàn)顯然是不可能的。為此，我們需要尋求一種基于區(qū)域范圍內(nèi)已知鉆孔場點數(shù)據(jù)獲取區(qū)域內(nèi)未知場點數(shù)據(jù)的估算方法，而地統(tǒng)計學(xué)方法能夠很好的實現(xiàn)上述目的。

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

地統(tǒng)計學(xué)是在傳統(tǒng)的統(tǒng)計學(xué)基礎(chǔ)上，以區(qū)域化變量為基礎(chǔ)，采用變異函數(shù)為主要工具，研究在空間分布上具有自相關(guān)性的自然現(xiàn)象的科學(xué)（李鐘山，1997）。本文綜合運用地統(tǒng)計學(xué)中的空間分析、探索性統(tǒng)計分析、Kriging插值、確定性插值的方法對大同盆地鉆孔測試數(shù)據(jù)進行分析，進而獲取鉆孔場點及未知場點的Vse20空間分布特征情況，現(xiàn)對其分析方法作簡要介紹：

2.1 空間分析

（1）全局空間自相關(guān)

全局空間自相關(guān)是用來描述目標(biāo)要素在大范圍區(qū)域內(nèi)的空間分布特征的參數(shù)，通常需要根據(jù)目標(biāo)要素的空間位置及其屬性值來進行衡量，進而給出二者的分布模式，常用的分布模式有聚類、離散及隨機模式。一般可以用Moran's I （莫蘭指數(shù)）來表示全局空間的自相關(guān)性（劉會，2017）。實際使用中Moran's I 指數(shù)需進行標(biāo)準化，其計算公式如下：

式中：n為研究對象的數(shù)目，Xi為觀測值，為Xi的平均值。當(dāng)I>0時表明目標(biāo)參數(shù)的空間位置及其屬性值存在空間正相關(guān)，當(dāng)I<0時表明目標(biāo)參數(shù)的空間位置及其屬性值存在空間負相關(guān)性，I=0時則表明二者之間不存在空間相關(guān)性，

（2）局域自相關(guān)檢驗

目標(biāo)要素在區(qū)域范圍內(nèi)既存在空間同質(zhì)性，又存在空間異質(zhì)性特征，而局域自相關(guān)檢驗就是為了衡量每個目標(biāo)空間要素屬性的空間異質(zhì)性而定義的參數(shù)。通?？捎肎 統(tǒng)計量（Getis'G）來衡量局域空間自相關(guān)性，其計算式為：

ZG得分計算式：

式中：

若G（d）取高值則為高值集聚，反之則為低值集聚。

2.2 探索性統(tǒng)計分析

探索性統(tǒng)計分析是指通過直方圖、正態(tài) QQ 圖、趨勢分析、Voronoi圖等方法來查找區(qū)域范圍內(nèi)的目標(biāo)要素異常值，檢查數(shù)據(jù)的局部變化，查找全局趨勢特征。其中直方圖主要用于顯示目標(biāo)參數(shù)的頻率分布;正態(tài) QQ 圖用于評估兩組數(shù)據(jù)之間的分布相似程度;趨勢分析主要用于識別輸入數(shù)據(jù)的集中分布趨勢;Voronoi地圖主要是用于計算多種局部統(tǒng)計量;半變異函數(shù)/協(xié)方差云可用來查找局部異常值及檢查數(shù)據(jù)集中空間自相關(guān)局部特征;交叉協(xié)方差云可用來檢查兩個數(shù)據(jù)集之間空間相關(guān)的局部特征。

2.3 Kriging插值

克里金（Kriging）插值法是從已有目標(biāo)要素出發(fā)，根據(jù)已有目標(biāo)點和預(yù)測點之間的空間位置關(guān)系來進行自相關(guān)分析，給出預(yù)測值的插值估計及方差估計現(xiàn)象（劉愛利，2012）。克里金方法能將概率模型與預(yù)測場點屬性值建立相關(guān)性，在地統(tǒng)計中預(yù)測變量可用下式來表示：

式中：Z（S）是預(yù)測變量，μ（S）是確定性趨勢，ε（S）是隨機自相關(guān)誤差。

克里金法是一種能實現(xiàn)精確或平滑插值的插值方法，通常用于克里金插值分析的數(shù)據(jù)集需符合正態(tài)分布，若數(shù)據(jù)不符合正太分布，則需通過探索性統(tǒng)計分析剔除異常偏離數(shù)據(jù)及采用對數(shù)變換的方法進行調(diào)整。

2.4 確定性插值

當(dāng)數(shù)據(jù)不滿足克里金假設(shè)條件，即不符合正態(tài)分布時，就需要用到確定性插值方法，確定性方法是一種不存在空間自相關(guān)性的方法，其并不符合隨機模型。常用的確定性方法包括全局多項式插值法、局部多項式插值法、反距離權(quán)重法、徑向基函數(shù)插值法及含障礙的插值法，本文主要采用徑向基函數(shù)插值法（RBF），RBF是一系列精確插值方法的組合，其可用于區(qū)域范圍內(nèi)由大量已知測試點生成區(qū)域內(nèi)分布平面的情況。

3 空間分析及結(jié)果

3.1 基于地統(tǒng)計分析的Vse20估計

（1）空間分析結(jié)果

對大同盆地所有鉆孔場點的測試數(shù)據(jù)進行空間自相關(guān)處理，結(jié)果如下：I= 0.0589 > 0，大于其期望值（E= -0.0008）;并且 Z = 44.1399，P值<0.0001，表明隨機分布結(jié)果的可能性小于 1%。Z得分和 P值表明可以拒絕零假設(shè)，這就意味著Vse20表現(xiàn)出統(tǒng)計意義上的顯著性聚類或離散模式，而不是隨機模式，即距離較近的場點剪切波速Vse20也相應(yīng)接近。對數(shù)據(jù)進行高低聚類處理獲得結(jié)果：G值為0.3631，大于其期望值（E=0.3585），Z值為4.9362，表示剪切波速Vse20有高值集聚情況。

（2）探索性統(tǒng)計分析結(jié)果

圖2（a）為剔除大同盆地334個異常鉆孔數(shù)據(jù)點后對剩余1154個剪切波速場點數(shù)據(jù)進行的原始直方圖分析結(jié)果，樣本數(shù)據(jù)最大為512.49m/s，最小為185.56m/s，均值為263.59m/s，標(biāo)準差為34.46m/s，偏度（Skewness）為2.04，峰度（Kurtosis）值為14.53，說明相比于標(biāo)準正態(tài)分布，該數(shù)據(jù)的分布偏高聳且狹窄，即屬于正偏態(tài);圖2（b）是經(jīng)過對數(shù)log變換后的數(shù)據(jù)直方圖，偏度（Skewness）為0.77，峰度（Kurtosis）值為7.27，樣本數(shù)據(jù)最大為6.24，最小為5.22，均值為5.57，標(biāo)準差為0.12，表明經(jīng)過變換后的數(shù)據(jù)比原始數(shù)據(jù)更符合正態(tài)分布。

正態(tài) QQ 圖上的點可指示數(shù)據(jù)集的單變量分布的正態(tài)性。如果數(shù)據(jù)是正態(tài)分布的，數(shù)據(jù)點將落在45°參考線上。如果數(shù)據(jù)不是正態(tài)分布的，數(shù)據(jù)點將會偏離45°參考線。圖3為大同盆地1154個鉆孔場點的剪切波速數(shù)據(jù)正態(tài) QQ 圖，圖3（a）為原始分析結(jié)果，圖3（b）為經(jīng)過對數(shù)變換后的結(jié)果，可以看出對數(shù)變換后的數(shù)據(jù)集偏離度更小，有利于后續(xù)克里金插值分析。

利用地統(tǒng)計分析模塊中趨勢分析功能，可以生成大同盆地Vse20的空間變化趨勢圖（圖4），圖中中間部分為散點，兩個側(cè)面分別表示樣本點在不同方位的投影，其中Y軸反映的是南北走向的變化趨勢，X軸為東西走向的變化趨勢，Z軸反映的是Vse20的大小。從圖4可以看出，SN軸與WE 軸的多項式曲線均呈現(xiàn)U字型趨勢，即Vse20從區(qū)域中心向各個邊緣呈現(xiàn)遞增趨勢，且南北走向的遞增趨勢比東西走向的遞增趨勢更明顯。

（3）Kriging插值分析結(jié)果

由于大同盆地Vse20數(shù)據(jù)進行對數(shù)變換后服從正態(tài)分布，因此本文選擇簡單克里金插值方法基于已有鉆孔場點測試數(shù)據(jù)進行空間插值分析，插值結(jié)果見圖5，分析可知大同盆地中部沖湖積平原區(qū)的Vse20基本位于188.4～249.5m/s，沖積平原區(qū)的Vse20基本位于249.4～269.3m/s，盆地邊緣洪積平原區(qū)的Vse20基本位于269.3～483.6m/s，盆地等效剪切波速數(shù)值與地質(zhì)單元類型具有明顯的關(guān)聯(lián)性，即每類地質(zhì)單元均有其優(yōu)勢剪切波速分布區(qū)間，對于無剪切波速測試數(shù)據(jù)的空間場點而言可根據(jù)地質(zhì)單元類型大致估算出該點的優(yōu)勢剪切波速分布數(shù)值，對于本文而言，給定任意場點地理位置，可直接由圖5求得該場點的Vse20，且預(yù)測精確，使用方便。

3.2基于確定性插值分析的場地類別估計

當(dāng)統(tǒng)計數(shù)據(jù)不滿足克里金假設(shè)條件，就需要采用確定性插值方法，確定性方法并不基于隨機空間過程模型，且數(shù)據(jù)中不存在空間自相關(guān)的顯式測量或建模。由于大同盆地任意場點的建筑場地類別不符合正態(tài)分布，我們選取了徑向基函數(shù)（RBF）插值法來對區(qū)域性場地點的場地類別進行預(yù)測分析，RBF是一個中等速度精確的確定性插值器，其可提供與克里金法的精確形式相當(dāng)?shù)念A(yù)測表面，更具自動化，且不對數(shù)據(jù)進行任何假設(shè)。

圖6為大同盆地徑向基函數(shù)法（RBF）空間插值分析結(jié)果，由圖6大同盆地建筑場地類別與構(gòu)造單元及地質(zhì)單元類型有很好的關(guān)聯(lián)性，其中僅在管涔山及采涼山基巖隆起山區(qū)前緣可見場地類別指數(shù)1.0～1.8的區(qū)間段出露，其分布特征與地質(zhì)單元及鉆孔測試結(jié)果相吻合，綜合其分布特征可將盆地邊緣1.0～1.8的區(qū)間段定義為建筑場地類別Ⅰ1類，同時本次統(tǒng)計未見盆地內(nèi)有Ⅰ0類場地出露;場地類別指數(shù)1.8～2.2的區(qū)間段對應(yīng)于盆地內(nèi)Qp3al及Qp3pl地質(zhì)單元，本次統(tǒng)計研究表明上述兩個地質(zhì)單元內(nèi)鉆孔場地類別均為Ⅱ類，故此可將1.8～2.2的區(qū)間段定義為Ⅱ0類區(qū)間段;場地類別指數(shù)2.2～2.5的區(qū)間段對應(yīng)盆地內(nèi)Qhal+pl地質(zhì)單元，本次統(tǒng)計研究表明該地質(zhì)單元內(nèi)絕大多數(shù)鉆孔場點的建筑場地類別為Ⅱ類，僅有少數(shù)鉆孔場點為Ⅲ類，本文將1.8～2.2的區(qū)間段定義為Ⅱ1類區(qū)間段;場地類別指數(shù)2.5～2.8的區(qū)間段對應(yīng)盆地中部外圍Qhal地質(zhì)單元，本次統(tǒng)計研究表明該地質(zhì)單元內(nèi)絕大多數(shù)鉆孔場點的建筑場地類別為Ⅲ類，僅有少數(shù)鉆孔場點為Ⅱ類，本文將2.5～2.8的區(qū)間段定義為Ⅲ0類區(qū)間段;場地類別指數(shù)2.8～3.0的區(qū)間段對應(yīng)盆地中部Qhal-l地質(zhì)單元，本次統(tǒng)計研究表明該地質(zhì)單元內(nèi)鉆孔場點的建筑場地類別均為Ⅲ類，故可將2.5～2.8的區(qū)間段定義為Ⅲ1類區(qū)間段，場地分類特征詳見表1?；诖朔诸惙椒梢钥焖俟浪愦笸璧氐慕ㄖ龅仡悇e特征。

4 結(jié)論

（1）分析了大同盆地的地質(zhì)構(gòu)造背景及地層巖性特征，對盆地進行了新構(gòu)造單元及地質(zhì)地貌單元劃分，將大同盆地劃分為5個構(gòu)造單元及五類地質(zhì)單元。

（2）搜集了整個大同盆地內(nèi)已經(jīng)實施完成的區(qū)域地質(zhì)鉆孔及其剪切波速測試數(shù)據(jù)，對已有場點的波速數(shù)據(jù)進行了簡要分析，為了估算整個大同盆地內(nèi)任意場點的等效剪切波速及場地類別數(shù)據(jù)，引入了地統(tǒng)計學(xué)方法。

（3）基于地統(tǒng)計分析，依據(jù)盆地內(nèi)已知的鉆孔場點等效剪切波速及場地類別數(shù)據(jù)，采用空間分析、探索性分析、克里金插值及確定性插值方法給出了大同盆地的Vse20及建筑場地類別的預(yù)測分布。

（4）大同盆地內(nèi)區(qū)域性場點的等效剪切波速及建筑場地類別數(shù)據(jù)與地質(zhì)單元類型具有很好的相關(guān)性，對于盆地內(nèi)無法實施鉆孔及原位測試試驗的場點來說，本文提供的預(yù)測分析方法，可為當(dāng)?shù)氐膰临Y源規(guī)劃及巖土勘察設(shè)計提供參考借鑒，對大同盆地的區(qū)域地質(zhì)條件分析、三維地質(zhì)模型及大震情景模型構(gòu)建具有重要的科學(xué)意義。

參考文獻：

陳勇， John M C， Dogan T， 2011. 基于特征價格模型的住宅需求價格彈性分析： 深圳住宅市場實證研究[J]. 城市發(fā)展研究， 18（2）：? 62-67.

晉銳， 李新， 閻保平，等， 2012. 黑河流域生態(tài)水文傳感器網(wǎng)絡(luò)設(shè)計[J]. 地球科學(xué)進展， 27（9）： 993-1005.

李鐘山， 1997. 地質(zhì)統(tǒng)計學(xué)中的區(qū)域化變量理論[J]. 世界地質(zhì)（2）： 85-93.

呂連宏， 張征， 遲志淼， 等， 2006. 地質(zhì)統(tǒng)計學(xué)在環(huán)境科學(xué)領(lǐng)域的應(yīng)用進展[J]. 地球科學(xué)與環(huán)境學(xué)報， 28（1）： 105-109.

劉穎， 張平宇， 李靜， 2011.長春市區(qū)新建住宅價格的空間格局分析[J]. 地理科學(xué)， 31（1）： 95-101.

劉愛利， 2012. 地統(tǒng)計學(xué)概論[M]. 北京： 科學(xué)出版社： 50-76.

劉會， 2017. 當(dāng)代中國農(nóng)村土地流轉(zhuǎn)的工業(yè)條件研究： 基于全局莫蘭指數(shù)與空間計量模型的研究[J]. 財經(jīng)理論研究（6）： 20-29.

郎藝超， 肖璐， George Christakos， 2018. 基于 SARIMA 模型和普通 Kriging 法對杭州市主城區(qū) PM_（2.5）短期預(yù)測和制圖[J]. 環(huán)境科學(xué)學(xué)報， 38（1）： 62-70.

梅志雄， 黎夏， 2008. 基于 ESDA 和 Kriging方法的東莞市住宅價格空間結(jié)構(gòu)[J]. 經(jīng)濟地理， 28（5）： 862-866.

彭艷菊， 呂悅軍， 黃雅虹， 等， 2009. 工程地震中的場地分類方法及適用性評述[J].地震地質(zhì)， 31（2）： 349-362.

熊林華， 張軍， 吳健平， 2015. 基于房產(chǎn)網(wǎng)站數(shù)據(jù)的商品住宅價格空間分布研究： 以成都市為例[J]. 測繪與空間地理信息， 38（9）： 150-154.

楊永川， 楊軻， 王志浩，等， 2012. 空間插值法在熱環(huán)境流動觀測中的應(yīng)用[J]. 中南大學(xué)學(xué)報（自然科學(xué)版）， 43（9）： 422-429.

徐占軍， 張媛， 張紹良， 等， 2018. 基于GIS與分區(qū)Kriging的采煤沉陷區(qū)土壤有機碳含量空間預(yù)測[J]. 農(nóng)業(yè)工程學(xué)報， 34（10）： 253-259.