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

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

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

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

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

1 目標(biāo)與方法

1.1 變異函數(shù)

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

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

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

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

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

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

1.2 研究方法

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

1.2.1 交叉驗(yàn)證

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

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

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

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

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

1.2.3 控制變量法

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

1.2.4 Sobol′方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

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

3.1.1 組樣長(zhǎng)度對(duì)實(shí)驗(yàn)變異函數(shù)的影響

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

圖4 樣品采樣長(zhǎng)度分布直方圖Fig.4 Histogram of sample length distribution

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

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

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

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

3.1.2 數(shù)據(jù)分布對(duì)實(shí)驗(yàn)變異函數(shù)的影響

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

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

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

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

3.1.3 滯后距對(duì)實(shí)驗(yàn)變異函數(shù)的影響

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

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

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

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

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

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

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

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

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

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

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

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

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

4 討 論

4.1 實(shí)驗(yàn)變異函數(shù)影響因素

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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