基于普通克里金法的同位素測氡探火數(shù)據(jù)優(yōu)化處理

劉 軒，王俊峰，周 斌，劉 碩

(太原理工大學(xué) a.安全與應(yīng)急管理工程學(xué)院，b.山西省礦井通風(fēng)與火災(zāi)防治工程技術(shù)研究中心，c.礦業(yè)工程學(xué)院，太原 030024)

礦井自燃火災(zāi)防治的關(guān)鍵在于火源位置的準(zhǔn)確探測，基于不同的探測原理，國內(nèi)外研究學(xué)者提出了不同探測方法，包括地表測溫法、遙感法、同位素測氡法等[1]。其中，同位素測氡法是利用煤巖介質(zhì)中天然放射性氡隨溫度升高析出率增加的特性，在地面探測氡及其子體的變化規(guī)律，并經(jīng)過一系列數(shù)據(jù)分析處理，從而給出火源位置、范圍及發(fā)展趨勢。大量的滅火工程實(shí)踐證明，該方法能夠克服遙感法分辨率低及地表測溫法反演計(jì)算復(fù)雜、不適合深部火區(qū)探測的缺點(diǎn)，是一種行之有效的非接觸隱蔽火源探測技術(shù)[2-3]。

等值線圖作為一種直觀的圖件，廣泛應(yīng)用于各種數(shù)據(jù)的分析處理之中，對于同位測氡法在地表有限測點(diǎn)所獲得的氡濃度數(shù)據(jù)，其生成等值線圖的關(guān)鍵是插值方法的選擇。目前，主流分析軟件中使用的插值方法是最小曲率法，該法通過構(gòu)造穿過平面每一個(gè)數(shù)據(jù)點(diǎn)并具有最小曲率的曲線，進(jìn)而組成氡濃度分布等值線圖，并將圖中的高值區(qū)域確定為火源位置。由于最小曲率法主要考慮曲線的光滑性加之測點(diǎn)數(shù)據(jù)量有限，因此形成的氡濃度分布失真，往往超過最大值和最小值的范圍，與實(shí)際相差較大，進(jìn)而影響火源位置的圈定[4]。事實(shí)上，地表氡具有空間分布的特點(diǎn)，即取值與所在區(qū)域內(nèi)的位置有關(guān)且相鄰測點(diǎn)間具有某種程度的自相關(guān)性，屬于區(qū)域化隨機(jī)變量的范疇，僅依靠最小曲率法無法在充分挖掘區(qū)域化隨機(jī)變量所蘊(yùn)含的結(jié)構(gòu)性信息的基礎(chǔ)上進(jìn)行插值[5]。

目前，常用的單變量克里金插值方法包括簡單克里金法、普通克里金法以及泛克里金法[6]。其中，簡單克里金法假設(shè)區(qū)域化變量的數(shù)學(xué)期望是已知的常數(shù)，在實(shí)際應(yīng)用中一般難以滿足[7]；泛克里金法適用于存在某種漂移的區(qū)域化變量，該方法將插值過程分為趨勢項(xiàng)和殘差項(xiàng)兩部分之和，但殘差變異函數(shù)的預(yù)測較為繁瑣，適用于煤礦可采儲量預(yù)測等問題[8]；普通克里金法要求所研究區(qū)域化變量的數(shù)學(xué)期望為未知常數(shù)，該方法是目前應(yīng)用最為廣泛的插值方法，具有線性無偏、最優(yōu)估計(jì)的特點(diǎn)，近年來已廣泛應(yīng)用于表層土壤磁學(xué)性質(zhì)的研究[9]、地下水位和土壤濕度的采樣分析[10]、淺源地震的重力變化特征延拓[11]以及礦井水文地質(zhì)建模[12]等領(lǐng)域的研究中，均顯示了良好的應(yīng)用效果。

地理信息系統(tǒng)(geographic information system，GIS)技術(shù)中地統(tǒng)計(jì)分析模塊集成普通克里金插值方法所需的探索性空間數(shù)據(jù)分析、變異函數(shù)建模等功能，是實(shí)現(xiàn)普通克里金插值方法的技術(shù)載體[13]。

筆者擬將地質(zhì)統(tǒng)計(jì)學(xué)中普通克里金插值方法應(yīng)用于同位素測氡探火數(shù)據(jù)的處理中，并將插值結(jié)果與最小曲率法所得結(jié)果進(jìn)行比較，最終通過鉆孔驗(yàn)證火源位置圈定的準(zhǔn)確性，以期優(yōu)化氡濃度數(shù)據(jù)的處理方法，最大程度提高同位素測氡法判斷火源位置的準(zhǔn)確度。

1 普通克里金法的理論基礎(chǔ)

普通克里金法是以區(qū)域化變量理論和空間連續(xù)性理論為基礎(chǔ)，以變異函數(shù)為基本工具，研究在空間分布上既有結(jié)構(gòu)性又有隨機(jī)性的地質(zhì)變量，建立符合地質(zhì)規(guī)律的統(tǒng)計(jì)模型來反映變量的變化規(guī)律，并對其空間分布進(jìn)行預(yù)測[14]。

1.1 區(qū)域化變量

區(qū)域化變量Z(x)是以空間點(diǎn)x的3個(gè)直角坐標(biāo)xu，xv，xw為自變量的隨機(jī)場Z(xu，xv，xw).從其定義可以看出，區(qū)域化變量描述的現(xiàn)象具有空間分布的特點(diǎn)，地表氡濃度因其取值是與位置有關(guān)的隨機(jī)函數(shù)，符合區(qū)域化變量的定義[15]。從另一個(gè)側(cè)面反映出僅通過純數(shù)學(xué)方法的最小曲率擬合，無法合理有效地預(yù)測地表氡濃度趨勢。而區(qū)域化變量所具有的結(jié)構(gòu)性和隨機(jī)性特征需要通過變差異函數(shù)來表征。

1.2 變異函數(shù)

變異函數(shù)是區(qū)域化變量增量的方差，將區(qū)域化變量Z(x)在x，x+h兩點(diǎn)處差值方差的一半定義為Z(x)在h方向上的變異函數(shù)γ(h)，其表達(dá)式為：

(1)

式中，Var代表方差。

而在實(shí)際工作中以試驗(yàn)變異函數(shù)來計(jì)算，其表達(dá)式見公式(2)：

(2)

式中：h為xi和xi+h兩點(diǎn)的距離；Z(xi)和Z(xi+h)分別為xi，xi+h兩點(diǎn)處的觀測+6值；n(h)為相距h的數(shù)據(jù)對數(shù)目；γ*(h)為實(shí)驗(yàn)變異函數(shù)。

圖1所示為理想的試驗(yàn)變異函數(shù)圖。

圖1 變異函數(shù)圖Fig.1 Variogram function

圖中有3個(gè)基本參數(shù)：

1) 變程a，用來度量空間相關(guān)性的最大距離，是變異函數(shù)達(dá)到穩(wěn)定值時(shí)的空間距離。變程越小，區(qū)域化變量空間規(guī)律性變化的范圍越小，該參數(shù)在該方向上變化越快，即非均質(zhì)性越強(qiáng)；相反，變程越大，區(qū)域化變量空間規(guī)律性變化的范圍越大，表明該參數(shù)在該方向上變化越慢，也就是非均質(zhì)性越弱。

2) 基臺值(C0+C)，是變異函數(shù)在變程處達(dá)到的平穩(wěn)值，反映了研究范圍內(nèi)的變異強(qiáng)度?；_值越大，該參數(shù)在該方向上變化幅度越大，也就是非均質(zhì)性越強(qiáng)；基臺值越小，該參數(shù)在該方向上變化幅度越小，即非均質(zhì)性越弱。

3) 塊金值C0，是h=0時(shí)的變異函數(shù)值，反映了區(qū)域化變量的隨機(jī)性大小[16]。

1.3 理論模型擬合

根據(jù)實(shí)測數(shù)據(jù)做出的試驗(yàn)變異函數(shù)圖，實(shí)際上只是一條非光滑曲線，還需用一條適當(dāng)?shù)膱A滑曲線對它進(jìn)行擬合，并用特定的函數(shù)來表示擬合結(jié)果，即變異函數(shù)的理論模型，該模型將直接參與克里金插值。以最常用的球狀模型為例，其表達(dá)式：

(3)

球狀模型曲線如圖2所示。

圖2 球狀模型曲線Fig.2 Spherical model curve

1.4 普通克里金法

普通克里金插值法對于任意待估點(diǎn)的估計(jì)值是通過該點(diǎn)影響范圍內(nèi)的n個(gè)有效樣品值的線性組合得到，即

(4)

1) 無偏估計(jì)，即估計(jì)誤差ε的期望值為0，可表示為：

(5)

2) 估計(jì)方差最小，見式(6).

(6)

在滿足以上條件的同時(shí)，為保證權(quán)重系數(shù)λi有解，需要對區(qū)域化變量Z(x)的數(shù)學(xué)期望m的取值做出假設(shè)。不同克里金法的區(qū)別就在于所做假設(shè)的不同，普通克里金法要求Z(x)滿足二階平穩(wěn)假設(shè)，即數(shù)學(xué)期望為未知常數(shù)m.整理得普通克里金方程組見公式(7).

(7)

該方程組有n+1個(gè)未知數(shù)(n個(gè)權(quán)重系數(shù)λi和一個(gè)拉格朗日乘數(shù)μ)和n+1個(gè)方程，可從方程組解出未知數(shù)值。綜合以上步驟，普通克里金插值法的流程如圖3所示。

圖3 普通克里金插值法流程圖Fig.3 Ordinary Kriging interpolation flowchart

2 應(yīng)用實(shí)例

2.1 研究區(qū)域概況

本文選取晉城煤業(yè)集團(tuán)古書院煤礦張嶺火區(qū)，著火煤層為3#煤層，該煤層較為穩(wěn)定，厚度為0.6～0.9 m，埋藏較淺，僅為15～30 m；上覆巖層為厚約6 m左右的粉砂質(zhì)泥巖，其上為中細(xì)粒長石英砂巖，厚度為15～19 m，砂巖上部覆蓋有第四系上更新黃土。根據(jù)現(xiàn)場踏勘結(jié)果及井下采動(dòng)影響區(qū)域，確定探測區(qū)域面積50 850 m2，采用HOLUX GR-260型手持GPS進(jìn)行測點(diǎn)定位，以15 m為間距共布置234個(gè)測點(diǎn)，并使其均勻分布于整個(gè)探測區(qū)域，并在各個(gè)測點(diǎn)處使用α杯法對自燃危險(xiǎn)區(qū)域的位置與范圍進(jìn)行探測。α杯法是將高吸附材料制成的α探杯倒置埋入土壤中，待杯內(nèi)氡濃度與土壤環(huán)境中的氡濃度達(dá)到放射性動(dòng)態(tài)平衡后(一般認(rèn)為4 h以上即可)將探杯取出并放入CD-1α杯探測儀測量，由于吸附于探杯內(nèi)壁上的氡子體與氡濃度存在線性關(guān)系，故可通過測量氡子體實(shí)現(xiàn)對土壤中氡濃度的測量。

探測原始數(shù)據(jù)為每3 min內(nèi)氡及其子體衰變個(gè)數(shù)的計(jì)數(shù)值(N/3 min)，本區(qū)域氡射氣底數(shù)為50 N/3 min.

2.2 異常值分析

對于測得的原始數(shù)據(jù)，需要對數(shù)據(jù)中存在的由于儀器震動(dòng)、人為操作等造成的異常值進(jìn)行剔除，并分析其產(chǎn)生的原因，如果不加剔除地把這些異常值包括進(jìn)數(shù)據(jù)的計(jì)算分析過程中，會(huì)影響插值結(jié)果。

箱形圖因其判斷異常值的標(biāo)準(zhǔn)以四分位數(shù)和四分位距為基礎(chǔ)，異常值不能對這個(gè)標(biāo)準(zhǔn)施加影響，所以常被用來觀察數(shù)據(jù)中存在的異常值。箱形圖中上、下邊緣之間的部分稱為內(nèi)限，處于內(nèi)限以外位置的點(diǎn)表示的數(shù)據(jù)為異常值，這里需要說明的是異常值是對數(shù)據(jù)中離群值點(diǎn)的統(tǒng)稱，其所代表的意義是由使用箱形圖的目的所決定的，在對測氡探火數(shù)據(jù)的分析中，異常值除顯示需要剔除的測點(diǎn)外，還可表征火區(qū)范圍內(nèi)的測點(diǎn)。

這樣的判斷是基于煤自燃演化規(guī)律，由于自燃煤層長時(shí)間氧化蓄熱致使溫度上升，溫度從火源中心點(diǎn)呈擴(kuò)散狀向周圍環(huán)境依次遞減，相應(yīng)地表火區(qū)對應(yīng)位置處氡濃度也呈依次遞減的擴(kuò)散分布，在箱形圖上則表現(xiàn)為排列緊密的數(shù)據(jù)點(diǎn)[17-18]。圖中孤立異常值點(diǎn)所處的氡濃度區(qū)間為1 441 N/3 min～1 817 N/3 min，共6個(gè)點(diǎn)，分析其產(chǎn)生的原因如下：

1) 將探測杯從埋坑中鏟出的過程，會(huì)在周圍產(chǎn)生一系列的沖擊裂隙，在增加了土壤透氣性的同時(shí)為氡的擴(kuò)散提供了更多的運(yùn)移通道。

2) 同樣在該過程中，鏟出的部分土壤會(huì)破壞其與大氣之間已形成的壓力平衡，造成空氣對流而導(dǎo)致更多的氡涌向地面。

3) 儀器的晃動(dòng)影響CD-1α杯測氡儀的穩(wěn)定性。

以上過程均可使瞬時(shí)氡濃度急劇升高，導(dǎo)致孤立異常值點(diǎn)的產(chǎn)生，對這些測點(diǎn)重新測量后所形成的箱形圖(圖4)連續(xù)異常值濃度區(qū)間為464 N/3 min～743 N/3 min，即火區(qū)范圍內(nèi)測點(diǎn)的氡濃度所在區(qū)間。

圖4 原始數(shù)據(jù)箱形圖Fig.4 Radon counts boxplot

2.3 探索性數(shù)據(jù)分析

使用普通克里金法進(jìn)行插值分析，為滿足二階平穩(wěn)假設(shè)，首先要確保數(shù)據(jù)呈正態(tài)分布，否則需要根據(jù)分布情況對數(shù)據(jù)進(jìn)行轉(zhuǎn)換。

頻率分布直方圖因其能清楚顯示各組頻數(shù)分布情況又易于顯示各組之間頻數(shù)的差別，故可直觀地反映出數(shù)據(jù)的分布狀態(tài)。調(diào)用ArcGIS統(tǒng)計(jì)分析模塊中的探索性數(shù)據(jù)分析功能，得到地表氡濃度數(shù)據(jù)的頻數(shù)分布直方圖。觀察直方圖呈“陡壁型”分布，如圖5所示。表1為ArcGIS探索性數(shù)據(jù)分析所得的各項(xiàng)統(tǒng)計(jì)指標(biāo)。

圖5 氡濃度計(jì)數(shù)直方圖Fig.5 Radon counts histogram

表1 氡濃度分布統(tǒng)計(jì)指標(biāo)Table 1 Statistical indicators for radon counts

中位數(shù)138 N/3 min向左偏離平均數(shù)182 N/3 min，頻數(shù)自左至右逐漸減少，峰度為2.733，偏度為1.807，數(shù)據(jù)分布呈現(xiàn)不對稱狀態(tài)。

同時(shí)，Shapiro-Wilk檢驗(yàn)得到P<0.001，在5%的顯著性水平下拒絕原假設(shè)H0(H0：該隨機(jī)變量服從正態(tài)分布)，即不服從正態(tài)分布，因而無法滿足插值要求。

所有理論分布中，擬合度最高的是對數(shù)正態(tài)分布，決定系數(shù)R2為0.967 8.同時(shí)Anderson-Darling檢驗(yàn)得到P值為0.404，在5%的顯著性水平下并未拒絕原假設(shè)H0(該隨機(jī)變量服從對數(shù)正態(tài)分布)，故該礦所測得氡濃度呈對數(shù)正態(tài)分布，需要經(jīng)過對數(shù)轉(zhuǎn)換才可進(jìn)行普通克里金插值。

2.4 變異函數(shù)建模

對經(jīng)對數(shù)轉(zhuǎn)換后的氡濃度數(shù)據(jù)進(jìn)行變異函數(shù)建模，球狀模型因其決定系數(shù)R2為0.78，成為理論模型中對實(shí)驗(yàn)變異函數(shù)擬合度最優(yōu)的模型，如圖6所示。表2為球狀模型變異函數(shù)曲線特征值。

圖6 變異函數(shù)模型Fig.6 Variogram function model

表2 球狀模型變異函數(shù)特征值Table 2 Variogram feature values of spherical models

根據(jù)普通克里金估計(jì)的無偏和方差最小條件，可求出普通克里金方程組的加權(quán)系數(shù)[19]，進(jìn)而得到普通克里金的估計(jì)值，即測場每一個(gè)數(shù)據(jù)點(diǎn)的氡濃度計(jì)算值(見圖7).

圖7 普通克里金插值分布Fig.7 Ordinary Kriging interpolation distribution

3 結(jié)果及討論

將測場區(qū)域內(nèi)西南角第一個(gè)測點(diǎn)作為坐標(biāo)(0,0)點(diǎn)，以測場東西方向的邊界連線作為x軸，以南北方向連線作為y軸，構(gòu)建測場坐標(biāo)系，如圖7、圖8中內(nèi)圈數(shù)據(jù)所示。兩圖中外圈數(shù)據(jù)代表通用橫墨卡托投影坐標(biāo)(universal transverse mercartor，UTM)，是由各測點(diǎn)經(jīng)緯度坐標(biāo)經(jīng)投影變換所形成(通過GPS定位的測點(diǎn)經(jīng)緯度坐標(biāo)屬于地理坐標(biāo)系統(tǒng)，無法構(gòu)建平面插值分布圖，需轉(zhuǎn)換為投影坐標(biāo)系)。

圖8 最小曲率插值分布Fig.8 Minimum curvature interpolation distribution

3.1 普通克里金插值分析

結(jié)合圖4中箱形圖分析所給出的異常值范圍(487 N/3 min～743 N/3 min)，在圖7所示區(qū)域內(nèi)圈定氡濃度異常區(qū)域2處：區(qū)域A形狀近似橢圓狀，呈西北-東南走向，面積約800 m2，域內(nèi)氡濃度極值點(diǎn)有2處，坐標(biāo)分別為A1(86，247)、A2(106，228)，氡計(jì)數(shù)為685 N/3 min、637 N/3 min.區(qū)域B呈東西走向，面積約850 m2，氡濃度極值點(diǎn)有2處，坐標(biāo)B1(83，42)、B2(112，46)，組成“雙峰”向四周逐漸遞減，氡計(jì)數(shù)分別為743 N/3 min、658 N/3 min.以上4點(diǎn)均呈現(xiàn)由中心極值向四周均勻遞減趨勢，基本符合煤自燃演化規(guī)律，可以確定為疑似著火點(diǎn)。需要特別說明的是，從4處疑似火源點(diǎn)的橫縱坐標(biāo)可以看出其并不局限于測點(diǎn)坐標(biāo)，這是因?yàn)橥ㄟ^插值在測點(diǎn)之間形成了若干數(shù)據(jù)點(diǎn)，均可作為火源點(diǎn)的位置。

3.2 最小曲率法插值結(jié)果分析

同樣結(jié)合箱形圖分析所給出的異常值范圍，在圖8所示最小曲率法插值分布圖中圈定4處疑似著火區(qū)域，較之普通克里金插值結(jié)果增加C、D兩處。區(qū)域C形狀近似圓形，面積約450 m2，域內(nèi)極值點(diǎn)1處，坐標(biāo)為C1(108,118)，氡計(jì)數(shù)為624 N/3 min.區(qū)域D呈東北-西南走向，面積約700 m2，域內(nèi)極值點(diǎn)1處，坐標(biāo)為D1(258,247)，氡計(jì)數(shù)為547 N/3 min，可確定C1、D1為疑似著火點(diǎn)。

3.3 鉆孔驗(yàn)證

分別對以上疑似著火點(diǎn)位置進(jìn)行鉆孔測溫驗(yàn)證，并對孔底CO濃度進(jìn)行測定，最終結(jié)果如表3所示。

表3 孔底數(shù)據(jù)表Table 3 Hole bottom data sheet

點(diǎn)A1、A2、B1、B2均超過了煤炭自燃的臨界溫度(60 ℃～85 ℃)且存在CO氣體，并且在鉆探過程中點(diǎn)A1、B1、B2孔口有劇烈的熱氣冒出，所取煤心有明顯燃燒過的特征；點(diǎn)A2在鉆探過程中孔口冒熱氣，所取煤心有烘烤現(xiàn)象，表明以上疑似著火點(diǎn)確定存在不同程度的煤自燃現(xiàn)象。點(diǎn)C1孔底溫度介于自燃臨界溫度之間，鉆探過程中孔口有少量熱氣冒出，所取煤心無燃燒或烘烤跡象，說明該點(diǎn)尚處于自熱階段。點(diǎn)D1孔底溫度較低，煤層無烘烤現(xiàn)象，頂板砂巖新鮮，且無CO氣體存在，可以認(rèn)為該鉆孔區(qū)域未發(fā)生煤炭自燃。

綜上所述，根據(jù)普通克里金法插值形成的氡濃度分布所圈定的火源位置，排除了最小曲率法所圈定的疑似著火點(diǎn)C1、D1，準(zhǔn)確度更高。

4 結(jié)論

1) 反應(yīng)氡濃度變異函數(shù)的數(shù)球狀模型特征值中，塊金值與基臺值的比值較小(0.340 7)，介于0.25～0.75之間，說明土壤氡這種變量具有中等的空間相關(guān)性，即一定范圍內(nèi)的相鄰測點(diǎn)之間能夠相互影響，也從側(cè)面證明了使用地統(tǒng)計(jì)學(xué)中克里金插值方法的必要性。

2) 使用普通克里金插值方法能夠?qū)⒆钚∏史ㄈΧǖ?處疑似火區(qū)(A、B、C、D)縮小為2處(A、B)，將疑似高溫火源點(diǎn)由6個(gè)(A1、A2、B1、B2、C1、D1)減小到4個(gè)(A1、A2、B1、B2)，鉆孔驗(yàn)證的高溫火源點(diǎn)位置與普通克里金插值方法所圈定的結(jié)果一致，說明使用普通克里金插值法在針對土壤氡濃度的區(qū)域化變量的分析中具有一定的優(yōu)越性，能夠避免最小曲率法造成的插值失真，可提高同位素測氡法對火源位置判斷的準(zhǔn)確度，對煤自燃火災(zāi)的防治具有指導(dǎo)作用。