杜世強,樂靖雯,周海麗,楊松杰

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

0 引言

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

利用卡方檢驗和相關性分析玻璃文物表面風化和類型、紋飾、顏色的關系以及所含化學成分之間的相關性.結合文獻[4],對化學成分的主要用途進行分類,再根據(jù)主要化學成分建立層次聚類和K-means聚類模型,得出玻璃文物分類、亞分類標準,最后構建隨機森林分類器模型,進而對未知類型玻璃進行分類.

1 相關因素分析

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

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

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

1.2 玻璃文物表面風化與其類型、紋飾和顏色的相關性分析

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

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

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

表1 表面風化與紋飾、玻璃類型、顏色相關性卡方檢驗表

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

1.3 玻璃文物中化學成分間的相關性分析模型建立與求解

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

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

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

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

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

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

2 玻璃文物中化學成分預測

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

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

2.1 風化前化學成分預測模型的建立與求解

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

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

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

(5)

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

將式(6)進行運算,可得:

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

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

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

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

表2 Logit回歸準確度表

由表2可知,回歸總體準確度為84.6%,預測精度較高,因此該模型可用于預測風化點風化前的強相關化學成分含量.

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

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

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

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

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

(9)

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

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

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

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

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

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

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

圖3 風化與未風化聚類結果圖

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

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

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

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

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

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

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

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

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

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

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

判斷|Ei-1-Ei|<δ是否成立,若成立則算法結束,否則轉入Step2進行二次迭代,直到|Ei-1-Ei|<δ,成立.

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

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

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

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

表4 高鉀玻璃分類結果

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

表5 鉛鋇玻璃分類結果

3.3 合理性分析

表6 高鉀方差分析

表7 鉛鋇方差分析

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

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

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

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

表8 分類結果(部分)

4 結語

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

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