秦飛龍,劉 劍,顏文勇,周仲禮,龔 灝,郭 科

(1.成都工業(yè)學(xué)院 大數(shù)據(jù)與人工智能學(xué)院，成都 611730；2.成都工業(yè)學(xué)院 汽車與交通學(xué)院，成都 611730；3.成都理工大學(xué) 數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室，成都 610059)

地球化學(xué)元素含量異常指的是在某個(gè)區(qū)域或空間內(nèi)，地球化學(xué)元素含量的空間分布與正常的地球化學(xué)含量分布模式產(chǎn)生偏離[1]。地質(zhì)環(huán)境、構(gòu)造背景、巖石類型等決定了地球化學(xué)元素異常的性質(zhì)，厘清地球化學(xué)異常分布規(guī)律及其何種元素異常形成礦床，就能夠有效進(jìn)行深部礦產(chǎn)預(yù)測(cè)[2]。選取地球化學(xué)異常提取方法需要充分考慮異常分布特性，才能符合實(shí)際地質(zhì)工作。因此，選取合理的方法提取地球化學(xué)異常至關(guān)重要。但由于地球化學(xué)異常形成的環(huán)境十分復(fù)雜，迄今并未有適用于各種地質(zhì)環(huán)境下的地球化學(xué)元素含量異常求取方法[3-4]。如傳統(tǒng)的三倍均方差法、概率格紙法、累計(jì)頻率方法等在異常提取時(shí)要求樣本數(shù)據(jù)滿足正態(tài)分布；然而采集到的地球化學(xué)元素?cái)?shù)據(jù)由于受到各種環(huán)境干擾，其分布規(guī)律并不呈正態(tài)分布[3]。分形方法近年來(lái)被廣泛應(yīng)用于地球化學(xué)異常提取，也有諸多成功的應(yīng)用實(shí)例；但該方法需要建立在礦化信息基礎(chǔ)上[5]。地球化學(xué)異常在地質(zhì)空間中是后尾分布[6]，而極值理論是對(duì)事件的超常大(或超常小)的隨機(jī)性狀進(jìn)行研究，估計(jì)在已有觀測(cè)水平上更為極端事件的概率[7]，趙鵬大院士曾將其描述為極值[8]。因此，地球化學(xué)異常的分布特性符合 極值理論(EVT)研究范疇，并且 EVT 中的廣義帕累托分布(GPD)正是對(duì)數(shù)據(jù)的后尾性特征進(jìn)行刻畫[9]。本文探討利用極值理論的 GPD 研究地球化學(xué)元素含量異?？臻g分布規(guī)律，進(jìn)行礦產(chǎn)預(yù)測(cè)，為深部礦產(chǎn)預(yù)測(cè)提供方法支撐。

1 地質(zhì)背景

本文的研究區(qū)為湖北大冶雞冠咀銅金礦床，礦床面積約為 2.18 km2，位于長(zhǎng)江中部和中下游的Cu、Au、Fe成礦帶以西，坐標(biāo)是114°54′42″～114°55′45″E，30°04′45″～30°05′50″N，為第四系隱伏礦床[1]。截止到 2019 年 10 月，在研究區(qū)內(nèi)共勘探到Ⅰ、Ⅱ、Ⅲ、Ⅶ主礦體群。礦體長(zhǎng)度約為 950 m，寬度為160～800 m，深度為 5～1 412 m；在水平方向上，礦體投影面積約為 0.58 km2，傾向?yàn)楸睎| 15°～72°，某些局部為北西西向[1]。Ⅰ、Ⅱ、Ⅲ、Ⅶ主礦體群為雁列狀排列，傾向北西，局部?jī)A向南[2]。主礦體主要分布在大理巖和侵入巖之間的斷裂接觸部位，并且還富集有大量的石榴石矽卡巖，在斷裂帶內(nèi)閃長(zhǎng)巖、角礫巖發(fā)育[1]。F1和 F3兩條主斷裂將研究區(qū)分為4個(gè)成礦部分(圖1)。Ⅶ號(hào)礦體為研究區(qū)的主礦體，在 22 號(hào)勘探線與 34 號(hào)勘探線之間，深度從 447 m 開(kāi)始，賦存十分豐富，具有巨大的找礦潛力[2]。因此，本文選?、魈?hào)礦體中 26 號(hào)勘探線 Mo、Ba、Zn 地球化學(xué)元素含量進(jìn)行異常識(shí)別研究，探索新的礦體靶位。

2 GPD模型設(shè)計(jì)

2.1 GPD理論基礎(chǔ)

GPD 分布是對(duì)大于某個(gè)觀測(cè)值的樣本數(shù)據(jù)進(jìn)行分布擬合，刻畫數(shù)據(jù)的后尾分布，也稱為閾值模型[10]。如果數(shù)據(jù)x1,x2,…,xn是服從同一分布的獨(dú)立隨機(jī)變量序列，在樣本數(shù)據(jù)中存在足夠大的正數(shù)數(shù)據(jù)μ，則超過(guò)μ的分布定義為

(1)

式中：μ為樣本數(shù)據(jù)的閾值；ξ為模型的形狀參數(shù)；β為模型的尺度參數(shù)。如果ξ<0 時(shí)，則μ≤x≤μ-β/ξ；ξ≥0 時(shí)，x≥μ。在實(shí)際數(shù)據(jù)處理中，(2)式應(yīng)用更為普遍。

(2)

如果β=1，(1)和(2)式稱為 GPD 的標(biāo)準(zhǔn)形式。對(duì)(2)式求導(dǎo)，可得出GPD分布的密度函數(shù)

(3)

在(2)和(3)式中，取尺度參數(shù)β=1和形狀參數(shù)ξ=-0.5,0,0.5，得出 GPD 的標(biāo)準(zhǔn)分布函數(shù)圖(圖2-A)和概率密度函數(shù)圖(圖2-B)。由圖2-A知，當(dāng)形狀參數(shù)ξ由正值向負(fù)值變換時(shí)，GPD 分布的尾部逐漸變細(xì)。如果ξ<0，GPD主要分布在[0,-1/ξ] 內(nèi)，形狀參數(shù)ξ越小，區(qū)間越窄；如果ξ>0，GPD后尾性越明顯，形狀參數(shù)ξ越大，區(qū)間越厚。由圖2-A知，GPD 的概率密度函數(shù)嚴(yán)格單調(diào)遞減。

由式(1)和(2)知，GPD模型的確定需要估計(jì)形狀參數(shù)、尺度參數(shù)和閾值。下面對(duì)對(duì)模型的閾值和參數(shù)進(jìn)行估計(jì)。

2.1.1 閾值估計(jì)

GPD的閾值估計(jì)目前并沒(méi)有精確的方法。在眾多方法中，平均超出量函數(shù)E(μ)(MEF)是最為有效且應(yīng)用較廣的閾值估計(jì)方法[2]，它是根據(jù)閾值μ和E(μ)的變化趨勢(shì)進(jìn)行閾值估計(jì)。

如果樣本數(shù)據(jù)x1,x2,…,xn為獨(dú)立同分布的隨機(jī)序列，μ為閾值，當(dāng)xi>μ，稱xi-μ為超閾值量，如果給定樣本數(shù)據(jù)的初始閾值為μ0，那么對(duì)任意閾值μ：μ>μ0，平均超出量函數(shù)E(μ)[11]定義如下

E(μ)=E(x-μ|x>μ)

(4)

進(jìn)行積分得

(5)

如果ξ<1，積分收斂，極限為0；

如果0<ξ<1，E(μ) 為

(6)

由式(6)知，E(μ)是關(guān)于μ的線性函數(shù)。

在實(shí)際數(shù)據(jù)處理中，E(μ)常常利用如下統(tǒng)計(jì)方法進(jìn)行計(jì)算[2]

(7)

式(7)中:n為樣本的個(gè)數(shù)，N+為>μ的個(gè)數(shù)；(xi-μ)+為高于μ的第i個(gè)樣品超過(guò)閾值量，也稱損失量。

對(duì)于一組數(shù)據(jù)，如果閾值μ0已知，那么，任意的x>μ0都符合GPD分布，并且由上面理論知{(x,E(x)):μ00。從而，通過(guò)E(μ)能夠有效估計(jì)數(shù)據(jù)的閾值。

2.1.2 估計(jì)形狀參數(shù)和尺度參數(shù)

在GPD分布中，估計(jì)參數(shù)的方法較多，但是精度不高，目前應(yīng)用最廣且較為有效的估計(jì)方法包括最大似然函數(shù)法[12]和矩法估計(jì)[1]，具體如下：

a.極大似然函數(shù)法

在GPD分布中，通過(guò)對(duì)方程(2)兩邊取導(dǎo)數(shù)，得GPD的概率密度函數(shù)為

(8)

同時(shí)對(duì)(8)式兩邊取自然對(duì)數(shù)，得對(duì)數(shù)似然函數(shù)

L(ξ,β:x)=-nlnβ-

(9)

在對(duì)數(shù)似然函數(shù)(9)式中，求取關(guān)于參數(shù)ξ和β的偏導(dǎo)數(shù)，得出如下的對(duì)數(shù)似然方程

(10)

b.矩法估計(jì)

在GPD分布中，利用矩法估計(jì)也可以得出形狀參數(shù)ξ和尺度參數(shù)β，并且其估計(jì)精度優(yōu)于極大似然函數(shù)方法。因?yàn)榫胤ü烙?jì)是利用數(shù)據(jù)的樣本矩反映數(shù)據(jù)的總體矩，它不會(huì)因?yàn)镚PD形式不一樣而得到不同的參數(shù)值[1]。矩法估計(jì)表達(dá)式如下

(11)

其中：p，δ分別是ti的均值與標(biāo)準(zhǔn)差；xi為樣本數(shù)據(jù)。

2.2 地球化學(xué)異常 GPD 模型構(gòu)建

由于地球化學(xué)數(shù)據(jù)X={x1,x2,…,xn} 的異常值位于觀察數(shù)據(jù)尾部，GPD研究的是高于某足夠大的閾值數(shù)據(jù)的后尾分布規(guī)律，本文利用GPD分布建立地球化學(xué)異常提取模型。由上述的相關(guān)理論可知，通過(guò)矩法估計(jì)求出地球化學(xué)數(shù)據(jù)的GPD分布的形狀參數(shù)和尺度參數(shù)，再利用平均超出量函數(shù)MEF估計(jì)數(shù)據(jù)的閾值，從而得出地球化學(xué)數(shù)據(jù)GPD分布。具體如下：

① 獲取數(shù)據(jù)X。

② 數(shù)據(jù)后尾性檢驗(yàn)：利用QQ圖或PP圖檢驗(yàn)地球化學(xué)數(shù)據(jù)是否具有后尾特性，利用 ADF (Augmented Dickey-Fuller) 檢驗(yàn)數(shù)據(jù)的平穩(wěn)性。

③ 參數(shù)估計(jì)：利用(11)式求出 GPD 分布中的參數(shù)。

④ 估計(jì)閾值：在X中選擇不相同的閾值μ，通過(guò)式(7)得出超額均值函數(shù)E(μi)，其中μi∈μ。繪制(μi,E(μi))的散點(diǎn)圖，如果出現(xiàn)某個(gè)閾值μ0，當(dāng)μ≥μ0，E(μ)關(guān)于μ是一個(gè)近似的線性函數(shù)，那么μ0就為異常值，亦為地球化學(xué)元素含量異常下限值。

⑤ 確定異常數(shù)據(jù)GPD分布：將得出的參數(shù)和閾值代入(2)或(3)式中，得出地球化學(xué)元素含量異常的GPD分布。

⑥ 結(jié)果檢驗(yàn)：模型確定后，需要進(jìn)行診斷性檢驗(yàn)來(lái)確定求出的模型分布與實(shí)際分布是否吻合，如果效果不佳，返回第④步，估計(jì)新的閾值。診斷性檢驗(yàn)常用的方法有PP圖和QQ圖。

在以上構(gòu)建的異常提取模型中，閾值選取必須要合理，如果模型的閾值估計(jì)較高，則后尾數(shù)據(jù)xi(xi>μ;i=1,2,…,n)變少，造成后尾數(shù)據(jù)整體較大，影響參數(shù)估計(jì)的敏感性，使得方差變大；反之，閾值較小，后尾數(shù)據(jù)量增多，有利于提高GPD參數(shù)估計(jì)的精度，但存在部分?jǐn)?shù)據(jù)并不符合GPD分布。并且，后尾數(shù)據(jù)量增多，數(shù)據(jù)的中心發(fā)生變化，使得估計(jì)值與理論值具有一定的差異。

盡管MEF在估計(jì)閾值μ上具有計(jì)算快、比較直觀等優(yōu)勢(shì)，但僅僅通過(guò)線性方式進(jìn)行閾值估計(jì)存在一定的缺陷。第一，該方法是在某個(gè)小區(qū)間范圍內(nèi)選取閾值，不能確定某一具體異常下限值，會(huì)造成地球化學(xué)元素含量異常分布區(qū)域過(guò)大及過(guò)小的問(wèn)題，給實(shí)際勘探帶來(lái)困難；第二，如果出現(xiàn)2個(gè)閾值μ1和μ2，使得平均超出量函數(shù)都呈現(xiàn)線性，可以選擇線性關(guān)系明顯的地方作為閾值。但是線性不明顯，選取就存在困難。因此，本文在MEF原理基礎(chǔ)上進(jìn)行了一定的修正，得出了一種更為穩(wěn)健的修正閾值法確定閾值，具體如下：

根據(jù)廣義帕累托分布原理，只要確立了某個(gè)閾值，對(duì)于所有閾值μ(μ>μ0)，樣本數(shù)據(jù)的 GPD分布的形狀參數(shù)和尺度參數(shù)是不會(huì)改變的，并由極值類型定量[13]可得

β(μ)GPD=δGEV+ξ(μ-μ0)

(12)

因此，對(duì)于μ>μ0，有

β(μ)=δ+ξ(μ0)或β(μ)=β(μ0)+ξ(μ-μ0)

(13)

如果令

β*(μ)=β(μ)-ξ(μ)

(14)

β*(μ) 不會(huì)隨閾值μ的變換發(fā)生變化，本文將β*(μ)稱為修正的閾值參數(shù)。因此，在某個(gè)閾值區(qū)間內(nèi)，經(jīng)修正函數(shù)擬合出β*(μ)不會(huì)隨閾值μ發(fā)生變化的拐點(diǎn)μ0作為閾值，并且大于μ0的所有超出量均滿足GPD分布。

2.3 GPD模型處理地球化學(xué)數(shù)據(jù)的有效性

本文所構(gòu)建的地球化學(xué)元素含量異常提取GPD模型充分考慮了異常在地質(zhì)中的實(shí)際分布特性。地球化學(xué)異常值相對(duì)于背景值來(lái)說(shuō)，它是位于整個(gè)分布的尾部，具有后尾特性，反映了數(shù)據(jù)在超常大(或小)的隨機(jī)形狀特征。EVT是研究數(shù)據(jù)的統(tǒng)計(jì)規(guī)律，是統(tǒng)計(jì)學(xué)的一個(gè)分支，用來(lái)刻畫數(shù)據(jù)的極值分布問(wèn)題，并且EVT中的GPD正是對(duì)數(shù)據(jù)的后尾性特征進(jìn)行刻畫。因此，利用GPD分布刻畫異常值的分布規(guī)律是完全可行的，并且本文還對(duì)閾值的選取設(shè)計(jì)了修正函數(shù)，能夠確定數(shù)據(jù)的某一精準(zhǔn)異常值，并已經(jīng)進(jìn)行了充分證明[1]。下面利用實(shí)際真實(shí)數(shù)據(jù)進(jìn)行驗(yàn)證說(shuō)明。

3 GPD模型分析原生暈地球化學(xué)元素異常分布規(guī)律

3.1 數(shù)據(jù)獲取

為了體現(xiàn)本文模型的有效性，將其應(yīng)用于研究區(qū)進(jìn)行數(shù)據(jù)處理，選?、魈?hào)礦體26號(hào)勘探線上的Mo、Ba、Zn成礦元素的含量。3種元素的面含量隨著深度的變化如表1所示。顯然，深度越深，Ba含量有逐漸降低的趨勢(shì)，Zn含量為逐漸升高的趨勢(shì)，Mo隨深度含量變化不大。從而，Ba為前緣暈元素，Zn為尾暈元素，Mo為近礦暈元素，也稱主成礦元素。但是這3種元素實(shí)際上并沒(méi)有表現(xiàn)為線性增加或減少(圖3)，而是出現(xiàn)了不同程度的中轉(zhuǎn)，說(shuō)明礦床表現(xiàn)出多期次的熱液疊加現(xiàn)象。尤其是主成礦元素表現(xiàn)為強(qiáng)富集時(shí)，元素的前緣暈為強(qiáng)異常、尾暈為弱異?，F(xiàn)象，指示深部還存在著盲礦體，為深部盲礦體預(yù)測(cè)提供了依據(jù)。

表1 元素的面含量Table 1 The surface content of metal elements

3.2 模型檢驗(yàn)

3.2.1 數(shù)據(jù)后尾性檢驗(yàn)

將Mo、Ba、Zn元素含量進(jìn)行統(tǒng)計(jì)分析(表2)，由統(tǒng)計(jì)結(jié)果知Mo、Ba、Zn元素的偏度分別為5.74、0.31、0.49，全部大于0，說(shuō)明元素不滿足正態(tài)分布，均為右偏；由變異系數(shù)知，Zn值最小、最穩(wěn)定，Mo值最大、最不穩(wěn)定；由峰度知，元素的峰度均高于正態(tài)分布值3，屬于尖峰分布。由此，3種元素含量分布均為非正態(tài)分布，并且呈現(xiàn)尖峰、右偏特征。采用 QQ 圖對(duì)數(shù)據(jù)進(jìn)行后尾性分析(圖4)，發(fā)現(xiàn)3種元素仍是右偏。

表2 元素含量統(tǒng)計(jì)結(jié)果Table 2 The statistical results of element contents

3.2.2 數(shù)據(jù)平穩(wěn)性檢驗(yàn)

平穩(wěn)性檢驗(yàn)可以檢驗(yàn)樣本數(shù)據(jù)的自相關(guān)性。利用 ADF 法得出元素的平穩(wěn)性結(jié)果如表3所示，顯然，3種元素的t統(tǒng)計(jì)值均比1%的臨界值(-3.43)低，無(wú)單位根，從而它們都是平穩(wěn)序列。

表3 元素含量的ADF檢驗(yàn)Table 3 ADF test of element contents

通過(guò)上面分析，Mo、Ba、Zn元素均是后尾、非正態(tài)的平穩(wěn)序列，完全滿足所構(gòu)建的GPD分布，從而利用所設(shè)計(jì)的模型進(jìn)行異常提取。

3.3 模型求解

利用(11)式的矩法估計(jì)求出 Mo、Ba、Zn 元素的GPD分布的尺度參數(shù)分別為6.62、49.11、3.50，形狀參數(shù)分別為0.37、0.17、0.31。根據(jù)(7)式繪出Mo、Ba、Zn元素的超額均值函數(shù)圖(圖5)，通過(guò)超額均值函數(shù)圖可以確定Mo、Ba、Zn元素分別在區(qū)間[35.86，38.23]、[378.41，388.73]、[135.21，140.18]以后，MEF為斜率>0的線性函數(shù)。因此，上述區(qū)間為地球化學(xué)元素含量的初始閾值區(qū)間?，F(xiàn)在利用本文設(shè)計(jì)的修正函數(shù)確定具體的異常下限值，將這3個(gè)區(qū)間分別均勻地插入N個(gè)數(shù)，插入個(gè)數(shù)越大，閾值越精確，一般選取50個(gè)點(diǎn)就可以達(dá)到精度要求[2]。本文在3個(gè)區(qū)間上分別均勻插入50個(gè)點(diǎn)，由修正尺度函數(shù)β*(μ)得出3種元素的含量分別在36.29、382.56、138.47后，函數(shù)的ξ和β*(μ)值不再變化(圖6、圖7)。由前面的敘述可知Mo的異常下限為36.29，Ba的異常下限為382.56，Zn的異常下限為138.47。

利用(1)式得出3種元素含量GPD分布為

當(dāng)元素含量的GPD分布確定后，利用診斷性檢驗(yàn)判別所求模型是否合理，通過(guò)PP圖得出元素異常含量實(shí)際數(shù)據(jù)與GPD分布擬合結(jié)果如圖8所示。由圖8知，元素含量的實(shí)際分布均擬合在直線周圍，擬合效果較好。元素含量的理論分布和實(shí)際分布吻合，建立的模型具有有效性，閾值選取準(zhǔn)確，能夠利用所求結(jié)果進(jìn)行元素含量異?？臻g分布和礦產(chǎn)預(yù)測(cè)。

3.4 原生暈地球化學(xué)元素異常分布

相應(yīng)的地質(zhì)剖面底圖由湖北省地質(zhì)局第一地質(zhì)大隊(duì)提供(圖9)。確定元素異常下限后，利用 MapGIS 得出3種元素在地質(zhì)體中的異常分布圖(圖10、圖11、圖12)。由這3張圖可知，本文構(gòu)建的GPD異常提取模型得出的元素含量異常的空間分布規(guī)律與礦體走勢(shì)吻合， GPD 模型確定的元素含量異常分布規(guī)律(圖10)和經(jīng)典的累計(jì)頻率方法得出結(jié)果吻合(圖13)，并且精度更高。說(shuō)明本文的模型合理。根據(jù)元素空間分布規(guī)律發(fā)現(xiàn)，前緣暈 Ba 元素不僅在礦體頂部富集，在礦體下部也存在富集現(xiàn)象(圖11)，該現(xiàn)象屬于前緣暈反分帶現(xiàn)象，說(shuō)明深部可能有礦體存在；尾暈 Zn 元素主要在礦體尾部富集，但是在礦體上方有部分與前緣暈存在疊加現(xiàn)象(圖12)，說(shuō)明 Ⅶ 號(hào)礦體具有多期成礦的地球化學(xué)特征，并且尾部元素分布和前緣暈共存，進(jìn)一步說(shuō)明深部盲礦體的存在；近礦暈 Mo 元素與礦體走勢(shì)吻合(圖10)，尤其在礦體深部出現(xiàn)延伸現(xiàn)象，說(shuō)明深部有盲礦體存在，特別是位于Ⅶ號(hào)礦體東南向下延伸端，鉆孔 ZK02620 地下1 300 m 處，Mo元素具有向下和向右的強(qiáng)異常趨勢(shì)，指示盲礦體存在，并經(jīng)實(shí)際工程驗(yàn)證在該處發(fā)現(xiàn)了鉬礦和銅金礦(2016 年鉆孔 ZK02620 終孔)[1]，進(jìn)一步證明 GPD 確定的結(jié)果具有有效性。根據(jù)分布規(guī)律發(fā)現(xiàn)近礦暈元素在鉆孔 KZK10 和鉆孔 KZK11 深度 1 100 m附近有強(qiáng)異常富集。由于該處未采集到地化數(shù)據(jù)，造成該現(xiàn)象的原因很可能是位于鉆孔 ZK2620 下的鉬礦體延伸所致[1]，再結(jié)合地質(zhì)成礦規(guī)律確定此處可劃為一個(gè)找礦靶區(qū)。

4 結(jié) 論

本文以實(shí)際地質(zhì)為背景，結(jié)合GPD原理和地球化學(xué)分布規(guī)律，構(gòu)建了地球化學(xué)異常提取模型，有利于深部地球化學(xué)找礦。主要結(jié)論如下：

a.本文在地化元素含量異常分布特性和GPD理論基礎(chǔ)上構(gòu)建了異常GPD提取模型，并且設(shè)計(jì)了一種修正閾值方法，避免了異常值過(guò)大或過(guò)小的影響，提高了異常提取精度。

b.將本文構(gòu)建的模型應(yīng)用于實(shí)際礦區(qū)進(jìn)行異常提取，結(jié)果該模型能夠有效提取元素含量異常下限值，并且通過(guò)元素含量異常分布規(guī)律得出了深部地球化學(xué)找礦標(biāo)志：①近礦暈Mo 元素含量異常作為直接地球化學(xué)找礦標(biāo)志；②前緣暈Ba元素不僅在礦體頂部富集，在礦體下部也存在富集現(xiàn)象，屬于前緣暈反分帶現(xiàn)象，說(shuō)明深部有盲礦體；③尾暈、前緣元素同時(shí)在礦體中部和下部富集，說(shuō)明在礦體深部具有向下延伸的可能性。

c.本文通過(guò)構(gòu)建模型確定的元素含量異常分帶和成礦地質(zhì)條件，預(yù)測(cè)到研究區(qū)26線位于鉆孔KZK10和KZK11深度 1 100 m下有找礦靶位。

本文研究成果得益于湖北省地質(zhì)局第一地質(zhì)大隊(duì)提供原始數(shù)據(jù)，以及成都理工大學(xué)劉洪軍教授、徐爭(zhēng)啟教授提出建設(shè)性的建議，作者借此向他們深表感謝。