杜世強(qiáng),樂靖雯,周海麗,楊松杰

(西北民族大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,甘肅 蘭州 730030)

0 引言

古代玻璃制品是考證絲綢之路歷史發(fā)展變化的重要物品,而研究它的埋藏環(huán)境與風(fēng)化之間的關(guān)系,對(duì)于充分挖掘中國古代對(duì)外商貿(mào)的歷史淵源具有重要意義.文獻(xiàn)[1]采用多元統(tǒng)計(jì)方法對(duì)四川、重慶、貴州、廣西、廣東地區(qū)出土的100余件古玻璃樣品化學(xué)成分?jǐn)?shù)據(jù)進(jìn)行分析,以期能夠提示我國漢代南方和西南地區(qū)的玻璃生產(chǎn)情況.文獻(xiàn)[2]對(duì)廣州出土的西漢早期至東漢的46件玻璃單色珠、耳珰等器物的化學(xué)成分進(jìn)行了成分分析,探討表面風(fēng)化對(duì)化學(xué)成分定量分析的影響,并依據(jù)CaO、Al2O3含量,結(jié)合微量元素Rb和Sr的比例,對(duì)所占比例最高的鉀硅酸鹽玻璃進(jìn)行了亞類劃分.文獻(xiàn)[3]對(duì)20多種較有代表性的中國古代鉛玻璃樣品進(jìn)行實(shí)物考察和分析研究,根據(jù)化學(xué)檢測(cè)結(jié)果,這批古玻璃可分為五個(gè)成分系統(tǒng),并可歸納為三類:鉛鋇玻璃類有PbO-BaO-SiO2和Na2O-PbO-BaO-SiO2,鉛硅玻璃類有PbO-SiO2,堿鉛硅玻璃類有K2O-PbO-SiO2和Na2O-PbO-SiO2.文獻(xiàn)[4]利用EDXRF探針無損分析技術(shù),分析了河南禹縣陽翟遺址出土的一批金元時(shí)期的玻璃樣品,結(jié)果表明,這批古玻璃種類較為復(fù)雜,主要包括鉀鈣玻璃、鈉鈣玻璃和鉛鉀鈣玻璃等.文獻(xiàn)[1-4]均采用化學(xué)方法對(duì)古玻璃進(jìn)行性能測(cè)定,其中文獻(xiàn)[2-3]主要通過主觀判斷來劃分古玻璃類別,文獻(xiàn)[1,4]則采用聚類分析和因子分析驗(yàn)證了分類結(jié)果.如果對(duì)大批量樣品的組成成分進(jìn)行分類時(shí),單憑主觀判斷易造成誤差,且缺乏說服力,所以本文在參考文獻(xiàn)[1-4]分類結(jié)果的基礎(chǔ)上,利用多種統(tǒng)計(jì)分析方法進(jìn)行數(shù)據(jù)分析,從而得出較為合理的對(duì)古玻璃分類的標(biāo)準(zhǔn).

利用卡方檢驗(yàn)和相關(guān)性分析玻璃文物表面風(fēng)化和類型、紋飾、顏色的關(guān)系以及所含化學(xué)成分之間的相關(guān)性.結(jié)合文獻(xiàn)[4],對(duì)化學(xué)成分的主要用途進(jìn)行分類,再根據(jù)主要化學(xué)成分建立層次聚類和K-means聚類模型,得出玻璃文物分類、亞分類標(biāo)準(zhǔn),最后構(gòu)建隨機(jī)森林分類器模型,進(jìn)而對(duì)未知類型玻璃進(jìn)行分類.

1 相關(guān)因素分析

1.1 數(shù)據(jù)預(yù)處理

分析2022年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽C題《古代玻璃制品的成分分析與鑒別》數(shù)據(jù),發(fā)現(xiàn)玻璃制品顏色列存在缺失值,檢測(cè)的化學(xué)成分存在較多空值.依據(jù)2022年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽C題所給背景條件,將化學(xué)成分比例累加和介于85%～105%之間的數(shù)據(jù)視為有效數(shù)據(jù),通過計(jì)算,編號(hào)15和編號(hào)17為無效數(shù)據(jù),把這兩個(gè)編號(hào)從數(shù)據(jù)表中剔除,對(duì)于未檢測(cè)到成分而缺失的空值,用“0”填充.對(duì)顏色列缺失,采用熱卡填充方法對(duì)含有缺失值元組的其他屬性與其他完整數(shù)據(jù)元組的屬性求歐氏距離,將距離最短的元組中對(duì)應(yīng)的值作為缺失值的估計(jì)量進(jìn)行填充,計(jì)算公式如式(1).

其中,L表示歐氏距離,xi表示待填充的數(shù)據(jù)屬性,yi表示完整的數(shù)據(jù)屬性.經(jīng)計(jì)算,“顏色”列中的空白處依次填充為“黑”“紫”“淺藍(lán)”“淺綠”.

1.2 玻璃文物表面風(fēng)化與其類型、紋飾和顏色的相關(guān)性分析

其中,A表示實(shí)際頻數(shù),T表示理論頻數(shù),且:

其中,R表示行,C表示列,TRC表示第R行和第C列的理論數(shù),nR,nC表示第R行和第C列的合計(jì)數(shù),N表示總的合計(jì)數(shù).最后可根據(jù)χ2值查閱卡方界值表,得到P值［5］.

以風(fēng)化情況作為變量X,以顏色、紋飾和類型作為變量Y,利用python編程實(shí)現(xiàn)上述模型,得卡方檢驗(yàn)表1.

表1 表面風(fēng)化與紋飾、玻璃類型、顏色相關(guān)性卡方檢驗(yàn)表

從表1可看出,玻璃風(fēng)化情況與紋飾、顏色不存在顯著差異,但與類型存在顯著差異.

1.3 玻璃文物中化學(xué)成分間的相關(guān)性分析模型建立與求解

由協(xié)方差方程Cor(X,Y)=E{[X-E(X)][Y-E(Y)]}可知其相關(guān)系數(shù),14種化學(xué)成分的相關(guān)性分析模型如式(4):

其中,i=1,2,…,14,j=1,2,…,14.

對(duì)模型進(jìn)行求解,得出不同類別玻璃文物化學(xué)成分之間的相關(guān)系數(shù),如圖1所示.

高鉀玻璃中氧化鍶與五氧化二磷間相關(guān)系數(shù)為0.78,呈顯著正相關(guān);氧化鐵和五氧化二磷間相關(guān)系數(shù)為0.81,呈顯著正相關(guān);二氧化硅和氧化鉀間相關(guān)系數(shù)為-0.86,呈顯著負(fù)相關(guān).而氧化鍶和氧化鈣,五氧化二磷和氧化鉛、氧化錫,氧化鈉和氧化銅、氧化鐵,相關(guān)性最低,獨(dú)立性最強(qiáng).

鉛鋇玻璃中氧化鉛與氧化鍶間相關(guān)系數(shù)為0.63,呈顯著正相關(guān);二氧化硫與氧化鉛間相關(guān)系數(shù)為-0.85,呈顯著負(fù)相關(guān).氧化銅和氧化鐵,氧化鎂和二氧化硫,氧化鈉和氧化鍶相關(guān)性最低,獨(dú)立性最強(qiáng).

圖1 高鉀玻璃和鉛鋇玻璃化學(xué)成分相關(guān)系數(shù)圖

2 玻璃文物中化學(xué)成分預(yù)測(cè)

文物在燒制過程中使用的原材料和助熔劑會(huì)直接與空氣發(fā)生快速氧化還原反應(yīng),使得文物表面本身的鉀元素或鉛元素含量在最開始時(shí)就達(dá)到最大值.文物經(jīng)過長時(shí)間埋藏,表面化學(xué)元素與周圍環(huán)境所含化學(xué)元素進(jìn)行大量交換,最終文物和埋藏環(huán)境的化學(xué)元素含量達(dá)到平均水平,即出土?xí)r間較長文物表面化學(xué)元素含量占比少于最開始出土?xí)r的化學(xué)元素含量,但仍具有可識(shí)別性.

鉛鋇玻璃前三類占比最高的化學(xué)成分為:氧化鉛(PbO)、氧化鋇(BaO)、氧化鍶(SrO);高鉀玻璃前兩類占比最高的化學(xué)成分為:氧化鉀(K2O)、氧化鐵(Fe2O3).高鉀玻璃預(yù)測(cè)數(shù)據(jù)部分選用這兩個(gè)強(qiáng)相關(guān)性化學(xué)元素,鉛鋇玻璃預(yù)測(cè)數(shù)據(jù)部分選用這三個(gè)強(qiáng)相關(guān)性化學(xué)元素.

2.1 風(fēng)化前化學(xué)成分預(yù)測(cè)模型的建立與求解

基于變量設(shè)定,構(gòu)建一個(gè)二項(xiàng)邏輯回歸模型,對(duì)影響文物類別劃分的因素進(jìn)行計(jì)量分析,規(guī)定因變量為0-1型二分變量,其中“1”表示為無風(fēng)化,“0”表示為風(fēng)化,解釋變量為所給的化學(xué)成分含量數(shù)據(jù)類別.

Logit模型的具體形式如式(5)所示:

Logit(p)=b0+b0x0+…+bpxp

(5)

根據(jù)Logit變換定義,有:

將式(6)進(jìn)行運(yùn)算,可得:

其中b0為常數(shù)項(xiàng),b1,b2,…,bp為偏回歸系數(shù).

由于因變量屬于二分變量,因此運(yùn)用極大似然法進(jìn)行估計(jì).文中使用了Logit模型變換,各自變量的偏回歸系數(shù)bi(i=1,2,…,p)表示自變量xi每改變一個(gè)單位,類別劃分為高鉀玻璃和劃分為鉛鋇玻璃的發(fā)生比自然對(duì)數(shù)值變化量.

利用所給數(shù)據(jù)采用極大似然函數(shù)進(jìn)行結(jié)果溯因分析,估計(jì)出使得目前結(jié)果可能性最大的參數(shù).根據(jù)該參數(shù)能夠確定任意未知類別樣品的劃分概率和后驗(yàn)概率,最后得到極大似然估計(jì)函數(shù)［6］為:

根據(jù)玻璃類型分組,使用SPSS進(jìn)行Logit回歸,回歸準(zhǔn)確度如表2所示.

表2 Logit回歸準(zhǔn)確度表

由表2可知,回歸總體準(zhǔn)確度為84.6%,預(yù)測(cè)精度較高,因此該模型可用于預(yù)測(cè)風(fēng)化點(diǎn)風(fēng)化前的強(qiáng)相關(guān)化學(xué)成分含量.

通過該模型得到強(qiáng)相關(guān)化學(xué)成分含量預(yù)測(cè)結(jié)果,利用SPSS畫折線圖分析兩類玻璃弱相關(guān)化學(xué)成分與強(qiáng)相關(guān)化學(xué)成分之間的關(guān)系,總結(jié)如下規(guī)律:

1)高鉀玻璃和鉛鋇玻璃:氧化鈣(CaO)與二氧化硅(SiO2)、氧化鎂(MgO)與氧化鐵(Fe2O3)存在某種函數(shù)關(guān)系,剩余弱相關(guān)化學(xué)成分與氧化鉀(K2O)存在某種函數(shù)關(guān)系.

2) 兩種玻璃:氧化鉀(K2O)、氧化銅(CuO)、二氧化硫(SO2)與二氧化硅(SiO2)存在某種函數(shù)關(guān)系,氧化鋁(Al2O3)與氧化鉛(PbO)存在某種函數(shù)關(guān)系,氧化鎂(MgO)、氧化鈣(CaO)、氧化鐵(Fe2O3)與氧化鋇(BaO)存在某種函數(shù)關(guān)系.

以高鉀玻璃中氧化鈣(CaO)與二氧化硅(SiO2)為例,建立式(9)所示關(guān)系式:

y=-0.002901x3+0.5583x2-35.54x+754.3

(9)

利用python擬合后函數(shù)圖像如圖2所示,由預(yù)測(cè)所得風(fēng)化前二氧化硅含量,進(jìn)而求得氧化鈣含量.鉛鋇玻璃文物類型風(fēng)化前化學(xué)成分含量也可求得.

圖2 高鉀玻璃文物氧化鈣(CaO)與二氧化硅(SiO2)擬合圖像

3 玻璃文物的分類和亞分類規(guī)則

3.1 玻璃文物分類規(guī)則選擇

對(duì)多次取樣的文物,計(jì)算其化學(xué)成分的平均值,用于表征文物的主要化學(xué)成分.查閱相關(guān)資料［4］,將各化學(xué)成分的主要用途進(jìn)行分類,如表3所示.

通過表3數(shù)據(jù),將著色劑和含量過低的其他氧化物剔除,判斷文物化學(xué)成分中二氧化硅、氧化鋁、氧化鉀、氧化鈣、氧化鉛和氧化鋇6種氧化物為影響古玻璃文物分類的主要因素.

將所給數(shù)據(jù)分成風(fēng)化和未風(fēng)化兩部分,去除其余8項(xiàng)特征,利用python軟件通過層次聚類［7］得到如圖3所示結(jié)果.表面風(fēng)化的文物除編號(hào)2和編號(hào)48外,其他分類均正確,表面未風(fēng)化文物除編號(hào)21外,其他分類均正確,該結(jié)果準(zhǔn)確度較高.

圖3 風(fēng)化與未風(fēng)化聚類結(jié)果圖

從聚類分析樹形圖可知,如從閥值λ約為0.7處 (圖3 b)劃分,可將所有樣品劃分為兩大類.

1)古玻璃樣品屬于高鉀玻璃:PbO、BaO含量較低,K2O、CaO含量較高.K2O是主要助熔劑氧化物,CaO可能是石灰石穩(wěn)定劑生成的產(chǎn)物或是其他助熔劑氧化物.

2)古玻璃樣品屬于鉛鋇玻璃:主要特點(diǎn)是PbO、BaO含量很高,K2O、CaO含量較低.PbO、BaO是主要助熔劑氧化物.

上述兩類玻璃的Al2O3都較高,但燒制過程中形成的中間氧化物不能作為分類依據(jù).

結(jié)合層次聚類結(jié)果,觀察可知兩種玻璃的分類規(guī)律.

1)高鉀玻璃分類規(guī)律:當(dāng)文物表面氧化鉀成分占比最高時(shí)可分為高鉀玻璃類.

2)鉛鋇玻璃分類規(guī)律:當(dāng)文物表面氧化鉛和氧化鋇成分占比最高時(shí)可分為鉛鋇玻璃.

3.2 玻璃文物亞分類規(guī)則選擇

采用K-means算法［8］分別對(duì)高鉀玻璃和鉛鋇玻璃進(jìn)行亞分類聚類,具體算法如下.

Step1:在樣本數(shù)據(jù)集D中選擇k個(gè)樣本點(diǎn),將k個(gè)樣本點(diǎn)的值分別賦值給初始聚類中心:

Step5:計(jì)算數(shù)據(jù)集D中所有點(diǎn)的平方差Ei,并且與前一次誤差Ei-1作比較:

判斷|Ei-1-Ei|<δ是否成立,若成立則算法結(jié)束,否則轉(zhuǎn)入Step2進(jìn)行二次迭代,直到|Ei-1-Ei|<δ,成立.

高鉀玻璃選取氧化鈉(Na2O)、氧化鈣(CaO)、氧化鎂(MgO)、氧化鋁(Al2O3)和五氧化二磷(P2O5)成分進(jìn)行分析,二氧化硅(SiO2)、氧化鉀(K2O)是高鉀玻璃的主要成分,氧化鐵(Fe2O3)、氧化銅(CuO)是著色劑,氧化鉛(PbO)、氧化鋇(BaO)、氧化鍶(SrO)、氧化錫(SnO2)、二氧化硫(SO2)等含量極低,可忽略不計(jì).

鉛鋇玻璃選取氧化鈣(CaO)、氧化鋁(Al2O3)成分進(jìn)行分析(此處刪去了另外三種化學(xué)成分,更有利于類別命名).

利用python求解模型,結(jié)合手肘法和輪廓系數(shù)法可知,高鉀玻璃分成四類較合理,鉛鋇玻璃分成三類較合理.

高鉀玻璃聚類結(jié)果如表4所示,根據(jù)高鉀玻璃中CaO 和 Al2O3的濃度以及是否含鈉元素［2］,分別將0、1、2、3命名為含鈉中鋁中鈣高鉀玻璃、低鈣低鋁高鉀玻璃、低鈣高鋁高鉀玻璃、不含鈉中鋁中鈣高鉀玻璃.

表4 高鉀玻璃分類結(jié)果

鉛鋇玻璃聚類結(jié)果如表5所示.根據(jù)鉛鋇玻璃中 CaO 和 Al2O3濃度,將0、1、2分別命名為低鈣低鋁鉛鋇玻璃、高鈣中鋁鉛鋇玻璃、中鈣高鋁鉛鋇玻璃.

表5 鉛鋇玻璃分類結(jié)果

3.3 合理性分析

表6 高鉀方差分析

表7 鉛鋇方差分析

從表6、表7可以看出,不同類別樣本選取的化學(xué)成分全部呈現(xiàn)出顯著性(P<0.05),方差分析對(duì)參與聚類分析的變量能很好地區(qū)分類別,類間差異足夠大.

3.4 玻璃文物分類模型和亞分類模型的建立與求解

隨機(jī)森林預(yù)測(cè)算法流程:基于bootsrap方法重采樣,隨機(jī)生成N個(gè)訓(xùn)練集,產(chǎn)生的訓(xùn)練集與原始訓(xùn)練集P樣本總數(shù)相等.利用訓(xùn)練集構(gòu)建相對(duì)應(yīng)的決策樹K1,K2,…,KN.在每一個(gè)內(nèi)部節(jié)點(diǎn)選擇分裂特征前,從訓(xùn)練集中的M個(gè)特征中隨機(jī)選取m(m≤M)個(gè)作為當(dāng)前內(nèi)部節(jié)點(diǎn)的分裂特征集,選擇對(duì)應(yīng)基尼值最小的特征作為分裂特征.每棵樹不進(jìn)行剪枝操作,都生長到底.對(duì)于測(cè)試集數(shù)據(jù)X,每個(gè)決策樹給出獨(dú)立的預(yù)測(cè)結(jié)果K1(X),K2(X),…,KN(X),即投出自己的一票.統(tǒng)計(jì)每一個(gè)決策樹的預(yù)測(cè)結(jié)果,將所得票數(shù)最多的預(yù)測(cè)值作為最終的預(yù)測(cè)結(jié)果［9］.

根據(jù)制定的分類規(guī)則,構(gòu)建兩個(gè)隨機(jī)森林分類器,將待分類數(shù)據(jù)分成風(fēng)化和未風(fēng)化,并選取相應(yīng)化學(xué)成分,分類結(jié)果如表8.

表8 分類結(jié)果(部分)

4 結(jié)語

不同文物因組成材料不同、埋藏地點(diǎn)不同,受到的侵蝕風(fēng)化程度也不相同.為了較好地修復(fù)文物,對(duì)古代玻璃制品成分進(jìn)行分析,可依靠數(shù)學(xué)分析工具分析出文物表面風(fēng)化情況與客觀因素之間的關(guān)系.文化風(fēng)化后,主成分含量均有所降低.根據(jù)檢測(cè)成分對(duì)未知玻璃制品進(jìn)行鑒別,結(jié)合k-means聚類算法建立隨機(jī)森林分類器,推導(dǎo)出文物本身具有的化學(xué)元素,對(duì)文物修復(fù)保護(hù)著指導(dǎo)作用,同時(shí)也對(duì)研究我國古代玻璃制品的風(fēng)化條件與玻璃制品表面特征和埋藏環(huán)境之間的關(guān)系具有積極意義.

本文研究成果獲全國大學(xué)生數(shù)學(xué)建模大賽甘肅賽區(qū)一等獎(jiǎng).