基于多重二階差商濾波的礦物油檢測自動判別方法

關鍵詞 礦物油測定；液相色譜-氣相色譜聯用法；鼓包峰輪廓；二階差商；濾波

由飽和烴（MOSH）和芳香烴（MOAH）組成的礦物油污染物在食用油和相關加工食品中普遍存在[1-3]。污染源可能來自食品包裝材料，也可能是食用油生產、運輸過程中工業油殘存導致的污染[4-8]。毒理學研究結果表明， MOSH 容易蓄積， MOAH 具有致癌和遺傳毒性。近年來，歐盟對礦物油的監管日趨嚴格。2022 年，歐盟植物、動物、食品和飼料委員會SCoPAFF 發布了《新型食品和食物鏈的毒理學安全-2022/04/21》[9]，規定了不同脂肪含量食品中MOAH 的限量，即低脂肪含量（≤4%脂肪）的干性食品中MOAH 限量為0.5 mg/kg，較高脂肪含量（4%～50%脂肪）的食品中MOAH 限量為1.00 mg/kg，脂肪和油脂（≥50%脂肪）中MOAH 限量為2 mg/kg[9]。我國國家標準GB 2760—2024《食品添加劑使用標準》允許食品級礦物油（白油和液體石蠟）作為加工助劑使用，在糖果和鮮蛋中的最大使用量為5.0 g/kg[10]。然而，GB/T 21309—2007《涂漬油脂或石蠟大米檢驗》和GB/T 37514—2019《動植物油油脂礦物油的檢測》兩項檢測標準對礦物油進行定性判斷的檢出限為百分級別，遠不能滿足歐盟限量標準要求[11-12]。

目前，礦物油定量分析主要采用液相-氣相色譜聯用技術（LC-GC-FID）[13-17]，分析得到的MOSH 和MOAH 譜圖均呈不規則的駝峰形狀。根據FID 的等碳響應，在駝峰外添加非礦物油組分的烴類-內標物，利用其峰面積即可實現礦物油的定量分析[18]。駝峰上方的尖峰或窄峰為天然烴類干擾物，通常不會影響駝峰輪廓的判別[1]。但是，當干擾物凈化不完全時，會在駝峰上方形成一個或數個豁口，干擾駝峰輪廓的準確判斷與勾勒，由于人們對駝峰的輪廓判斷存在主觀差異，導致礦物油分析結果的重現性較差。Bauwens 等[19]提出了一種去除窄尖峰、保留寬駝峰的方法，在每個二維（2D）點重建一維（1D）色譜圖（類似于1D 軸上的投影），并根據設置的“基線寬度”消除窄峰。在保留寬駝峰的同時，作為駝峰頂部的峰值出現的信號會自動移除。但是，該方法可實現的功能較少，而且對某些較寬的尖峰去除效果不理想。

本研究針對現有問題設計了一套基于多重二階差商濾波的駝峰曲線形態的自動判別算法。由于尖峰的下極值點高于駝峰曲線，而駝峰曲線較平滑，相對駝峰上的點，尖峰所處的下極值點的二階差商通常為絕對值較大的負數。因此，通過過濾二階差商較大的點可以達到去除尖峰的目的。此外，由于豁口是由尖峰簇組成的，通過多次過濾二階差商較大的點，還可以去除豁口。結果表明，對駝峰曲線進行正向、反向、正向三次二階差商過濾，即可去除色譜數據中的豁口。根據上述原理設計了相應的礦物油自動分析軟件，在使用軟件處理駝峰數據時，用戶可以設置不同的過濾“閾值”，也可在自動形成的駝峰折線上手動刪除不必要的數據點。在進行積分運算時，用戶可以選擇使用“自動基線”或“基線數據”，也可以選擇“按駝峰積分”或“碳譜積分”以確定不同的積分區間。經實驗驗證，本算法及軟件可以有效地處理礦物油檢測的色譜數據，在保證計算精度的基礎上大幅提高計算效率，縮短數據處理的時間。

1 實驗部分

1.1 儀器與試劑

液相色譜-氣相色譜聯用儀（LC-GC，日本島津公司）：配備二元泵和紫外檢測器的LC20A 液相色譜儀和配備氫火焰離子化檢測器（FID）的GC 2010 plus 氣相色譜儀；LC-GC 接口（德國Axel Semrau 公司）和PAL 自動進樣器（瑞士CTC 公司）；Biofuge primo R 低溫高速離心機（美國賽默飛世爾科技公司）；R-100 旋轉蒸發儀（瑞士步琦有限公司）；TUS-200P 恒溫金屬浴（上海一恒科技有限公司）。

正己烷、二氯甲烷和無水乙醇（農殘級，德國默克公司）；異丙醇（色譜級，美國賽默飛世爾科技公司）；3-氯過氧苯甲酸（m-CPBA， 70%，日本東京化學工業公司）；硫代硫酸鈉（分析純，國藥集團化學試劑有限公司）。MOSH/MOAH 標準溶液（北京壇墨質檢科技有限公司），包括：正十一烷（n-C11）、正十三烷（n-C1₃）、環己基環已烷（Cycy）、5α-膽甾烷（Cho）、正戊基苯（5B）、1，3，5-三叔丁基苯（TBB）、1-甲基萘（1-MN）、2-甲基萘（2-MN）和苝（Per）。其中， n-C11、Cycy、5B、TBB、1-MN、2-MN 的濃度為300 mg/L， Cho 和Per 的濃度為600 mg/L， n-C₁₃ 的濃度為150mg/L， Cycy 為MOSH 定量內標， TBB 或2-MN 為MOAH 定量內標，其余內標和C₇～C₄₀的正構烷烴混合物（德國默克公司）用于評估LC-GC-FID 儀器性能[13]。硅膠60（0.063～0.200 mm，德國默克公司）使用前在400 ℃活化16 h，冷卻后置于干燥器內，7 d 內使用。BAM-K009 潤滑油標準品（德國柏林聯邦材料研究與測試研究所）。

1.2 樣品

2023年在電商平臺購買了42 個食品樣品，包括32 個植物油樣品、6 個植物油基魚罐頭樣品、2 個嬰幼兒配方奶粉樣品和2 個奶酪樣品。所有樣品在保質期內進行分析。參考樣品RVEP233004 椰子油和RVEP200612 奶粉從德國標物局（DRRR）獲得；參考樣品葵花籽油和橄欖油為法國油脂研究所（ITERG）在ISO 20122 標準國際驗證試驗中提供的樣品。

1.3 實驗方法

1.3.1皂化和萃取

根據ISO 20122 標準對食用油樣品進行前處理[20]。稱取3 g 食用油樣品至100 mL 離心管中，加入20 μL MOSH/MOAH 標準溶液和30mL 正己烷-乙醇（1∶1， V/V）混合溶劑，充分混合后移取10 mL 至另一個離心管中，加入3mL KOH 溶液（33%， m/V）， 60 ℃恒溫振蕩皂化30 min 直至溶液澄清；冷卻至室溫后，向皂化液中加入5 mL 正己烷和5 mL 乙醇-水（1∶1， V/V）混合溶劑，振蕩，離心10 min；將正己烷層轉移至梨形燒瓶中，再次用5 mL 正己烷萃取殘留物；合并2 次正己烷提取液，轉移至3 g 硅膠柱，用15 mL正己烷-二氯甲烷（7∶3， V/V）混合溶劑洗脫；收集所有洗脫液，濃縮至1～2 mL。

參照文獻[21]的方法進行奶粉前處理。稱取5 g 樣品至離心管，加入10 μL MOSH/MOAH 標準溶液和10 mL 50 ℃溫水，振蕩2～3 min；加入10 mL KOH 溶液（0.6 g/mL）和20 mL 乙醇， 60 ℃恒溫振蕩皂化30 min；冷卻至室溫后，加入15 mL 正己烷，振蕩5 min，加入10 mL 去離子水，振蕩、離心，取上清液；向殘留物中加入10 mL 正己烷重復提取，合并正己烷上清液，減壓濃縮至2～3 mL；加入5 mL 乙醇-水（1∶1， V/V）混合溶劑，振蕩、離心；將所有上清液轉移到3 g 硅膠柱中，用15 mL 正己烷-二氯甲烷（7∶3，V/V）混合溶劑洗脫，收集所有淋洗液，濃縮至1～2 mL。

其它食品基質的前處理參照文獻[22]的方法進行。將樣品均質后，稱取5 g 樣品至離心管中，加入10 μL MOSH/MOAH 標準溶液和20 mL KOH-水-乙醇溶液（25∶50∶90， m/V/V）， 60 ℃恒溫振蕩皂化30 min后，提取步驟與奶粉樣品一致。

1.3.2 環氧化

樣品中的生物烯烴會干擾MOAH 分析，必須通過環氧化將其去除[23]。移取1.3.1 節所得提取物1 mL 至離心管中，加入0.5 mL m-CPBA-乙醇溶液（200 g/L， m/V）， 40 ℃振蕩15 min 后，加入0.5 mL 乙醇和2 mL 100 g/L 硫代硫酸鈉溶液，繼續振蕩30 s，終止反應，離心，取上清液進行后續分析。

1.4 LC-GC分析

MOSH 和MOAH 的分離分析通過LC-GC-FID 聯用系統進行。液相色譜的色譜柱為Allure Si HPLC色譜柱（250 mm×2.1 mm， 5 μm， 60 ?，美國瑞思泰克公司），流動相為正己烷和二氯甲烷。將50～100 μL樣品注入液相色譜儀，在進樣1.5 min 內，以300 μL/min 的流速將流動相從最初的100%正己烷轉變為70%正己烷，保持至6.2 min；采用500 μL/min 的流速儀用100%二氯甲烷反沖9 min；切換至100%正己烷，以500 μL/min 的流速正向平衡色譜柱10 min，降低流速至300 μL/min，直至程序結束。MOSH 餾分的洗脫時間為2.0～3.5 min， MOAH 餾分的洗脫時間為4.5～6.0 min，通過紫外檢測器（230 nm）監測MOSH 和MOAH 的流出情況， MOSH 和MOAH 餾分通過HPLC-GC 接口從HPLC 全部轉移至GC，隨后，通過保留間隙技術和部分并行溶劑蒸發技術將大量的溶劑蒸汽排出，溶質保留在分析柱頭。保留間隙系統包括MXT Siltek 無涂層的預柱（10 m×0.53 mm，美國瑞斯泰克公司）、MXT-1 100%二甲基聚硅氧烷分析柱（15 m×0.25 mm×0.25 μm，美國瑞斯泰克公司）和溶劑排除出口（SVE）。GC 在MOSH 和MOAH 餾分轉移時壓力分別為84 和90 kPa（氫氣），大量溶劑蒸汽從SVE 出口排出， SVE 出口在餾分轉移前0.5 min 打開，在轉移結束0.3 min 后關閉。SVE 出口關閉后， MOSH 和MOAH 通道的分析壓力升至150 kPa 并保持恒定。柱溫箱升溫程序設定為：初始溫度60 ℃保持6 min，以15 ℃/min 升至120 ℃，再以25 ℃/min 升至370 ℃，保持6 min。FID 溫度為380 ℃。氮氣、氫氣（燃燒氣）和空氣（助燃氣）的流速分別設置為30、40 和400 mL/min。

1.5 數據處理

人工積分使用Lab solutions 軟件（Version 5.92，日本島津公司）進行數據處理。分別以Cycy、TBB 或2-MN 作為MOSH 或MOAH 的定量內標對MOSH 和MOAH 進行定量分析。通過積分計算鼓包峰與基線之間的空白面積與內標峰的面積比獲得礦物油的含量。

基于多重二階差商濾波的駝峰曲線形態的自動判別算法，將原始譜圖數據導入軟件后，自動積分得到MOSH 或MOAH 含量。

2 結果與討論

2.1 礦物油含量自動計算算法

在礦物油檢測中，由于儀器配套的數據分析軟件不具備自動建立礦物油譜線的功能，需要檢測人員手工建立礦物油譜線，這給礦物油檢測工作增加了難度和工作量。Bauwens等[19]提出通過在每個2D點重建1D色譜圖（類似于1D軸上的投影）來消除尖峰的方法。該方法可以提取出粗略的礦物油色譜圖，但對有些較寬的尖峰去除效果不佳。例如，在文獻[19]中，在1D 時間1300s附近，原始色譜數據（紅色曲線）包含一個較寬的尖峰，應被平滑去除，但該方法得到的礦物油色譜圖（綠色曲線）未能將其去除。

本研究設計了一種思路簡單、功能完備且能準確去除多種尖峰的算法。首先，通過對礦物油色譜數據進行谷值點搜索建立駝峰的輪廓；然后，對谷值點進行三次二階差商過濾，去除無效的谷值點（包括尖峰和豁口），得到色譜數據正確的下包絡線，即礦物油譜線；最后，通過對礦物油譜線和基線之間進行面積積分，根據用戶輸入的內標質量、樣品質量和礦物油譜線面積即可以計算出樣品中的礦物油的含量。本方法不僅能夠自動建立礦物油譜線，也能自動建立基線，還能根據需要對自動建立的譜線對參數進行人工微調。實驗結果表明，本方法能很好地提高礦物油含量分析的效率和精確度。

2.2 礦物油譜線的自動生成算法

礦物油譜線的生成主要分為2步：（1）尋找有效區域內的谷值點，建立礦物油駝峰的基本輪廓；（2）對基本輪廓中的無效點進行過濾，確定準確的礦物油譜線。

2.2.2 去除無效的谷值點

雖然谷值點基本反映了礦物油駝峰的形狀輪廓，但有一些數據點明顯不屬于礦物油在色譜圖中形成的駝峰。如圖3 中遠高于大部分綠色數據點的尖峰數據點，單個的峰以尖峰形式出現，簇狀的峰則以豁口的形式出現，這些都不屬于礦物油在色譜圖中的駝峰。觀察這些無效的谷值點可發現一個共同特征，即二階差商是絕對值較大的負數，說明這些無效谷值點在谷值點序列中既是局部極大值點，又是數據系列的拐點。因此，可以采用對谷值點序列進行二階差商過濾的方法去除這些無效的谷值點。序列的二階差商計算方法如下：

其中， D₀是根據經驗設置的判別閾值，為負數。當谷值點序列中一個點的二階差商小于D₀時，即可認為該點為無效的谷值點，將其從序列中刪除。通過調節閾值D₀可設置算法對豁口的屏蔽程度：閾值的絕對值越小，屏蔽程度越大；閾值的絕對值越大，屏蔽程度越小。

從谷值點序列中過濾二階差商絕對值較大的點，可以去除單個的尖峰。但是，豁口是由一簇尖峰形成的，一次差商過濾難以完全去掉整個豁口。實驗發現，對谷值點序列進行正向、反向、正向3次差商過濾，即可基本去除色譜數據中的所有豁口。對于一些特殊的豁口，如果采用默認參數去除效果不佳，還可通過適度調節閾值D0，增大或減小對豁口的去除程度，優化礦物油譜線。

圖3中紫色數據點及其形成的折線即為礦物油氣相色譜圖中的駝峰，也是色譜圖數據對應的正確的礦物油譜線。對于較寬的尖峰，本算法也能很好地去除豁口，獲得準確的駝峰譜圖。因此，本方法能夠有效地去除多種尖峰。

在實際處理數據時，若通過調整閾值得到的礦物油譜線未達到最佳，本程序還提供手動刪除無效數據點的功能，從而使礦物油譜線更符合實際。

2.3 基線的確定和積分計算

計算礦物油譜線對應的積分數值，還需確定積分使用的基線，通過計算礦物油譜線和基線之間包絡的面積換算出礦物油檢測的含量。

2.3.1 完全自動生成基線

大部分情況下，積分使用的基線就是包絡線序列的起始點和終止點之間連接的直線。但是，礦物油譜線駝峰對應的中間谷值點有時會低于包絡線序列的起始點和終止點之間連接的直線，這時基線將變為包絡線序列的起始點、終止點以及駝峰對應的中間谷值點3 點構成的折線，如圖4所示。

2.3.2 處理加載的基線數據

如圖5A 所示，實驗數據形成的基線存在干擾，不利于積分計算，為了確保檢測結果的可靠性，有必要對基線數據進行曲線擬合。通過觀察可知，二次曲線擬合基本可以滿足形狀上的要求（圖5B）。

在計算數值積分時，由于礦物油譜線和自動生成的基線數據均表現為折線形狀，需要在具體的積分位置計算折線上對應的y 值。

其中， j 是折線上對應的數據點的編號。對于加載的數據基線，可以記錄二次曲線擬合的3 個參數，直接采用二次曲線公式計算數值積分對應的y值：

2.3.3 計算積分面積

確定了礦物油譜線和基線即確定了數值積分需要的上、下邊界，就可以計算礦物油譜線與基線間包含的面積，依據積分面積可換算出礦物油檢測的含量。通常， MOSH 選擇Cycy 作為定量內標， MOAH 選擇2MN 或TBB作為定量內標。在加載檢測數據時，軟件可直接讀入色譜數據中的內標峰面積。在處理色譜數據時，手動輸入內標質量M_ISTD（ng）和樣品質量M_smp（g）。根據礦物油譜線面積A_smp 和標準物質的峰面積A_ISTD，利用式（8）計算出礦物油在樣品中的含量。

根據上述原理自動生成基線，獲得礦物油譜線和基線數據，對二者之間包含的區域進行自動數值積分，再根據所得到的積分面積換算出礦物油的含量。與文獻[19]的方法相比，本方法功能更完備，更利于提高用戶操作效率。

2.4 算法流程

（1）加載數據。首先按固定格式將色譜數據轉存為ASCII 碼，采用軟件加載ASCII 碼數據文件。在加載了色譜數據文件的基礎上，繼續加載基線數據和碳數譜圖的數據。（2）計算有效區間。利用式（2）和式（3）計算礦物油檢測數據的有效區間，提取對應的數據編號。（3）計算谷值點。利用式（4）計算礦物油色譜數據在有效區間內包含的谷值點。（4）計算礦物油譜線。設置合適的閾值，利用式（5）刪除谷值點序列中的無效數據點，必要時手動刪除谷值點序列中無效數據點。（5）確定基線數據。可以根據礦物油譜線的起始點、終止點和中間谷值點構建基線，也可由用戶直接加載基線數據。基線數據需要在有效區間內采用最小二乘法進行二次曲線擬合，同時保留3 個參數。（6）計算數值積分。缺省狀態下，積分區間為礦物油譜線的全部有效區間，可根據碳數譜圖數據選擇相應的積分區間。數值積分時，利用式（6）和式（7）計算y 值。（7）計算含量。根據用戶輸入的內標質量M_ISTD（ng）、樣品質量M_smp（g）、礦物油譜線面積A_smp 和標準物質的峰面積A_ISTD，利用式（8）計算樣品中的礦物油含量。（8）顯示并保存數據。計算結果直接在窗口中顯示，用戶可以選擇保存計算結果到數據文件。本方法已全部用軟件實現并投入使用中，軟件操作界面如圖6所示。

2.5 算法驗證與應用

測定了德國標物局（DRRR）國際能力驗證的食用油（RVEP233004）和ISO 20122 國際標準驗證的葵花籽油和橄欖油樣品，相關數據分別使用傳統的色譜工作站人工處理和本研究建立的計算機軟件自動計算，結果見表1。采用z′比分數衡量實驗室測定值與參考值的偏離程度，人工處理和軟件自動計算的z′比分數的絕對值均小于2，表明兩種數據處理方式均能滿足能力驗證要求。

以市售米粉樣品為實際樣品，采用潤滑油標準品（MOSH∶MOAH=89∶11）制備2 個濃度水平的加標樣品，用于加標回收實驗。按照1.3節方法處理后，注入HPLC-GC-FID 進行測定，采用所建立的計算機軟件自動計算，結果見表2，MOSH 和MOAH 的加標回收率分別為82.4%～108.7%和92.9%～95.7%，RSD≤3.5%。

選擇食用油、奶粉、奶酪和魚罐頭等真實樣品基質共42 個樣品，進行相應的前處理操作，將供試液注入HPLC-GC 聯用儀器，分別進行人工積分和軟件自動積分，再剔除低于定量限的數據，得到42 組MOSH 數據和28 組MOAH 數據。采用SPSS 軟件對人工積分和軟件積分的數據進行處理，結果見表3，MOSH 和MOAH 不同處理方式的差值均符合正態分布，采用配對t 檢驗方法進行顯著性差異分析，結果見表4， pgt;0.05，表明兩種數據處理方式不存在顯著性差異。

檢測結果表明，本算法適用于食用油、乳酪和奶粉等基質樣品中定量限以上的MOSH和MOAH檢測。

3結論

建立了一種礦物油檢測色譜數據分析的自動判別方法。分析結果表明，本方法不僅原理簡單，而且比現有方法更準確，功能更完備，在很大程度上提高了礦物油色譜數據分析的準確性和工作效率。（1）構建的駝峰生成方法符合礦物油數據的實際形態，可準確地描繪礦物油在色譜數據中的實際區域。通過構建礦物油色譜數據的“谷點”序列，在“谷點”序列中排除二階差商較小的“凸點”，有效消除了干擾數據，最終建立了符合實際的礦物油譜線。（2）開發的數據處理軟件工具可以在礦物油色譜數據中自動形成積分計算的有效區間和基線，還可以針對加載的基線數據，通過最小二乘法在有效區域局部形成基線的二次曲線。本軟件可由使用者根據碳數譜圖確定積分區間，方便確定碳數區間對應的物質含量。本方法是按照傳統數據分析算法設計的，還存在少量需要手動刪除部分干擾峰的情況。在積累大量的實驗數據后，通過機器學習等人工智能方法，開發更準確高效的數據分析模型，為礦物油的數據分析提供更科學高效的分析工具。