玻璃文物成分分析與鑒別的數學模型

張夢婷,游慧玲,張和濤,全 然

(河南工業大學a.理學院; b.化學化工學院; c.管理學院,鄭州 450001)

0 引言

在玻璃文物成分分析與鑒別研究中,大部分文獻資料是通過某些特定指標的含量直接鑒別玻璃文物所屬類型[1-4]。但基于數學模型鑒別玻璃文物類型的研究較少,只有極少數學者進行研究。如崔劍鋒等[5]采用多元統計分析方法進行探索。本研究提出并建立多種有效的數學模型對玻璃文物的化學成分進行分析及預測,鑒別未知玻璃文物類型,這對玻璃文物的鑒別及保護具有重要的現實意義[6]。

現有一批我國古代玻璃文物,考古工作者依據這些文物樣品的化學成分及其他檢測手段已將其分為高鉀玻璃和鉛鋇玻璃兩種類型。基于這批文物樣品的相關數據構建多種數學模型以解決以下4個問題。

問題一：分析玻璃文物的表面風化與玻璃類型、紋飾及顏色的關系,探究玻璃文物表面風化前后化學成分含量的統計規律,預測玻璃文物風化前的化學成分含量。

問題二：分析高鉀玻璃和鉛鋇玻璃的分類規律,對玻璃文物進行亞類劃分。

問題三：鑒別未知類別玻璃文物所屬類型,對分類結果的敏感性進行分析。

問題四：分析不同類別玻璃文物化學成分之間的關聯關系,對不同類別玻璃文物關聯關系的差異性進行判斷。

1 問題分析

對于問題一,將玻璃文物的相關數據進行預處理后,通過方差分析探究玻璃文物的表面風化與玻璃類型、紋飾及顏色的關系。構建化學成分模型,分類討論玻璃文物風化前后的化學成分含量統計規律。構建差分整合移動平均自回歸模型[7](Autoregressive integrated moving average model,ARIMA),對玻璃文物風化前的化學成分含量進行預測。

對于問題二,在問題一數據處理結果的基礎上基于主成分分析探究玻璃文物的分類規律構建線性分類規律模型,利用系統聚類模型對玻璃文物進行亞類劃分,給出具體的劃分方法及結果。

對于問題三,將問題二中的聚類中心作為特征指標,使用歐式距離衡量待檢測成分與特征指標的匹配度,鑒別未知類別玻璃文物所屬類型,分析結果的敏感性及穩定性。

對于問題四,基于K-Means++聚類模型選出初始聚類中心,基于卡林斯基-哈拉巴斯指數(Calinski-Harabazindex,CHI)和輪廓系數確定最優聚類個數,通過CHI、戴維森堡丁指數(Davies-Bouldin index,DBI)和輪廓系數3種評估指標橫向對比聚類效果,通過分類標簽的形式得到不同類別玻璃文物化學成分關聯關系的差異性。

2 模型建立與求解

2.1 玻璃文物的表面風化與玻璃類型、紋飾及顏色的關系

2.1.1 數據分析與處理

針對附件表單1,刪掉顏色數據缺失部分所在的行。針對附件表單2,將化學成分占比數據缺失的位置填充數字0,表示未檢測到該成分。由于檢測結果中各化學成分比例之和介于85%～105%之外的數據視為無效數據,故刪掉15號和17號數據。

為討論方便,對玻璃的紋飾x1、類型x2及顏色x3進行賦值量化,具體如表1所示。

表1 文本數據賦值量化結果

2.1.2 基于單因素方差分析的關聯性分析

方差分析[8]是通過對實驗數據進行分析,檢驗方差相等的多個正態總體均值是否相等,進而判斷各因素對實驗指標的影響是否顯著的一種方法。

分別建立玻璃紋飾x1、類型x2及顏色x3對玻璃文物表面風化的單因素方差分析表,計算玻璃紋飾x1、類型x2及顏色x3這3個因素的組內平均值及總均值,計算離差平方和[9],用其表示檢測結果間的差異性,再進行F檢驗。結果如表2所示。其中,SS為離差平方和,df為自由度,MS為均方,F為組間均方與組內均方之比,P-value為影響因素對玻璃文物表面風化無顯著影響的概率,F-crit為顯著性水平為0.05時的F臨界值。

表2 單因素方差分析結果

當F>F-crit時,表示該因素對玻璃文物表面風化有顯著影響。由表2可見,對于玻璃紋飾x1,由于FF-crit,故玻璃類型對玻璃文物的表面風化有顯著影響。對于玻璃顏色x3,由于F

2.1.3 風化前后玻璃文物化學成分的總體統計規律及可視化分析

將不同類型的玻璃文物按照風化前和風化后分成高鉀風化、高鉀未風化、鉛鋇風化、鉛鋇未風化4類,篩選數據分別統計討論得出：高鉀玻璃文物風化后,化學成分K2O、CaO、MgO、Al2O3、Fe2O3、CuO明顯減少;鉛鋇玻璃文物風化后,化學成分PbO、BaO明顯減少。

為使結果更加精確,通過加權平均構造玻璃文物風化前后各個化學成分含量變化的數學模型,并進行詳細討論。記aij為風化前不同類型玻璃文物第i次檢驗到第j種化學成分的含量,bij為風化后不同類型玻璃文物第i次檢驗到第j種化學成分的含量,m表示化學成分的個數,Cj為第j種化學成分風化前后變化的加權平均百分數結果。建立如下風化前后各個化學成分變化的數學模型：

(i=1, 2,…,m;j=1, 2, …,m)

(1)

由式(1)可計算出高鉀玻璃和鉛鋇玻璃在風化前后各個化學成分的變化百分比,具體如圖1和圖2所示。可見,在高鉀玻璃文物中,化學成分含量變化率較大的是SO2、CuO、BaO(按變化率從高到低的次序排列,下同);在鉛鋇玻璃文物中,化學成分含量變化率較大的是P2O5、MgO、CaO、Al2O3、Na2O,這與前述通過篩選統計得到的結果一致。

圖1 高鉀玻璃文物風化前后各化學成分含量的變化百分比

圖2 鉛鋇玻璃文物風化前后化學成分含量的變化百分比

2.1.4 基于ARIMA模型預測風化前玻璃文物的化學成分含量

為預測風化前玻璃文物的化學成分含量,結合上述分析及圖1和圖2的數據變化趨勢,識別化學成分變化為非平穩性,構建ARIMA模型進行分析。

通過做一階差分得到平穩時間序列,通過計算并分析自相關函數ACF和偏自相關函數PACF[10],分別選取最佳的ARIMA模型階層和階數,構建模型為：

Xt=c+α1Xt-1+α2Xt-2+…+αpXt-p+εt

+β1εt-1+…+βqεt-q

(2)

通過Matlab運算[11],預測出風化后的玻璃文物在風化前的化學成分,部分化學成分預測值如表3所示。

表3 部分化學成分預測值

2.2 玻璃文物的亞類劃分

2.2.1 基于PCA的玻璃文物分類特征提取

由于玻璃文物的化學成分種類多,各化學成分對分類結果的影響不盡相同,為了排除非主要成分對分類結果的影響,采用主成分分析(Principalcomponent analysis,PCA)對數據進行降維處理[12]。

針對附件表單2中14種化學成分的檢測數據,令主成分信息保留率T=0.70,特征值大于1為提取標準,共得到5個主成分指標。取5個主成分指標的貢獻率[13]為權重,構建玻璃文物的主成分分類規律模型：

F=0.30028z1+0.17211z2+0.11837z3+0.08031z4

+0.07607z5

(3)

其中,14種化學成分的主成分系數及5個主成分的貢獻率如表4所示。化學成分的主成分系數越大表明該化學成分越能代表主成分的大小;主成分的貢獻率越大表明該主成分保留的數據信息越多。由表4可以看出,主成分z1的主要化學成分有BaO、Na2O、CaO等,貢獻率為30.0278%,其保留的數據信息最多;主成分z2的主要化學成分有SrO、MgO、Al2O3等,貢獻率為17.211%;主成分z3的主要化學成分有Al2O3等,貢獻率為11.837%;主成分z4的主要化學成分有K2O、Fe2O3等,貢獻率為8.0308%;主成分z5的主要化學成分有SiO2等,貢獻率為7.607%。

表4 14種化學成分的主成分系數及5個主成分的貢獻率

2.2.2 基于聚類算法的玻璃文物亞類劃分

由高鉀玻璃文物和鉛鋇玻璃文物相關數據的散點圖可以看出,該問題適合采用系統聚類算法[10]。沿用2.1節中的數據及均值公式,分別求出風化前后高鉀玻璃文物和鉛鋇玻璃文物化學成分的均值;通過Matlab循環計算輪廓系數及聚類數[7],繪制風化高鉀玻璃文物、未風化高鉀玻璃文物、風化鉛鋇玻璃文物及未風化玻璃文物的樹狀圖,如圖3所示。

圖3 樹狀聚類圖

根據圖3得到聚類分析結果,對玻璃文物進行亞類劃分,結果如表5所示。

表5 亞類劃分結果匯總

2.3 基于聚類分析的化學成分分析與鑒別

針對附件表單3中未知類別玻璃文物的化學成分,將表5中亞類聚類中心的結果作為特征指標,使用歐氏距離[9]衡量待檢測成分與特征指標的匹配度,從而實現對未知玻璃文物的分類鑒定,結果如表6所示。

表6 未知文物的類別鑒定結果

對分類結果的敏感性進行分析。由于數據量適中,采取描述性分析方法進行檢驗,結果如表7所示。可見,數據波動幅度較大,所建模型敏感性較強。

表7 敏感性分析

為測試模型結果的穩定性,對原始未分類玻璃文物的各項化學成分進行1%、2%、5%、10%、20%、25%的擾動,結果準確,故穩定性良好。

2.4 基于K-Means++聚類算法的玻璃文物化學成分關聯度分析

2.4.1 基于K-Means++的初始聚類中心選取

K-Means算法中的k個初始聚類中心的選擇對最終聚類結果及運行時間有極大的影響。傳統K-Means算法的聚類效果依賴于聚類中心的初始化,而K-Means++聚類算法[14]以“使初始聚類中心間的相互距離盡可能的遠”為基本原則,對初始聚類中心的選擇方法進行優化,從而提高聚類精度,加快收斂速度。

給出K-Means++聚類算法選取初始聚類中心的算法步驟：

步驟1：在數據集X中隨機令一個樣本點作為初始聚類中心C1。

步驟2：計算每個樣品點與初始聚類中心C1的距離D(xm)。

步驟3：從數據集X中隨機選擇下一個聚類中心Cj,每個樣品點被選擇的概率為：

(4)

其中,M(xj)為Cj與距樣本點最近的選擇點之間的距離。令P(x)最大時對應的樣本點為新的初始聚類中心。

步驟4：重復步驟2和步驟3,直至選出k個初始聚類中心。

2.4.2 最優聚類數的確定

采用CHI和輪廓系數確定最優聚類數。CHI的本質為類間距離與類內距離的比值。將容量為N的數據集X聚成K類,用各個類中心點與數據集X中心點的距離平方和來度量類間距離Bk,用類內各點與類中心點的距離平方和來度量類內距離Wk。CHI的計算公式為：

(5)

其中,cq為類q的中心點,ce為數據集X的中心點,nq為某一類q中的數據量,Qq為類q的數據集。

輪廓系數可用于評價聚類效果,其范圍在[-1, 1]。同類樣本距離越近,不同類樣本距離越遠輪廓系數越大,說明聚類效果越合理。令某個樣本的輪廓系數為S,聚類總的輪廓系數為SC,具體如下：

(6)

其中,a表示一個樣本與其所在簇內其他樣本的平均距離,b代表一個樣本與其他簇樣本的平均距離。

運用Matlab計算并制成簇計算標準值,如圖4所示。圖中分類所代表的柱最高時聚類的個數最優。得出：高鉀未風化玻璃文物在k=3時最優,高鉀風化玻璃文物在k=4時最優,鉛鋇未風化玻璃文物在k=5時最優,鉛鋇風化玻璃文物在k=4時最優。

圖4 簇計算標準值

2.4.3 基于K-Means++聚類算法的玻璃文物化學成分關聯度分析模型

為對數據進行可視化分析,構建了玻璃文物化學成分的關聯度分析模型。運用Matlab制圖,結果如圖5所示。圖5中每一行、每一列的聚類均較為分散,故聚類效果較好。

圖5 最優聚類可視化分析

通過表8對比可得：高鉀玻璃文物風化前后的分類只有K2O、CaO、CuO發生改變,而鉛鋇玻璃文物全部化學成分都發生了改變。

表8 分類標簽

采用CHI、DBI、輪廓系數[15]3個評價指標分析聚類效果,如表9所示。發現4個分類下的輪廓系數均接近1,說明簇內距離小,簇間距離大,CHI同理。各項DBI都較小,說明其類內距離較小,類間距離較大,即不同類別之間的相似度較小,說明聚類結果的質量較好,各類別之間差異較明顯。

表9 聚類效果分析

2.4.4 不同類別之間化學成分關聯關系的差異性判斷

為了得到更加準確的聚類結果,利用Matlab運行算法4次后,將不同類別玻璃文物風化前后化學成分的4種聚類結果以分類標簽的形式匯總,結果如表10所示。

表10 聚類結果的分類標簽形式匯總

通過表10可以看出,只有極少次數的聚類會使聚類結果發生明顯變化。其中,鉛鋇未風化玻璃文物聚類數為5與高鉀風化玻璃文物聚類數為4時對應的聚類結果相較其他結果差異性較大,而其他大多數聚類結果都圍繞著眾數上下波動或保持不變,但一般會高于正常值,說明不同距離造成的影響具有一定的同質性,但弱于正常值。

3 結論

針對玻璃文物表面風化問題,通過單因素方差分析研究發現：玻璃類型顯著影響玻璃文物的表面風化程度,而玻璃紋飾和顏色影響不顯著。建化學成分模型分析得到：風化前后,高鉀玻璃文物SO2、CuO、BaO等化學成分含量變化率較大,鉛鋇玻璃文物P2O5、MgO、CaO、Al2O3、Na2O等化學成分含量變化率較大。通過ARIMA模型預測得到風化前多個玻璃文物的化學成分含量。

針對玻璃文物的分類問題,通過主成分分析確定化學成分的分類規律,通過系統聚類算法對4種類型的玻璃文物進行了亞類劃分,具有良好的合理性及穩定性。

針對未知類別玻璃文物的化學成分分析與類別鑒定問題,基于亞類劃分結果構建化學成分識別模型,確定了8種未知類別文物的類別,對匹配度進行描述統計并對化學成分數據進行5%～30%的擾動,結果表明,該模型具有較好的穩定性及敏感性。

針對不同類別玻璃文物化學成分的關聯性及差異性問題,通過K-Means++聚類算法描述玻璃文物化學成分間的關聯性,通過多次聚類揭示4類玻璃文物間的差異。

通過以上分析可知,所建數學模型合理,能實現對玻璃文物成分的有效分析,可為其他玻璃文物成分研究提供參考。