三維地質(zhì)統(tǒng)計學(xué)礦產(chǎn)資源量估算變異函數(shù)影響因素與參數(shù)敏感性研究

李曉暉 江健華 馬 良 祁軼宏 潘 楊

(1.合肥工業(yè)大學(xué)資源與環(huán)境工程學(xué)院,安徽 合肥 230009;2.安徽省公益性地質(zhì)調(diào)查管理中心(安徽省地質(zhì)調(diào)查與環(huán)境監(jiān)測中心),安徽 合肥 230051;3.安徽省礦產(chǎn)資源儲量評審中心,安徽 合肥 230051)

礦產(chǎn)資源量估算對于評估礦石品位分布、礦床資源總量、礦山邊界與壽命等具有十分重要的作用[1]。隨著計算機技術(shù)與數(shù)字地質(zhì)的飛速發(fā)展,基于三維地質(zhì)信息化技術(shù)的三維地質(zhì)統(tǒng)計學(xué)方法在礦產(chǎn)資源量估算領(lǐng)域的應(yīng)用越發(fā)廣泛,國外礦業(yè)大國已將該方法作為礦產(chǎn)資源量估算的重要手段[2]。隨著國際合作不斷加深,近年來三維地質(zhì)統(tǒng)計學(xué)礦產(chǎn)資源量估算方法在國內(nèi)應(yīng)用也日益增多[3-4]。

變異函數(shù)作為地質(zhì)統(tǒng)計學(xué)的重要工具,是克里格估值的基礎(chǔ),確定合理的變異函數(shù)對于礦產(chǎn)資源量估算至關(guān)重要。確定變異函數(shù)過程主要包括構(gòu)建穩(wěn)健的實驗變異函數(shù)模型以及對理論變異函數(shù)模型進行擬合兩個環(huán)節(jié)。但由于成礦地質(zhì)條件、勘查類型、礦床模型、品位數(shù)據(jù)分布等均存在差異,實驗變異函數(shù)模型構(gòu)建與理論變異函數(shù)模型擬合在實際應(yīng)用中較難把握,通常伴隨有礦產(chǎn)資源量估算者的主觀性,且大多借助計算機進行操作,估算過程難以全面和直觀展現(xiàn),有學(xué)者稱地質(zhì)統(tǒng)計學(xué)法資源量估算為“暗箱”操作,可靠性有待進一步提升[5]。因此,業(yè)內(nèi)學(xué)者針對變異函數(shù)模型構(gòu)建中的不同方面開展了相關(guān)研究,如塊金效應(yīng)[6-7]、特高品位影響[8]、各向異性結(jié)構(gòu)套[9]、算法擬合[10-11]等。變異函數(shù)模型構(gòu)建是一個整體性工作,由多要素共同決定,目前研究主要針對變異函數(shù)中某一影響因素或參數(shù)進行探討,缺少不同要素對變異函數(shù)影響程度的對比分析,并且鮮有針對其涉及的諸多參數(shù)開展系統(tǒng)性研究。

基于安徽省安慶—貴池礦集區(qū)某大型鐵銅礦主礦體數(shù)據(jù),本研究對三維地質(zhì)統(tǒng)計學(xué)礦產(chǎn)資源量估算中實驗變異函數(shù)穩(wěn)健性的影響因素、理論變異函數(shù)模型各參數(shù)的敏感性分別開展分析,系統(tǒng)評價不同影響因素及參數(shù)變化對變異函數(shù)的影響,以期明確三維變異函數(shù)構(gòu)建過程中的風(fēng)險性要素,從而更好地推動三維地質(zhì)統(tǒng)計學(xué)礦產(chǎn)資源量估算方法的進一步實踐應(yīng)用。

1 目標(biāo)與方法

1.1 變異函數(shù)

1.1.1 實驗變異函數(shù)穩(wěn)健性影響因素

綜合前人研究,影響實驗變異函數(shù)穩(wěn)健性的因素包括組樣長度、數(shù)據(jù)分布、滯后距等[12]。上述因素會對實驗變異函數(shù)穩(wěn)健性產(chǎn)生一定的影響,例如:使實驗變異函數(shù)曲線變得復(fù)雜,或使實驗變異函數(shù)曲線丟失相關(guān)信息點,又或使實驗變異函數(shù)曲線波動增大導(dǎo)致無法擬合理論變異函數(shù)模型等[13-14]。因此,本研究對組樣長度、數(shù)據(jù)分布、滯后距的影響因素開展分析,構(gòu)建更為穩(wěn)健的實驗變異函數(shù),以便獲得更加可靠的理論變異函數(shù)模型。

1.1.2 理論變異函數(shù)擬合參數(shù)

在理論變異函數(shù)擬合過程中,需要確定的相關(guān)參數(shù)包括塊金值、基臺值、變程與3個垂直方向上的各向異性比率等參數(shù)[15]。上述參數(shù)共同決定了3個方向理論變異函數(shù)的形態(tài),雖然各參數(shù)可以通過交叉驗證等方式調(diào)整,但在一定程度上還依賴于礦產(chǎn)資源量估算者的經(jīng)驗,需要估算人員根據(jù)自身知識儲備開展具體的分析和研究。如圖1所示,相同的實驗變異函數(shù)數(shù)據(jù)在不同的滯后距條件下,可以擬合獲得不同形態(tài)的理論變異函數(shù)模型,雖然趨勢大致相同,但塊金值、基臺值等關(guān)鍵參數(shù)存在顯著差異,進而影響三維資源量估算結(jié)果。

圖1 理論變異函數(shù)擬合結(jié)果Fig.1 Fitting results of theoretical variogram

因此,本研究針對上述擬合參數(shù)開展敏感性分析,以便更好地確定理論變異函數(shù)模型,提高礦產(chǎn)資源量估算精度。

1.2 研究方法

針對實驗變異函數(shù)穩(wěn)健性的相關(guān)影響因素,本研究采用數(shù)學(xué)統(tǒng)計方法中方差、標(biāo)準(zhǔn)差以及變異系數(shù)等統(tǒng)計量開展分析,定量評價實驗變異函數(shù)的穩(wěn)健性。針對理論變異函數(shù)模型的諸多擬合參數(shù)的敏感性,基于交叉驗證結(jié)果,本研究采用Sobol′敏感性分析方法與控制變量法開展分析。

1.2.1 交叉驗證

在構(gòu)建更為合理的理論變異函數(shù)模型時,通常會采用交叉驗證方法開展分析論證。該方法原理是將每個實測點抽出,用其周圍點對該點進行克里格估值,然后用實測值與克里格估值結(jié)果進行比較,最終評判理論變異函數(shù)擬合的優(yōu)劣[16]。理論上,合理的變異函數(shù)模型其交叉驗證指標(biāo)應(yīng)滿足的條件為[17]:① 實際克里格誤差均值(ME)趨近于0;② 實際克里格誤差方差與理論克里格方差趨于相等,兩者相對誤差(RE)在15%以內(nèi);③ 2倍標(biāo)準(zhǔn)差范圍內(nèi)誤差(TSD)所占比例大于95%。

1.2.2 數(shù)學(xué)統(tǒng)計方法

數(shù)學(xué)統(tǒng)計方法中方差、標(biāo)準(zhǔn)差、變異系數(shù)等統(tǒng)計量可用于反映數(shù)據(jù)的離散程度、變異程度,以及定量分析實驗變異函數(shù)的穩(wěn)健性。方差是單個樣本值與全體樣本值平均值的離差平方的平均數(shù);標(biāo)準(zhǔn)差為方差的算術(shù)平方根,用S表示。在實際應(yīng)用中,方差計算公式為

變異系數(shù)可用于比較量綱不同的數(shù)據(jù)的離散程度,計算公式為

式中,為全體樣本的平均值,n為樣本數(shù)。

1.2.3 控制變量法

在物理學(xué)中,常采用控制變量法探究多因素或多變量的相關(guān)問題?？刂谱兞糠看沃桓淖兤渲心骋粋€變量,而保持其他變量不變,從而分析被控制的變量對輸出結(jié)果的影響[18]。理論變異函數(shù)模型是一個含多參數(shù)的模型,通過控制變量法可以有效地探究單一參數(shù)值變化對交叉驗證結(jié)果的影響,厘定理論變異函數(shù)模型中的風(fēng)險性參數(shù)。

1.2.4 Sobol′方法

Sobol′方法是由俄羅斯數(shù)學(xué)家I.M.Sobol提出的一種基于方差方法的全局敏感性分析方法,完全基于模型運行后得到的輸出值進行計算[19]。該方法假定輸入的參數(shù)能夠獨立且均勻地在d維單位超立方體Ωd中取值,f(x)為平方可積函數(shù),即Ωd=,d為參數(shù)個數(shù)。那么,根據(jù)HDMR(高維模型擬合)方法,f(x)可進行如下表示

對式(3)進行平方積分運算,可得到方差V(Y)的分解形式:

式中,Vi為參數(shù)Xi產(chǎn)生的方差;Vij為參數(shù)Xi和Xj共同作用產(chǎn)生的方差;V1,2,…,d為所有參數(shù)共同作用產(chǎn)生的方差。

對式(4)兩邊同時除以V(Y),并定義:

式中,Si為一階敏感度系數(shù),被稱為“重要性指標(biāo)”,可以用來衡量參數(shù)Xi對輸出結(jié)果的貢獻大小,通?？梢杂糜谧R別風(fēng)險性參數(shù),Si越大表明參數(shù)Xi在合理范圍內(nèi)變化會對輸出結(jié)果產(chǎn)生較大影響;Sij,…,S1,2,…,d為與參數(shù)Xi有關(guān)的高階敏感性指數(shù),表示多個參數(shù)組合相互作用時的敏感度。

為了衡量參數(shù)Xi對輸出結(jié)果所能產(chǎn)生的全部影響,定義總階敏感性系數(shù)STi為

式中,-i表示除Xi外所有參數(shù),如果STi趨近于0,則表明參數(shù)Xi在合理范圍內(nèi)會隨意變動,而輸出結(jié)果不發(fā)生明顯改變?？傠A敏感性系數(shù)用來表示參數(shù)Xi對輸出結(jié)果方差的總貢獻,包括Xi的一階指數(shù)加上所有由于相互作用而產(chǎn)生的高階指數(shù)。

因此,總階敏感性系數(shù)與一階敏感度系數(shù)的差,可以用來分析第i個參數(shù)與其他所有參數(shù)的相互作用影響,當(dāng)差值越大,表示該參數(shù)與其他參數(shù)之間的相互作用越強;反之,亦然。

已有研究表明,Sobol′方法能夠更全面地反映多參數(shù)共同變化時模型輸出結(jié)果的響應(yīng)[20]。因此,該方法非常適用于分析多參數(shù)共同作用的理論變異函數(shù)模型,通過定性定量分析能夠較好地識別出敏感性參數(shù),輔助理論變異函數(shù)模型構(gòu)建,從而提高礦產(chǎn)資源量估算精度。

運用Sobol′方法對理論變異函數(shù)模型參數(shù)開展敏感性分析,步驟為:① 確定各參數(shù)合理取值范圍,并據(jù)此進行Monte Carlo采樣,得到多組用于擬合理論變異函數(shù)模型的參數(shù)樣本;② 對采樣得到的參數(shù)樣本進行交叉驗證,得到交叉驗證評價指標(biāo),用于評價理論變異函數(shù)擬合效果;③ 根據(jù)采樣得到的參數(shù)樣本及其對應(yīng)的交叉驗證結(jié)果,運用Sobol′方法計算理論變異函數(shù)參數(shù)的一階敏感度系數(shù)和總階敏感度系數(shù),以定量分析其參數(shù)的敏感性,具體流程如圖2所示[21]。

圖2 Sobol′方法敏感性分析流程Fig.2 Sensitivity analysis process of Sobol′ method

2 實例數(shù)據(jù)

本研究以安徽省安慶—貴池礦集區(qū)某大型鐵銅礦床控制的主礦體為例進行分析,主礦體數(shù)據(jù)來源于該礦床經(jīng)評審備案的勘探報告中的相關(guān)成果資料。該礦床為典型的矽卡巖型鐵銅礦床,礦體主體賦存于閃長巖體與三疊系中統(tǒng)月山組灰?guī)r的接觸帶附近。該礦床勘查類型為Ⅱ類,矽卡巖型鐵礦勘查基本網(wǎng)度為200 m×200 m,并加密至100 m×100 m網(wǎng)度。控制礦床主礦體的鉆孔數(shù)量共計32個,主礦體呈似層狀、透鏡狀產(chǎn)出,礦體產(chǎn)狀平緩,無較大的巖脈破壞切割礦體,也未見斷裂構(gòu)造切割破壞礦體。主礦體近SN走向,厚度變化趨勢為中部厚向南北兩端減薄。

本研究選擇TFe品位作為區(qū)域化變量,以全鐵資源量作為估算目標(biāo)開展三維地質(zhì)統(tǒng)計學(xué)資源量估算中變異函數(shù)的相關(guān)研究。TFe品位最小值為3.1%,最大值為67.03%,無特高品位,變化系數(shù)為23%,變化系數(shù)小于 50%,表明礦化均勻連續(xù),主礦體樣品數(shù)據(jù)統(tǒng)計分析結(jié)果見表1。

表1 礦體樣品統(tǒng)計結(jié)果Table 1 Statistical results of orebody samples

基于礦床勘探數(shù)據(jù),本研究首先在Surpac三維地質(zhì)軟件中建立了該礦床地質(zhì)鉆孔數(shù)據(jù)庫;然后依據(jù)《礦產(chǎn)地質(zhì)勘查規(guī)范 鐵、錳、鉻》(DZ/T 0200—2020)對礦體進行解譯;最后采用人機交互方式,基于Surpac的三角網(wǎng)化工具構(gòu)建了該礦體的三維實體模型(圖3)。該模型將作為估算域參與三維地質(zhì)統(tǒng)計學(xué)資源量估算研究。

3 結(jié)果與分析

3.1 實驗變異函數(shù)穩(wěn)健性影響因素分析

3.1.1 組樣長度對實驗變異函數(shù)的影響

主礦體內(nèi)所有樣品采樣長度分布特征如圖4所示。

圖4 樣品采樣長度分布直方圖Fig.4 Histogram of sample length distribution

為分析組樣長度對實驗變異函數(shù)穩(wěn)健性的影響,對主礦體樣品分別按1.0、1.2、1.4、1.6、1.8、2.0 m長度進行組合,并計算全方位變異函數(shù)。為了消除滯后距這一影響因素,分別針對6組不同組樣長度下組合成的數(shù)據(jù),計算在60、70、80 m 3個不同滯后距下的全方位變異函數(shù),結(jié)果如圖5所示。

圖5 不同組樣長度下全方位變異函數(shù)Fig.5 Omnidirectional variograms with different sample composited lengths

通過對比實驗變異函數(shù)曲線可知:在不同滯后距下,6組實驗變異函數(shù)曲線的變異性差距均較小,變化趨勢相同。雷達圖是一種表現(xiàn)多維數(shù)據(jù)的圖表,其各軸代表不同的組樣長度,軸坐標(biāo)代表變異系數(shù),連接各坐標(biāo)軸上的數(shù)據(jù)點可以形成一個多邊形域,可通過多邊形域形狀直觀了解實驗變異函數(shù)曲線的變異程度。本研究計算上述實驗變異函數(shù)曲線的變異系數(shù),并通過“雷達圖”展示(圖6)。由圖6可知:3個多邊形域的形狀基本均構(gòu)成“正六邊形”,說明在不同滯后距下,基于不同組樣長度的實驗變異函數(shù)曲線其變異系數(shù)趨于相同。由此表明:組樣長度在一定區(qū)間內(nèi)對于實驗變異函數(shù)的穩(wěn)健性不會產(chǎn)生實質(zhì)變化,故取均值1.50 m作為組樣長度對樣品進行組合,進行后續(xù)分析。

圖6 不同組樣長度下實驗變異函數(shù)曲線變異系數(shù)雷達圖Fig.6 Radar chart of variation coefficient of experimental variogram curve under different sample composited lengths

3.1.2 數(shù)據(jù)分布對實驗變異函數(shù)的影響

對主礦體組合樣數(shù)據(jù)進行了統(tǒng)計分析,通過Q-Q圖(圖7(a))可以看出,原數(shù)據(jù)呈現(xiàn)明顯的負(fù)偏態(tài)分布。對組合樣數(shù)據(jù)進行對數(shù)轉(zhuǎn)換,轉(zhuǎn)換后的數(shù)據(jù)基本呈現(xiàn)正態(tài)分布(圖7(b))。

圖7 組合樣對數(shù)轉(zhuǎn)換前后Q-Q圖Fig.7 Q-Q diagram before and after logarithmic transformation of composite samples

為探究數(shù)據(jù)分布對實驗變異函數(shù)的影響,本研究對主礦體未經(jīng)對數(shù)轉(zhuǎn)換與經(jīng)過對數(shù)轉(zhuǎn)換的組合樣數(shù)據(jù)構(gòu)建3個方向?qū)嶒炞儺惡瘮?shù),結(jié)果如圖8所示。由圖8可知:數(shù)據(jù)分布會對實驗變異函數(shù)穩(wěn)健性產(chǎn)生一定的影響,未經(jīng)對數(shù)變換之前,3個方向?qū)嶒炞儺惡瘮?shù)曲線差異較大,且波動幅度較大,沒有一個明顯趨于平緩的趨勢,基臺值與變程確定困難;經(jīng)過對數(shù)變換后,隨著距離增加,基臺值逐漸趨近到一個較為穩(wěn)定的數(shù)值,3個方向?qū)嶒炞儺惡瘮?shù)曲線的穩(wěn)定性得到了很大提高,為后期理論變異函數(shù)擬合提供了較好條件。

圖8 組合樣對數(shù)轉(zhuǎn)換前后主礦體三方向?qū)嶒炞儺惡瘮?shù)Fig.8 Three direction experiment variogram of main orebody before and after logarithmic conversion of composite samples

3.1.3 滯后距對實驗變異函數(shù)的影響

為分析滯后距對實驗變異函數(shù)穩(wěn)健性的影響,本研究以10 m為初始滯后距,并以10 m為增量,計算不同滯后距條件下礦體全方位的實驗變異函數(shù),共計獲得12種不同滯后距下的實驗變異函數(shù),如圖9所示。由圖9可知:隨著滯后距增大,實驗變異函數(shù)的穩(wěn)健性逐漸提高,但實驗變異函數(shù)曲線包含的結(jié)構(gòu)點減少,因此,對于滯后距的選取應(yīng)充分考察變異函數(shù)的特征結(jié)構(gòu),進行綜合研判。

圖9 不同滯后距下全方位實驗變異函數(shù)Fig.9 Variogram of omnidirectional experiment under different lag distance

為更好地選取合理的滯后距,以構(gòu)建穩(wěn)健且有效的實驗變異函數(shù)模型,本研究進一步對實驗變異函數(shù)值的方差、標(biāo)準(zhǔn)差、變異系數(shù)等統(tǒng)計量進行了分析,結(jié)果如圖10所示。由圖10可知:隨著滯后距增加,實驗變異函數(shù)值的方差、標(biāo)準(zhǔn)差和變異系數(shù)總體上呈現(xiàn)下降趨勢。滯后距為10～60 m時,方差、標(biāo)準(zhǔn)差和變異系數(shù)迅速下降,曲線穩(wěn)定性有所上升;滯后距為60～80 m時,實驗變異函數(shù)曲線較為平穩(wěn);滯后距為80～120 m時,方差、標(biāo)準(zhǔn)差和變異雖有小幅波動,但總體已相對穩(wěn)定?？紤]該礦區(qū)勘查網(wǎng)度為100 m×100 m,部分鉆孔間距約為80 m,當(dāng)滯后距超過80 m時,可能會導(dǎo)致平滑效應(yīng)出現(xiàn),消弭了礦體局部變異特征。因此,本研究結(jié)合勘查網(wǎng)度,最終選取在60～80 m滯后距范圍構(gòu)建實驗變異函數(shù),以開展進一步的分析討論。

圖10 方差、標(biāo)準(zhǔn)差、變異系數(shù)變化曲線Fig.10 Variation curves of variance,standard deviation and coefficient of variation

3.2 理論變異函數(shù)模型參數(shù)敏感性分析

本研究基于Sobol′方法對理論變異函數(shù)模型參數(shù)開展敏感性分析。結(jié)合礦床地質(zhì)信息,依據(jù)上述構(gòu)建的實驗變異函數(shù)形態(tài),采用球狀模型對上述構(gòu)建的實驗變異函數(shù)模型進行擬合。通過人機交互方式,確定理論變異函數(shù)模型參數(shù)合理擬合范圍(表2)。本研究進一步在合理范圍內(nèi)利用Monte Carlo方法對各參數(shù)進行隨機模擬采樣,采樣得到325組理論變異函數(shù)擬合參數(shù),并對其進行交叉驗證,評價理論變異函數(shù)參數(shù)取值優(yōu)劣。

表2 理論變異函數(shù)模型參數(shù)合理擬合范圍Table 2 Reasonable fitting range of theoretical variogram model parameters

本研究計算了塊金值、基臺值、變程、主軸/次主軸比率、主軸/次軸比率等5個理論變異函數(shù)參數(shù)對ME、RE、TSD3個交叉驗證評價標(biāo)準(zhǔn)的一階敏感度系數(shù)(S1)與總階敏感度系數(shù)(ST),結(jié)果如圖11所示。圖11中填充部分代表一階敏感度(S1),白色部分代表一階和全階敏感度的差值(ST-S1)。由圖11可知:塊金值對于ME、RE、TSD的一階敏感度系數(shù)均超過0.8,表現(xiàn)了強烈的顯著敏感;變程對于RE、TSD的一階敏感度系數(shù)接近于0.1,在5個參數(shù)中為第2敏感參數(shù);基臺值、主軸/次主軸比率和主軸/次軸比率對于ME、RE、TSD的一階敏感度系數(shù)趨于0,表現(xiàn)不敏感。

圖11 交叉驗證評價指標(biāo)的一階敏感度系數(shù)與總階敏感度系數(shù)Fig.11 First order sensitivity coefficient and total order sensitivity coefficient of cross-validation evaluation index

為研究各參數(shù)變化對交叉驗證評價指標(biāo)的具體影響,采用控制變量法對單一參數(shù)進行敏感性分析,進一步探討單一參數(shù)值的變化與交叉驗證評價指標(biāo)是否存在線性關(guān)系。為不失一般性,本研究不強調(diào)參數(shù)具體數(shù)值大小,將合理擬合范圍的五分位點作為各參數(shù)的變化值,并將合理擬合范圍的中值確定為參數(shù)的固定值(表3)。基于表3,可組成25組理論變異函數(shù)擬合參數(shù)樣本,利用擬合參數(shù)樣本進行普通克里格估值,并進行交叉驗證,其交叉驗證結(jié)果中的ME、RE、TSD的變化趨勢如圖12所示。

表3 單一參數(shù)敏感性分析參數(shù)擬合值Table 3 Parameter fitting value of single parameter sensitivity analysis

圖12 單一參數(shù)敏感性分析交叉驗證結(jié)果Fig.12 Cross-validation results of single parameter sensitivity analysis

整體上來看,塊金值對ME、RE、TSD的影響最大,變程次之,其余3個參數(shù)的影響較小,與Sobol′方法敏感度分析結(jié)果大致相同。同時,從圖12中可以發(fā)現(xiàn):各參數(shù)值的變化與交叉驗證結(jié)果評價指標(biāo)的優(yōu)劣呈一定的線性相關(guān)。隨著塊金值不斷提高,ME雖有所擴大但仍保持一個近于0的水平,絕對值小于0.005;而在塊金值不斷提高后,RE與TSD都得到明顯改善,在塊金值達到0.075后,RE與TSD均基本滿足交叉驗證結(jié)果的評價標(biāo)準(zhǔn)。對于變程的擴大,ME沒有明顯變化,其絕對值維持在0.003左右,但RE略有提高且TSD略微降低。

4 討 論

4.1 實驗變異函數(shù)影響因素

通過實驗變異函數(shù)穩(wěn)健性影響因素分析發(fā)現(xiàn),組樣長度對于實驗變異函數(shù)穩(wěn)健性影響較小,而數(shù)據(jù)分布、滯后距等因素會對實驗變異函數(shù)穩(wěn)健性產(chǎn)生較大影響。

數(shù)據(jù)分布方面,由圖8可知:經(jīng)過對數(shù)變化后數(shù)據(jù)計算的實驗變異函數(shù)曲線穩(wěn)健性有所提高,主要是由于對數(shù)變換使得數(shù)據(jù)方差降低,減小了偏態(tài)數(shù)據(jù)對實驗變異函數(shù)曲線穩(wěn)健性的影響。目前,多數(shù)研究認(rèn)為當(dāng)樣品品位呈偏態(tài)分布或非正態(tài)分布時,應(yīng)采用的地質(zhì)統(tǒng)計學(xué)資源量估算方法主要有對數(shù)克里格法或指示克里格法[22],因為這些方法其核心均是通過一定方式降低或避免非正態(tài)數(shù)據(jù)方差過大對實驗變異函數(shù)曲線帶來的影響。

關(guān)于滯后距,圖9和圖10表明:實驗變異函數(shù)的穩(wěn)健性會隨滯后距的擴大而提高,但當(dāng)滯后距大于某一值時,實驗變異函數(shù)曲線將不再發(fā)生明顯變化。在過大的滯后距下,雖然實驗變異函數(shù)曲線趨于平穩(wěn),但其結(jié)構(gòu)點不斷減少,導(dǎo)致平滑效應(yīng)出現(xiàn),丟失了部分地質(zhì)信息。本研究通過對實驗變異函數(shù)值的方差、標(biāo)準(zhǔn)差、變異系數(shù)等統(tǒng)計量進行分析,結(jié)合礦區(qū)勘查網(wǎng)度,確定了合理的滯后距范圍。該處理方式既避免了過大滯后距消弭礦體局部特征,又充分保證了實驗變異函數(shù)的穩(wěn)健性,可以有效指導(dǎo)實驗變異函數(shù)滯后距選取,將其應(yīng)用于進一步的理論變異函數(shù)擬合工作。

4.2 理論變異函數(shù)擬合參數(shù)

通過Sobol′方法敏感性分析(圖11),反映出塊金值對交叉驗證結(jié)果表現(xiàn)出了強烈的顯著敏感性;同時,就塊金值的敏感度而言,一階敏感度占比高,一階和全階敏感度的差值較小,說明塊金值與其他所有參數(shù)之間的相互作用并不強烈,塊金值在理論變異函數(shù)模型中單獨作用,因此進一步可以看出塊金值是影響估值結(jié)果最為重要的參數(shù)之一。在地質(zhì)統(tǒng)計學(xué)估值過程中,隨著塊金值不斷增大,待估塊段較近的樣品點對較遠樣品點的屏蔽效應(yīng)減弱[7],進而導(dǎo)致較遠的樣品點也參與到估值中,使得參與估值的樣品點權(quán)重趨于相同,因此會對估算結(jié)果產(chǎn)生較大影響。

在剩余參數(shù)中,變程對RE和TSD而言表現(xiàn)較為敏感,變程不斷擴大對RE與TSD指標(biāo)有不利作用(圖11、圖12)。而基臺值、主軸/次主軸比率和主軸/次軸比率等參數(shù)表現(xiàn)不敏感,在合理范圍內(nèi)會隨意變動,對于輸出結(jié)果不會發(fā)生明顯作用(圖11、圖12)。主要原因是:變程表示數(shù)據(jù)間具有空間相關(guān)性的范圍,選取較大的變程將會使得與待估點無空間相關(guān)性的樣品點用于估值計算,導(dǎo)致估算精度降低。

4.3 理論變異函數(shù)模型最優(yōu)性檢驗

在擬合理論變異函數(shù)模型時,通常會采用交叉驗證方法對模型進行檢驗,開展擬合參數(shù)的分析論證[16]。單一參數(shù)敏感性分析結(jié)果也表明了各參數(shù)值的變化與交叉驗證結(jié)果評價指標(biāo)的優(yōu)劣呈一定的線性相關(guān)。但需要注意,隨著塊金值不斷增大,交叉驗證評價指標(biāo)ME隨之?dāng)U大但同時RE不斷減小,不同指標(biāo)間出現(xiàn)背離,對于理論變異函數(shù)模型產(chǎn)生不同的評價結(jié)果(圖12)。

基于本研究對于塊金值的探討,認(rèn)為ME增大是由于塊金值增大導(dǎo)致平滑效應(yīng)增強,估值結(jié)果趨于均勻化且與真實值偏離;RE減小是由于塊金值增大使得用于估計的已知數(shù)據(jù)點影響范圍擴大,屏蔽效應(yīng)減弱,從而減小了實際克里格誤差方差與理論克里格方差,使其逐漸趨于相等,但也使得理論變異函數(shù)模型反映的礦體空間結(jié)構(gòu)失去真實性。因此,在本質(zhì)上這些指標(biāo)無法作為衡量理論變異函數(shù)參數(shù)優(yōu)劣的絕對標(biāo)準(zhǔn),僅通過交叉驗證結(jié)果評價指標(biāo)的變化調(diào)整參數(shù)值,可能會導(dǎo)致變異函數(shù)無法反映礦體的空間變異性,確定參數(shù)值還需綜合多方面信息。

5 結(jié) 論

(1)針對三維地質(zhì)統(tǒng)計學(xué)礦產(chǎn)資源量估算中實驗變異函數(shù)的影響因素與理論變異函數(shù)參數(shù)的敏感性開展了系統(tǒng)性研究,解析了變異函數(shù)構(gòu)建過程中的風(fēng)險性要素。研究表明:數(shù)據(jù)分布、滯后距選取等因素會對實驗變異函數(shù)穩(wěn)健性產(chǎn)生較大影響,組樣長度對于實驗變異函數(shù)穩(wěn)健性影響較小;理論變異函數(shù)模型中的塊金值對于交叉驗證結(jié)果影響最為顯著,變程次之,其余參數(shù)表現(xiàn)不敏感,塊金值和變程應(yīng)視為風(fēng)險性參數(shù),在資源量估算過程中應(yīng)進行深入研究;理論變異函數(shù)模型參數(shù)值的變化雖然與交叉驗證結(jié)果指標(biāo)的優(yōu)劣呈一定的線性相關(guān)性,但交叉驗證結(jié)果并不能作為調(diào)整參數(shù)值的絕對標(biāo)準(zhǔn),因此需要結(jié)合礦體的空間變異性,建立更為合理的評價方法與指標(biāo)。

(2)由于礦床類型繁多、地質(zhì)特征和成礦條件各異,未來可基于本研究方法針對不同種類礦床與數(shù)據(jù)開展更為綜合與全面的研究,獲取更為合理的變異函數(shù)模型和更為準(zhǔn)確的資源量估算結(jié)果,進一步推進三維地質(zhì)統(tǒng)計學(xué)礦產(chǎn)資源量估算方法的實踐與應(yīng)用。