基于K-Means聚類與灰色關(guān)聯(lián)分析對玻璃文物成分的研究

徐 惠，胡 珊，姚旭敏，儲昭順

（安徽財經(jīng)大學(xué)統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030）

針對古代玻璃制品的成分分析與鑒別，李青會等人于2006年對戰(zhàn)國時期的玻璃珠進行研究，利用現(xiàn)代技術(shù)，得出戰(zhàn)國時期中國境內(nèi)同時存在三種硅酸鹽玻璃，中國古代的PbO、BaO、SiO2和K2O、SiO2玻璃在技術(shù)發(fā)展上應(yīng)該具有密切聯(lián)系的結(jié)論［1］；又于2007年利用相關(guān)技術(shù)對中國古代玻璃的化學(xué)成分進行分析，指出應(yīng)加強對西周到戰(zhàn)國時期中國出土的釉砂、玻砂、鑲嵌玻璃珠，以及伴隨古代玻璃同時出土的綠松石等文物的科技研究［2］。

對于近代玻璃的研究，趙娟等人于2002年在相關(guān)組分制造硒硫化鎘顏色玻璃中的作用研究中，指出在制造硒硫化鎘顏色玻璃時，除了加入玻璃的基本組成SiO2和著色劑CdS、Se 外，還必須加入ZnO、堿金屬氧化物、冰晶石、B2O3和重金屬硫化物等次要組分［3］；王承遇等人于2003年綜述了影響浮法玻璃等的耐風(fēng)化性因素，風(fēng)化產(chǎn)物的形貌和風(fēng)化過程，得出了玻璃表面風(fēng)化析堿量隨風(fēng)化溫度、濕度和時間而增加的結(jié)論［4］。

文章針對古代玻璃制品的研究，借鑒近代玻璃的研究成果，立足于新的數(shù)據(jù)，從成分分析與鑒別方面進行研究，擬運用K-均值聚類、灰色關(guān)聯(lián)分析、BP神經(jīng)網(wǎng)絡(luò)等方法建立模型，研究高鉀類型與鉛鋇類型玻璃的分類依據(jù)，以及兩類玻璃文物風(fēng)化的難易程度。

1 數(shù)據(jù)來源和假設(shè)

文章數(shù)據(jù)來源于2022年全國大學(xué)生數(shù)學(xué)建模競賽。為便于解決問題，研究過程和結(jié)論均建立在以下假設(shè)成立的基礎(chǔ)上：（1）顏色空白數(shù)據(jù)的劃分整合對關(guān)系判斷無影響；（2）剔除風(fēng)化玻璃選出未風(fēng)化（未風(fēng)化玻璃選出風(fēng)化）數(shù)據(jù)對結(jié)論無影響；（3）含量預(yù)測時單一變量預(yù)測不會影響最終預(yù)測結(jié)果；（4）預(yù)測值不真實，為考慮真實性，使用原始數(shù)據(jù)進行BP神經(jīng)網(wǎng)絡(luò)模型對化學(xué)成分分析無影響。

2 基于隨機森林插補的K-Means聚類分析模型

2.1 隨機森林插補法數(shù)據(jù)處理

根據(jù)玻璃制品的相關(guān)數(shù)據(jù)，分析表面風(fēng)化與其他三個變量的關(guān)系，比較其出現(xiàn)頻率，進行四個變量之間的卡方檢驗，分析變量之間的差異性。研究玻璃文物表面有無風(fēng)化以及其化學(xué)成分含量的統(tǒng)計規(guī)律。首先，根據(jù)表單數(shù)據(jù)將玻璃類型分為高鉀和鉛鋇類，控制變量分析其他數(shù)據(jù)，改變前幾個變量，分析元素的平均值，多個控制變量比較分析。其次，使用缺失森林插補法進行模擬預(yù)測，預(yù)測風(fēng)化前化學(xué)成分含量［5］，發(fā)現(xiàn)未風(fēng)化的高鉀與鉛鋇玻璃數(shù)量相同，而風(fēng)化的高鉀玻璃明顯少于鉛鋇玻璃；所以高鉀類型玻璃相對較難風(fēng)化，鉛鋇類型玻璃相對易于風(fēng)化。對所有數(shù)據(jù)按照風(fēng)化、無風(fēng)化進行劃分，統(tǒng)計不同類型、紋飾、顏色的成分總含量。對類型、紋飾與表面風(fēng)化進行分類，計算均值，發(fā)現(xiàn)SnO2只存在于鉛鋇風(fēng)化紋飾C 中，無風(fēng)化玻璃中含量較多的都是活潑型較低的成分（K2O 除外），含量相似玻璃中鉛鋇紋飾A 和C 所含的氧化物種類與數(shù)量大致相同；高鉀紋飾B 型不含有的氧化物（等于0）較多，說明在風(fēng)化過程中高鉀玻璃與環(huán)境元素交換時，交換的元素較少，較穩(wěn)定。高鉀玻璃在風(fēng)化后的主要化學(xué)成分含量減少；鉛鋇玻璃在風(fēng)化后主要化學(xué)成分含量增加（SiO2的含量變化相反）。最后，將總體數(shù)據(jù)分為高鉀和鉛鋇類型，視高鉀和鉛鋇風(fēng)化時期的數(shù)據(jù)為缺失，根據(jù)隨機森林算法對缺失值進行填充，得到補充完整的數(shù)據(jù)。利用補充的高鉀和鉛鋇風(fēng)化數(shù)據(jù)預(yù)測未風(fēng)化的數(shù)據(jù)。

2.2 K-Means聚類研究

分析高鉀和鉛鋇玻璃的分類規(guī)律，使用K-Means分析方法顯現(xiàn)分類標準。對高鉀和鉛鋇玻璃亞類進行劃分，建立了分層聚類模型。檢查兩類數(shù)據(jù)的缺失值情況，初步判斷聚類的可信率；根據(jù)聚類樹狀圖進行分析，對聚類情況和數(shù)據(jù)進行探測，進行亞類劃分。

依據(jù)單因素方差進行合理性分析，分析高鉀與鉛鋇類的氧化物，將標準差較大數(shù)值的氧化物與亞分類的氧化物進行比較，觀察是否都存在；改變其中一項氧化物的數(shù)值重新分層聚類，劃分亞類與之前的比較，若改變數(shù)值的氧化物不出現(xiàn)在亞類中，說明檢驗的敏感性。

分析表單數(shù)據(jù)，確定類別1為高鉀，類別2為鉛鋇；對數(shù)據(jù)進行K-Means聚類分析，得到字段差異分析和聚類中心點坐標，以聚類1 的中心點坐標作為橫坐標，建立成分中心散點圖，如圖1所示?？傮w與y=x進行對比，形成成分比較分析散點圖，取出與y=x的斜率相差較大的點進行分析。

圖1 成分中心散點圖

利用單因素方差分析得出分析圖，如圖2所示。分析含量較大的成分，二氧化硅（SiO2）、氧化鉀（K2O）的含量在高鉀類含量較多，氧化鉛（PbO）、氧化鋇（BaO）在鉛鋇中含量較多。這四種氧化物的標準差較大，數(shù)據(jù)較為分散。

圖2 單因素方差分析圖

2.3 分層聚類分析

采用歐氏距離作為指標，對關(guān)聯(lián)度較大的前k種氣象因素進行分層K-Means 聚類（k的值取決于樣本的多少以及精度的需求），聚類完成后得到多個化學(xué)成分樣本子集。在確定了該相似性距離為聚類指標后，采用結(jié)合歐氏距離和ρ(X,Y)為標準的多重聚類方法［6］。

在X和Y歸一化之后，其之間的協(xié)方差ρ(X,Y)表示為

建立分層聚類模型為

其中，yi代表在一系列觀測變量中的個體值，K代表群的個數(shù)，πk表示一個個體值屬于k群（或k群的大?。┑南闰灨怕?，θ表示模型的參數(shù)，f(yi|θ)是指當(dāng)特定集群的混合密度θ作為模型參數(shù)時，yi的分布情況。

同樣地，弗萊利和拉夫特里提出可以將上式中的模型表現(xiàn)為下述的相似形式：

其中，f(yi|θ)符合多元正態(tài)（高斯）分布φk，參數(shù)包括平均值μk和方差矩陣∑k.

總體數(shù)據(jù)分為高鉀和鉛鋇類型，分別運用Matlab 對其進行分層聚類；檢查數(shù)據(jù)的缺失值情況，初步判斷聚類的可信率；根據(jù)聚類表了解樣本之間的距離和聚類情況；高鉀和鉛鋇類型的聚類樹狀圖如圖3 所示，分析聚類樹狀圖，對聚類情況和數(shù)據(jù)進行探測，發(fā)現(xiàn)各類之間的層次關(guān)系，對數(shù)據(jù)進行亞類劃分。

圖3 高鉀（左）和鉛鋇（右）類型的聚類樹狀圖

根據(jù)聚類系數(shù)與樹狀圖進行亞類劃分，高鉀類型：高鉀主要分為二氧化硅（SiO2）與類別1，類別1分為氧化銅（CuO）與類別2，類別2 分為氧化鉀（K2O）與類別3，類別3 分為氧化鈣（CaO）與氧化鋁（Al2O3）。鉛鋇類型：鉛鋇主要分為二氧化硅（SiO2）、氧化鉛（PbO）與類別1，類別1分為氧化鋇（BaO）與類別2，類別2分為氧化鋁（Al2O3）與類別3，類別3分為氧化銅（CuO）與氧化鈣（CaO）。

2.4 基于單因素方差分析的模型檢驗

將分類數(shù)據(jù)進行單因素方差分析，得到方差齊性檢驗，如表1所示。觀測變量標準差，如果標準差較大，說明觀測變量的變動主要是由控制變量引起的，可以主要對觀測變量進行解釋，控制變量給觀測變量帶來了顯著影響，即控制變量對于觀測變量具有敏感性。

表1 方差齊性檢驗

由表1 可知，高鉀類標準差較大的有二氧化硅（SiO2）、氧化鉀（K2O）、氧化鈣（CaO）、氧化鋁（Al2O3）、氧化鐵（Fe2O3），相對于高鉀分類主要結(jié)果多出氧化鐵（Fe2O3）；鉛鋇類標準差較大的有二氧化硅（SiO2）、氧化鈣（CaO）、氧化鋁（Al2O3）、氧化銅（CuO）、氧化鉛（PbO）、氧化鋇（BaO）、五氧化二磷（P2O5），相對于鉛鋇分類主要結(jié)果多出五氧化二磷（P2O5），但所有分類結(jié)果都完整，展示出分類結(jié)果的合理性。

3 基于主成分分析的玻璃文物成分分析

3.1 研究思路

為便于分析，將數(shù)據(jù)表單分為風(fēng)化與未風(fēng)化兩個部分，通過主成分分析提取主成分，得出所有主成分貢獻率均小于30%，除第一主成分大于20%外，其余主成分均小于20%，剔除其他成分形成數(shù)據(jù)表，利用BP 神經(jīng)網(wǎng)絡(luò)模型進行二次鑒別；在BP 神經(jīng)網(wǎng)絡(luò)模型中，去除變量氧化鈉的影響，研究其他變量的特征，得到鑒別玻璃類別結(jié)果。

3.2 研究方法

3.2.1 主成分分析模型

對表單三數(shù)據(jù)進行主成分分析，提取主成分。使用主成分分析法計算各個主成分的貢獻率，得到特征值、貢獻率和累計貢獻率，得出所有主成分貢獻率均小于30%，除第一主成分大于20%外，其余主成分均小于20%，如表2所示。

重點研究第一主成分各個變量的特征向量。通過對第一主成分表達式的分析，可以得出氧化鈉變量的系數(shù)，相較于其他變量，小于1個數(shù)量級；同時氧化鈉的數(shù)值大多數(shù)為空，故剔除變量氧化鈉的數(shù)據(jù)。

3.2.2 BP神經(jīng)網(wǎng)絡(luò)模型

在BP 神經(jīng)網(wǎng)絡(luò)中，相鄰層的神經(jīng)元進行全連接，每層各個神經(jīng)元之間無連接，網(wǎng)絡(luò)按照有監(jiān)督方式學(xué)習(xí)，當(dāng)一對學(xué)習(xí)模式提供網(wǎng)絡(luò)后，各神經(jīng)元獲得網(wǎng)絡(luò)的輸入響應(yīng)產(chǎn)生連接權(quán)值。然后按減少希望輸出與實際輸出誤差的方向，從輸出層經(jīng)各中間層逐層修正各連接權(quán)值，回到輸入層［7］。

在BP神經(jīng)網(wǎng)絡(luò)模型中，基于Python去除變量氧化鈉的影響［8］，研究其他變量的特征，從而達到鑒別玻璃類別的目的。為便于推斷玻璃類別，將數(shù)據(jù)分為風(fēng)化、未風(fēng)化兩個部分數(shù)據(jù)，分別進行BP 神經(jīng)網(wǎng)絡(luò)模型檢驗；得到鑒別結(jié)果，其中A1、A6、A7為高鉀玻璃，A2、A3、A4、A5、A8為鉛鋇玻璃。

3.3 模型檢驗

對于表中數(shù)據(jù)運用SPSS 進行量化分析，分別判斷樣本所屬大類（類別1：高鉀；類別2：鉛鋇），得到效應(yīng)量化分析表，如表3所示。進而對不同氧化物類別進行單因素方差分析，用于檢驗高鉀和鉛鋇類別數(shù)據(jù)是否存在顯著性差異。

表3 效應(yīng)量化分析表

效應(yīng)量化分析的結(jié)果顯示，基于二氧化硅（SiO2）、氧化鈉（Na2O）、氧化鈣（CaO）、氧化鎂（MgO）、氧化鋁（Al2O3）、氧化鐵（Fe2O3）、氧化鉛（PbO）、氧化鋇（BaO）、五氧化二磷（P2O5）、氧化鍶（SrO）、氧化錫（SnO2），對于Eta方（η2值）較大，數(shù)據(jù)的差異來源于不同組別間的差異；Cohen′s f值較大，即數(shù)據(jù)的效應(yīng)量化的差異程度為大程度差異，分析效應(yīng)量化數(shù)據(jù)Eta方（η2值），與Cohen′s f值，對于分類結(jié)果氧化物的敏感程度較高。

4 灰色關(guān)聯(lián)分析模型

4.1 研究思路

通過建立灰色關(guān)聯(lián)度模型，代入亞類劃分出的高鉀和鉛鋇3個類別中的氧化物數(shù)值，進行灰色關(guān)聯(lián)性分析，得到氧化物與高鉀、鉛鋇之間的關(guān)聯(lián)度；然后對兩種類型的化學(xué)成分進行多配對樣本Friedman 檢驗，具有顯著性的使用Friedman 檢驗，不具有顯著性的使用正態(tài)性檢驗直方圖，得到氧化物之間的關(guān)聯(lián)與差異幅度；最后進行灰色關(guān)聯(lián)分析，得到氧化物之間的灰色系數(shù)表，分析高鉀與鉛鋇類型中氧化物之間的線性與非線性關(guān)系。對得到的關(guān)聯(lián)關(guān)系進行高鉀類型與鉛鋇類型的比較，得出差異性。

4.2 研究方法

4.2.1 灰色關(guān)聯(lián)分析模型

針對數(shù)據(jù)進行無量綱化處理（均值化、初值化），求解母序列（對比序列）和特征序列之間的灰色關(guān)聯(lián)系數(shù)值和灰色關(guān)聯(lián)度值，對灰色關(guān)聯(lián)度值進行排序，得出結(jié)論，即設(shè)系統(tǒng)特征行為序列為

系統(tǒng)的相關(guān)因素行為序列為

記折線

為(Xi)°

令

則灰色絕對關(guān)聯(lián)度為

計算關(guān)聯(lián)系數(shù)為

則關(guān)聯(lián)度為

灰色系統(tǒng)理論著重考慮點點之間的距離遠近對關(guān)聯(lián)度的影響。其中，ρ稱為分辨系數(shù)，一般情況下，取ρ= 0.5.對于ρ取值的一般原則，避免了系統(tǒng)因子觀測序列的異常值支配整個系統(tǒng)關(guān)聯(lián)度取值的情況，能夠使關(guān)聯(lián)度更好地體現(xiàn)系統(tǒng)的整體性［9］。根據(jù)觀測值動態(tài)變化選取分辨系數(shù)ρ的值，使其取值具有一定的客觀基礎(chǔ)，具體取值規(guī)律如下：

記?v為所有差值絕對值的均值，即

記ε?=，則ρ的取值為

當(dāng)?max > 3?v時，ε?≤ρ≤1.5ε?；

當(dāng)?max ≤3?v時，1.5ε?≤ρ≤2ε?.

注：分辨系數(shù)ρ∈（0，∞），ρ越小，分辨力越大，一般ρ的取值區(qū)間為（0，1），具體取值可視情況而定。當(dāng)ρ≤0.5463時，分辨力最好，通常取ρ=0.5.

關(guān)聯(lián)度表示各評價項與“參考值”（母序列）之間的相似關(guān)聯(lián)程度，其是由關(guān)聯(lián)系數(shù)進行計算平均值得出，關(guān)聯(lián)度值介于0～1之間，該值越大表示評價項與“參考值”（母序列）相關(guān)性越強，關(guān)聯(lián)度越高，意味著評價項與“參考值”（母序列）之間關(guān)系越緊密，因而其評價越高。結(jié)合關(guān)聯(lián)度值，針對所有評價項進行排序，得到各評價項排名［10］。

結(jié)合上述關(guān)聯(lián)系數(shù)結(jié)果進行加權(quán)處理，最終得出關(guān)聯(lián)度值，使用關(guān)聯(lián)度值針對14 個評價對象進行評價排序，如表4所示。

表4 灰色關(guān)聯(lián)度分析

從高鉀關(guān)聯(lián)度表可以看出：針對高鉀類別內(nèi)部氧化物本次5個評價項，分析可得評價氧化鉛（PbO）最高（關(guān)聯(lián)度為：0.807），其次是二氧化硅（SiO2）（關(guān)聯(lián)度為：0.614）；從鉛鋇關(guān)聯(lián)度表可以看出：針對鉛鋇類別內(nèi)部氧化物本次6 個評價項，分析可得二氧化硅（SiO2）評價最高（關(guān)聯(lián)度為：0.868），其次是氧化鋁（Al2O3）（關(guān)聯(lián)度為：0.860）。

4.2.2 多配對樣本Friedman檢驗?zāi)Ｐ?/p>

高鉀類型二氧化硅（SiO2）、氧化鉀（K2O）、氧化鈣（CaO）顯著性P 值<0.05，水平上呈現(xiàn)顯著性，拒絕原假設(shè)，因此數(shù)據(jù)不滿足正態(tài)分布，可以進行Friedman 檢驗［11］，如表5 所示。氧化鋁（Al2O3）、氧化銅（CuO）顯著性P值>0.05，水平上不呈現(xiàn)顯著性，不能拒絕原假設(shè)，因此數(shù)據(jù)滿足正態(tài)分布，建議采用方差分析［12］。

表5 高鉀Friedman檢驗分析結(jié)果表

通過分析可知，顯著性P 值為0.000***，因此統(tǒng)計結(jié)果顯著，說明二氧化硅（SiO2）、氧化鉀（K2O）、氧化鈣（CaO）之間存在顯著差異；其差異幅度Cohen′s f值為：4.13，差異幅度非常大。

鉛鋇類型二氧化硅（SiO2）、氧化鈣（CaO）、氧化鋁（Al2O3）、氧化銅（CuO）、氧化鋇（BaO）顯著性P值<0.05，水平上呈現(xiàn)顯著性，拒絕原假設(shè)，因此數(shù)據(jù)不滿足正態(tài)分布，可以進行Friedman 檢驗，如表6 所示。氧化鉛（PbO）顯著性P值>0.05，水平上不呈現(xiàn)顯著性，不能拒絕原假設(shè)，因此數(shù)據(jù)滿足正態(tài)分布，建議采用方差分析。

表6 鉛鋇Friedman檢驗分析結(jié)果表

顯著性P值為0.000***，因此統(tǒng)計結(jié)果顯著，說明二氧化硅（SiO2）、氧化鈣（CaO）、氧化鋁（Al2O3）、氧化銅（CuO）、氧化鋇（BaO）之間存在顯著差異；其差異幅度Cohen′s f值為：1.55，差異幅度較大。

鉛鋇類型氧化鉛（PbO）正態(tài)圖基本上呈現(xiàn)中型，數(shù)據(jù)基本可接受為正態(tài)，所以與其他化學(xué)成分之間差異性較小。

4.2.3 結(jié)果分析

根據(jù)高鉀玻璃和鉛鋇玻璃亞類劃分后的氧化物數(shù)據(jù)，進行灰色關(guān)聯(lián)度相關(guān)性分析，得到高鉀和鉛鋇相關(guān)系數(shù)，如表7所示。

表7 高鉀和鉛鋇相關(guān)系數(shù)

根據(jù)高鉀關(guān)聯(lián)度相關(guān)性分析，得到氧化物二氧化硅（SiO2）與氧化鉀（K2O）、氧化鈣（CaO）、氧化鋁（Al2O3）、氧化銅（CuO）呈負指數(shù)相關(guān)；氧化物氧化鉀（K2O）與氧化鈣（CaO）、氧化鋁（Al2O3）、氧化銅（CuO）呈正指數(shù)相關(guān)；氧化物氧化鈣（CaO）與氧化鋁（Al2O3）、氧化銅（CuO）呈正指數(shù)相關(guān)；氧化物氧化鋁（Al2O3）與氧化銅（CuO）呈正指數(shù)相關(guān)。

根據(jù)鉛鋇關(guān)聯(lián)度相關(guān)性分析，得到氧化物二氧化硅（SiO2）與氧化鈣（CaO）、氧化鉛（PbO）、氧化鋇（BaO）、氧化銅（CuO）呈負指數(shù)相關(guān)，與氧化鋁（Al2O3）呈正指數(shù)相關(guān)；氧化物氧化鈣（CaO）與氧化鋇（BaO）、氧化銅（CuO）呈負指數(shù)相關(guān)，與氧化鋁（Al2O3）、氧化鉛（PbO）呈正指數(shù)相關(guān)；氧化物氧化鋁（Al2O3）與氧化鉛（PbO）、氧化鋇（BaO）、氧化銅（CuO）呈負指數(shù)相關(guān)；氧化物氧化鉛（PbO）與氧化鋇（BaO）、氧化銅（CuO）呈負指數(shù)相關(guān)；氧化物氧化鋇（BaO）與氧化銅（CuO）呈正指數(shù)相關(guān)。

5 結(jié)論

通過分析得到高鉀與鉛鋇的灰色關(guān)聯(lián)度，進行相關(guān)性分析得到各類氧化物的正負相關(guān)指數(shù)，分析得到高鉀與鉛鋇大類的氧化物之間正負相關(guān)指數(shù)交替出現(xiàn)。對于高鉀玻璃類別氧化物正相關(guān)指數(shù)偏多，一類氧化物往往會帶動另類氧化物的同向含量變化，氧化物總體呈現(xiàn)方驂并路趨勢；對于鉛鋇玻璃類別氧化物負指數(shù)相關(guān)偏多，一類氧化物往往會帶動另類氧化物的異向含量變化，氧化物總體呈現(xiàn)此消彼長趨勢。