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

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

（內江師范學院，四川 內江 641000）

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

1 理論基礎

1.1 加性對數比變換（ALR）

加性對數比變換（ALR）是目前在成分分析中運用較為普遍的數據轉換方法。為了消除定和約束的影響，通過計算各個組別所占成分比例與其均值的比值的對數，將成分數據從單形空間變換到歐式空間，從而能夠應用經典統計方法解決成分數據問題。

其具體計算公式如下：

式中：i表示數據的組別，xi表示各組數據的成分比例，xp表示數據比例的平均數，yi則為新產生的數據變量。

1.2 獨立樣本Mann-Whitney U檢驗

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

圖1 獨立樣本M ann-W hitneyU 檢驗流程圖

1.3 最小二乘法

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

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

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

2 風化前化學成分預測模型

由于古代玻璃易受風化作用影響而改變其各化學成分含量，因此根據樣本數據建立以下風化前化學成分預測模型。

2.1 選取重要性指標

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

定義玻璃的風化程度函數為：

其中a為常量，Ci為重要性指標在玻璃制品中的含量比例，bi為重要性指標的對應系數，n為重要性指標的數量。

將樣本數據進行預處理，將無法檢測出化學成分的部分填充為0，去除成分比例累加總和處于85%～105%區間范圍外的數據，以便于進行后續數據分析。

由于n元成分數據受到定和約束，故必然存在一組數據與其余數據呈現負相關情況，易與實際情況產生偏差，故針對普通數據的傳統統計學分析方法難以適用于成分數據分析。為了進一步對數據進行分析，對數據進行加性對數比變換（ALR）處理，將原始數據比例關系轉化為加性關系，從而可以更加充分的體現成分特性，使成分數據中的可解釋性更強。

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

表1 M ann-W hitneyU 檢驗分析結果表

2.2 建立多元線性回歸模型

對Cohen′s d系數進行排序，選取排名靠前的部分化學成分，認為它們有顯著差異，并將其作為函數的自變量。高鉀玻璃風化程度函數的自變量為玻璃中二氧化硅、氧化鉀、氧化鈣、氧化鋁的對應含量比例。

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

擬合所得的高鉀玻璃風化程度函數為：

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

圖2 高鉀玻璃風化程度函數擬合結果圖

同理，其余體系的玻璃風化程度函數也可按照上述方法得出。

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

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

圖3 鉛鋇玻璃風化程度函數擬合結果圖

計算函數預測結果與實際結果的差值，得到預測數據的偏差情況，所得結果如圖4所示。

圖4 樣本預測數據偏差情況圖

2.3 風化前化學成分預測模型

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

分別探究玻璃在不同玻璃體系下的風化程度與變化程度最大的化學成分含量之間的線性關系，得到高鉀玻璃和鉛鋇玻璃的風化程度與二氧化硅含量的函數分別為：

通過反解該線性函數，可以得到相應風化程度下參與反應的二氧化硅含量。將獲得的數據乘以各化學成分反應比例，即可得對應樣本在風化過程中參與反應的各化學成分含量。風化前的數據可以通過從風化后的數據中減去變化來獲得。

高鉀玻璃與鉛鋇玻璃各化學成分含量計算公式分別為：

其中Ci為風化后對應化學成分含量，為所求風化前對應化學成分含量，qi為對應反應比例。

將預測函數所得的數據代入擬合的風化程度函數，發現5.1%以內的計算偏差比例為94.12%，其中嚴重風化情況下的樣本計算偏差相對較大，故認為預測函數較為有效。

3 玻璃化學成分體系判別模型

3.1 選取重要性指標

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

通過對化學成分進行正態性檢驗，發現大部分數據并不滿足正態分布。故同樣采用獨立樣本Mann-Whitney U 檢驗，計算各個化學成分的Cohen′s d 系數并對其進行排名，所得結果見表2。

表2 M ann-W hitneyU 檢驗分析結果表

3.2 建立玻璃體系判別模型

選取Cohen′s d 系數大于1 的化學成分，這些化學成分被認為具有較大顯著性差異。因此，可以得出結論，對玻璃體系的判別能力具有顯著性差異從大到小分別為氧化鉀、氧化鉛、氧化鈣、氧化鋇，氧化硅和氧化鍶，它們在玻璃中相應化學成分的比例即為判別函數的自變量。

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

函數擬合程度如圖5所示。

圖5 玻璃體系判別函數擬合結果圖

選取部分樣本數據代入判別函數，經過多次測試，發現當判斷偏差設置為±0.4時，該測試的準確度約為93.75%，因此認為多元線性回歸所得的判別函數具有較高的準確性。

將預測結果與實際數據的偏差情況可視化，其結果如圖6所示。

圖6 判別函數預測結果偏差圖

4 玻璃類型判別方法

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

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

若兩者皆不成立，則將該玻璃樣品進行單獨討論，計算結果與判定標準之間的差距，以差距較小的一方作為判別結果。

5 實例計算

基于2022 年全國大學生數學建模競賽C 題所提供的數據，對表單中未知體系的風化玻璃文物的化學成分進行分析，鑒別其所屬體系，具體數據見表3。

表3 玻璃樣品及其化學成分含量表

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

表4 風化前化學成分含量預測與分類結果表

從表4中觀察所得，發現編號為08與26的兩份樣品在兩種假設下的偏差均高于設定標準，通過計算并比較其在不同玻璃體系情況下的偏差程度，并取較小一方作為判斷結果，可得編號08 與26 的玻璃體系為鉛鋇玻璃。

6 總結

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