基于多元線性回歸的古代玻璃體系判別方法

錢鴻羽，岳思情，許 雷，鄭心怡

（內(nèi)江師范學(xué)院，四川 內(nèi)江 641000）

玻璃作為人類最早期發(fā)明的人造材料之一[1]，在世界各國之間的貿(mào)易與交流活動(dòng)中占據(jù)重要位置。由于其相對(duì)穩(wěn)定的化學(xué)性質(zhì)，迄今為止保存了大量古代玻璃制品，這是研究國家間早期貿(mào)易的寶貴物證。對(duì)古代玻璃制品的探究，對(duì)進(jìn)一步研究各國間的交流活動(dòng)具有重要意義。當(dāng)前，學(xué)術(shù)界對(duì)古玻璃的研究主要集中在古代玻璃的成分體系研究、玻璃起源與產(chǎn)地研究、玻璃制作工藝研究等領(lǐng)域[2]，而對(duì)古代玻璃化學(xué)成分體系的研究為進(jìn)一步研究玻璃生產(chǎn)技術(shù)及其起源奠定了基礎(chǔ)，因而成為古代玻璃研究的重要一環(huán)。目前研究古代玻璃化學(xué)成分的方法主要有張欣睿[3]采用的X射線熒光分析、李玲等[4]采用的PIXE 技術(shù)、李穎等[5]運(yùn)用拉曼光譜、紅外光譜和ED-XRF 等技術(shù)，以及李曼[6]運(yùn)用激光剝蝕電感耦合等離子體發(fā)射光譜進(jìn)行化學(xué)成分分析的方法。由于古代玻璃易受外界影響而產(chǎn)生風(fēng)化，從而導(dǎo)致內(nèi)部化學(xué)成分發(fā)生變化，影響對(duì)古代玻璃化學(xué)成分體系的判別。在2022 年全國大學(xué)生數(shù)學(xué)建模競賽C題的基礎(chǔ)上，引入風(fēng)化前玻璃各化學(xué)成分預(yù)測函數(shù)，通過風(fēng)化后古代玻璃的化學(xué)組成含量預(yù)測相應(yīng)的風(fēng)化前化學(xué)成分，從而確定玻璃樣品的成分。

1 理論基礎(chǔ)

1.1 加性對(duì)數(shù)比變換（ALR）

加性對(duì)數(shù)比變換（ALR）是目前在成分分析中運(yùn)用較為普遍的數(shù)據(jù)轉(zhuǎn)換方法。為了消除定和約束的影響，通過計(jì)算各個(gè)組別所占成分比例與其均值的比值的對(duì)數(shù)，將成分?jǐn)?shù)據(jù)從單形空間變換到歐式空間，從而能夠應(yīng)用經(jīng)典統(tǒng)計(jì)方法解決成分?jǐn)?shù)據(jù)問題。

其具體計(jì)算公式如下：

式中：i表示數(shù)據(jù)的組別，xi表示各組數(shù)據(jù)的成分比例，xp表示數(shù)據(jù)比例的平均數(shù)，yi則為新產(chǎn)生的數(shù)據(jù)變量。

1.2 獨(dú)立樣本Mann-Whitney U檢驗(yàn)

Mann-Whitney U檢驗(yàn)是用得最廣泛的兩獨(dú)立樣本秩和檢驗(yàn)方法，其假設(shè)基礎(chǔ)是：若兩個(gè)樣本有差異，則它們的中心位置將不同。具體步驟如圖1 所示。

圖1 獨(dú)立樣本M ann-W hitneyU 檢驗(yàn)流程圖

1.3 最小二乘法

定義2 個(gè)變量的n對(duì)樣本數(shù)據(jù)為（x1,y1），（x2,y2），…，（xn,yn）設(shè)回歸直線為：y=a+bx。

由yi=bxi+a+ei得代表樣本點(diǎn)（xi,yi）到直線上點(diǎn)（xi,bxi+a）的距離。

當(dāng) |ei|=0 時(shí)，則表示點(diǎn)（xi,yi）在直線上，利用來刻畫各樣本數(shù)據(jù)與直線y=a+bx的“整體接近程度”；為了便于計(jì)算，通常利用來刻畫“整體接近程度”。

2 風(fēng)化前化學(xué)成分預(yù)測模型

由于古代玻璃易受風(fēng)化作用影響而改變其各化學(xué)成分含量，因此根據(jù)樣本數(shù)據(jù)建立以下風(fēng)化前化學(xué)成分預(yù)測模型。

2.1 選取重要性指標(biāo)

根據(jù)干福熹研究所得，我國早期古代玻璃體系可分為鉛鋇玻璃、高鉀玻璃等[7]，故可按玻璃體系分別建立玻璃風(fēng)化程度函數(shù)與風(fēng)化前化學(xué)成分預(yù)測函數(shù)。

定義玻璃的風(fēng)化程度函數(shù)為：

其中a為常量，Ci為重要性指標(biāo)在玻璃制品中的含量比例，bi為重要性指標(biāo)的對(duì)應(yīng)系數(shù)，n為重要性指標(biāo)的數(shù)量。

將樣本數(shù)據(jù)進(jìn)行預(yù)處理，將無法檢測出化學(xué)成分的部分填充為0，去除成分比例累加總和處于85%～105%區(qū)間范圍外的數(shù)據(jù)，以便于進(jìn)行后續(xù)數(shù)據(jù)分析。

由于n元成分?jǐn)?shù)據(jù)受到定和約束，故必然存在一組數(shù)據(jù)與其余數(shù)據(jù)呈現(xiàn)負(fù)相關(guān)情況，易與實(shí)際情況產(chǎn)生偏差，故針對(duì)普通數(shù)據(jù)的傳統(tǒng)統(tǒng)計(jì)學(xué)分析方法難以適用于成分?jǐn)?shù)據(jù)分析。為了進(jìn)一步對(duì)數(shù)據(jù)進(jìn)行分析，對(duì)數(shù)據(jù)進(jìn)行加性對(duì)數(shù)比變換（ALR）處理，將原始數(shù)據(jù)比例關(guān)系轉(zhuǎn)化為加性關(guān)系，從而可以更加充分的體現(xiàn)成分特性，使成分?jǐn)?shù)據(jù)中的可解釋性更強(qiáng)。

根據(jù)玻璃體系劃分樣品數(shù)據(jù)，提取處理后數(shù)據(jù)中玻璃體系屬于高鉀的樣本，并對(duì)化學(xué)成分含量進(jìn)行正態(tài)性檢驗(yàn)。發(fā)現(xiàn)總體不符合正態(tài)分布。因此，選擇非參數(shù)檢驗(yàn)中的獨(dú)立樣本Mann-Whitney U 檢驗(yàn)[8]，按照Cohen′s d 系數(shù)排序分析具有顯著性差異的化學(xué)成分，從而得到文物樣品表面是否存在風(fēng)化與化學(xué)成分含量之間的統(tǒng)計(jì)規(guī)律，確定了判斷玻璃風(fēng)化程度的重要指標(biāo)，檢驗(yàn)結(jié)果見表1。

表1 M ann-W hitneyU 檢驗(yàn)分析結(jié)果表

2.2 建立多元線性回歸模型

對(duì)Cohen′s d系數(shù)進(jìn)行排序，選取排名靠前的部分化學(xué)成分，認(rèn)為它們有顯著差異，并將其作為函數(shù)的自變量。高鉀玻璃風(fēng)化程度函數(shù)的自變量為玻璃中二氧化硅、氧化鉀、氧化鈣、氧化鋁的對(duì)應(yīng)含量比例。

根據(jù)玻璃樣本大體的風(fēng)化情況對(duì)風(fēng)化程度進(jìn)行賦值。根據(jù)多次擬合結(jié)果，確定未風(fēng)化樣品的值為0，一般風(fēng)化樣品賦值為0.8。使用最小二乘法進(jìn)行函數(shù)擬合，以確定相應(yīng)化學(xué)成分的系數(shù)與未知常數(shù)。對(duì)獲得的殘差進(jìn)行檢驗(yàn)，去除異常值，并進(jìn)行多次線性擬合，直至總體殘差的置信區(qū)間均包含零點(diǎn)，表明回歸模型可以更好地符合原始數(shù)據(jù)，從而得到最優(yōu)的多元線性回歸函數(shù)。

擬合所得的高鉀玻璃風(fēng)化程度函數(shù)為：

其中YK代表高鉀玻璃的風(fēng)化比例，CSiO2代表SiO2在玻璃制品中的含量比例，CK2O，CCaO，CAl2O3同理。函數(shù)擬合程度如圖2所示。

圖2 高鉀玻璃風(fēng)化程度函數(shù)擬合結(jié)果圖

同理，其余體系的玻璃風(fēng)化程度函數(shù)也可按照上述方法得出。

通過對(duì)鉛鋇玻璃樣本中的數(shù)據(jù)進(jìn)行獨(dú)立樣本Mann-Whitney U 檢驗(yàn)，按照Cohen′s d 系數(shù)進(jìn)行排序，可得鉛鋇玻璃風(fēng)化程度函數(shù)的自變量為玻璃中二氧化硅、氧化鈣、氧化鋁、氧化鉛、五氧化二磷的對(duì)應(yīng)含量比例。根據(jù)多次擬合結(jié)果，確定將未風(fēng)化樣本賦值為0，一般風(fēng)化樣品為0.75，嚴(yán)重風(fēng)化樣品賦值為1。經(jīng)過多次擬合與剔除異常值，擬合得到的鉛鋇玻璃風(fēng)化程度函數(shù)為：

其中YPb代表高鉀玻璃的風(fēng)化比例，代表SiO2在玻璃制品中的含量比例同理。函數(shù)擬合程度如圖3所示。

圖3 鉛鋇玻璃風(fēng)化程度函數(shù)擬合結(jié)果圖

計(jì)算函數(shù)預(yù)測結(jié)果與實(shí)際結(jié)果的差值，得到預(yù)測數(shù)據(jù)的偏差情況，所得結(jié)果如圖4所示。

圖4 樣本預(yù)測數(shù)據(jù)偏差情況圖

2.3 風(fēng)化前化學(xué)成分預(yù)測模型

處理后的樣品數(shù)據(jù)應(yīng)根據(jù)玻璃系統(tǒng)進(jìn)行分類，并根據(jù)風(fēng)化情況進(jìn)行二次分類。去除異常數(shù)據(jù)后，計(jì)算不同風(fēng)化情況下各化學(xué)成分自身的差值。由于在化學(xué)反應(yīng)中涉及的化學(xué)成分是成比例的，因此可以通過簡化得到不同玻璃體系下風(fēng)化過程中各化學(xué)成分反應(yīng)比例。高鉀玻璃與鉛鋇玻璃風(fēng)化反應(yīng)比例分別為：

分別探究玻璃在不同玻璃體系下的風(fēng)化程度與變化程度最大的化學(xué)成分含量之間的線性關(guān)系，得到高鉀玻璃和鉛鋇玻璃的風(fēng)化程度與二氧化硅含量的函數(shù)分別為：

通過反解該線性函數(shù)，可以得到相應(yīng)風(fēng)化程度下參與反應(yīng)的二氧化硅含量。將獲得的數(shù)據(jù)乘以各化學(xué)成分反應(yīng)比例，即可得對(duì)應(yīng)樣本在風(fēng)化過程中參與反應(yīng)的各化學(xué)成分含量。風(fēng)化前的數(shù)據(jù)可以通過從風(fēng)化后的數(shù)據(jù)中減去變化來獲得。

高鉀玻璃與鉛鋇玻璃各化學(xué)成分含量計(jì)算公式分別為：

其中Ci為風(fēng)化后對(duì)應(yīng)化學(xué)成分含量，為所求風(fēng)化前對(duì)應(yīng)化學(xué)成分含量，qi為對(duì)應(yīng)反應(yīng)比例。

將預(yù)測函數(shù)所得的數(shù)據(jù)代入擬合的風(fēng)化程度函數(shù)，發(fā)現(xiàn)5.1%以內(nèi)的計(jì)算偏差比例為94.12%，其中嚴(yán)重風(fēng)化情況下的樣本計(jì)算偏差相對(duì)較大，故認(rèn)為預(yù)測函數(shù)較為有效。

3 玻璃化學(xué)成分體系判別模型

3.1 選取重要性指標(biāo)

為避免風(fēng)化引起的化學(xué)成分含量變化對(duì)判定玻璃體系的影響，本次從樣本數(shù)據(jù)中選取無風(fēng)化樣本作為數(shù)據(jù)擬合的內(nèi)容。篩選后，對(duì)類別進(jìn)行賦值處理，將“-1”賦值給高鉀玻璃，將“1”賦值給鉛鋇玻璃，使得定類變量“類別”能夠轉(zhuǎn)化為定量變量，便于后續(xù)判斷和篩選對(duì)玻璃體系影響較大的化學(xué)成分。

通過對(duì)化學(xué)成分進(jìn)行正態(tài)性檢驗(yàn)，發(fā)現(xiàn)大部分?jǐn)?shù)據(jù)并不滿足正態(tài)分布。故同樣采用獨(dú)立樣本Mann-Whitney U 檢驗(yàn)，計(jì)算各個(gè)化學(xué)成分的Cohen′s d 系數(shù)并對(duì)其進(jìn)行排名，所得結(jié)果見表2。

表2 M ann-W hitneyU 檢驗(yàn)分析結(jié)果表

3.2 建立玻璃體系判別模型

選取Cohen′s d 系數(shù)大于1 的化學(xué)成分，這些化學(xué)成分被認(rèn)為具有較大顯著性差異。因此，可以得出結(jié)論，對(duì)玻璃體系的判別能力具有顯著性差異從大到小分別為氧化鉀、氧化鉛、氧化鈣、氧化鋇，氧化硅和氧化鍶，它們?cè)诓Ａе邢鄳?yīng)化學(xué)成分的比例即為判別函數(shù)的自變量。

將量化后的玻璃體系作為函數(shù)因變量，用最小二乘法求解相應(yīng)自變量的系數(shù)。經(jīng)過多次擬合與異常值處理，可獲得玻璃體系判別函數(shù)為：

函數(shù)擬合程度如圖5所示。

圖5 玻璃體系判別函數(shù)擬合結(jié)果圖

選取部分樣本數(shù)據(jù)代入判別函數(shù)，經(jīng)過多次測試，發(fā)現(xiàn)當(dāng)判斷偏差設(shè)置為±0.4時(shí)，該測試的準(zhǔn)確度約為93.75%，因此認(rèn)為多元線性回歸所得的判別函數(shù)具有較高的準(zhǔn)確性。

將預(yù)測結(jié)果與實(shí)際數(shù)據(jù)的偏差情況可視化，其結(jié)果如圖6所示。

圖6 判別函數(shù)預(yù)測結(jié)果偏差圖

4 玻璃類型判別方法

古代玻璃易受到外界影響而風(fēng)化，從而導(dǎo)致內(nèi)部化學(xué)成分含量發(fā)生較大變化，故難以直接通過玻璃化學(xué)成分體系判別模型得出結(jié)果，所以首先對(duì)該玻璃樣本風(fēng)化前的化學(xué)成分進(jìn)行預(yù)測。根據(jù)探究得到高鉀玻璃和鉛鋇玻璃的風(fēng)化程度與二氧化硅含量的關(guān)系分別如公式（5）、（6）所示，即鉛鋇玻璃與高鉀玻璃在風(fēng)化后其二氧化硅和其他化學(xué)成分含量的變化趨勢相反，故若高鉀玻璃被判別為鉛鋇玻璃代入其風(fēng)化前化學(xué)成分預(yù)測模型時(shí)，會(huì)與其風(fēng)化前數(shù)據(jù)產(chǎn)生更大的差異。此時(shí)代入玻璃體系判別后所得結(jié)果必然與假設(shè)相矛盾。綜上可得，當(dāng)且僅當(dāng)對(duì)玻璃類型的假設(shè)正確時(shí)，代入玻璃體系判別模型后所得的結(jié)果才與假設(shè)相一致。

因此，可以先對(duì)玻璃樣品類型進(jìn)行假設(shè)。首先，假設(shè)所有風(fēng)化玻璃均為高鉀玻璃。將數(shù)據(jù)代入高鉀玻璃風(fēng)化前后化學(xué)元素含量的變化函數(shù)，計(jì)算風(fēng)化前各化學(xué)成分的含量。將獲得的數(shù)據(jù)代入玻璃體系判別函數(shù)，判斷所得結(jié)果是否在可接受范圍內(nèi)。如果是，則認(rèn)為該假設(shè)成立，文物屬于高鉀玻璃；如果不是，則假設(shè)與結(jié)果相矛盾，假設(shè)不成立，再次假設(shè)該文物屬于鉛鋇玻璃，重復(fù)上述步驟。

若兩者皆不成立，則將該玻璃樣品進(jìn)行單獨(dú)討論，計(jì)算結(jié)果與判定標(biāo)準(zhǔn)之間的差距，以差距較小的一方作為判別結(jié)果。

5 實(shí)例計(jì)算

基于2022 年全國大學(xué)生數(shù)學(xué)建模競賽C 題所提供的數(shù)據(jù)，對(duì)表單中未知體系的風(fēng)化玻璃文物的化學(xué)成分進(jìn)行分析，鑒別其所屬體系，具體數(shù)據(jù)見表3。

表3 玻璃樣品及其化學(xué)成分含量表

將提供的玻璃樣本均假設(shè)為高鉀玻璃，并代入對(duì)應(yīng)的風(fēng)化前各化學(xué)成分預(yù)測函數(shù)，得到風(fēng)化前各化學(xué)成分的含量。將所得數(shù)據(jù)代入玻璃體系判別函數(shù)，觀察是否與假設(shè)相一致。若不一致，則假設(shè)樣本的體例為鉛鋇玻璃并重復(fù)以上步驟。所得數(shù)據(jù)與判斷結(jié)果見表4。

表4 風(fēng)化前化學(xué)成分含量預(yù)測與分類結(jié)果表

從表4中觀察所得，發(fā)現(xiàn)編號(hào)為08與26的兩份樣品在兩種假設(shè)下的偏差均高于設(shè)定標(biāo)準(zhǔn)，通過計(jì)算并比較其在不同玻璃體系情況下的偏差程度，并取較小一方作為判斷結(jié)果，可得編號(hào)08 與26 的玻璃體系為鉛鋇玻璃。

6 總結(jié)

研究玻璃制品的化學(xué)成分及鑒別對(duì)保護(hù)古代玻璃文物具有重要意義。本文從成分?jǐn)?shù)據(jù)視角出發(fā)，采用加性對(duì)數(shù)比變換以消除定和限制對(duì)后續(xù)分析的影響。引入了風(fēng)化前玻璃各化學(xué)成分預(yù)測函數(shù)，并通過風(fēng)化后古玻璃的化學(xué)組成含量預(yù)測相應(yīng)的風(fēng)化前化學(xué)成分含量，從而確定玻璃樣品的成分體系，得到古代玻璃體系判別模型。在應(yīng)用本模型時(shí)，先對(duì)玻璃樣品的類別進(jìn)行假設(shè)，在假設(shè)的基礎(chǔ)上代入相應(yīng)的風(fēng)化前化學(xué)成分預(yù)測模型，得到在該假設(shè)下風(fēng)化前各化學(xué)成分比例。隨后將此數(shù)據(jù)代入古代玻璃體系判別模型，根據(jù)所得結(jié)果與假設(shè)的一致性確定玻璃樣品所屬的玻璃類別。模型擬合過程中所運(yùn)用的多元線性回歸方法降低了計(jì)算的復(fù)雜性，同時(shí)確保了結(jié)果的準(zhǔn)確性。最后，通過實(shí)例說明了該方法的有效性和適用性，為古代玻璃化學(xué)成分體系的分析與判別提供了新的思路。