基于GLS算法的氣相色譜儀檢測偏差原因分析

王欣瑋，趙 欣，宋 悅

(北京市燃?xì)饧瘓F(tuán)有限責(zé)任公司， 北京 100035)

2020年，中國天然氣消費(fèi)量為3 281×108m3，較2019年增長7.8%[1]，國內(nèi)天然氣市場上并存著頁巖氣、煤層氣、液化天然氣和管道進(jìn)口天然氣等多種氣源[2]。這些氣源的熱值不盡相同，目前，在天然氣交接計(jì)量方面，國際上主要以能量計(jì)量為主，而我國通常采用體積計(jì)量，由于體積計(jì)量會將天然氣中無效的組分也考慮在內(nèi)，不能反映天然氣作為能源的最終價值，也不便與其他類型的能源進(jìn)行定量比較[3-4]。在能量計(jì)量中發(fā)熱量的獲取至關(guān)重要，然而直接測量法不易在現(xiàn)場實(shí)施，國內(nèi)天然氣計(jì)量站現(xiàn)場天然氣單位發(fā)熱量的獲取方式主要是用氣相色譜儀檢測天然氣的組成，然后通過GB/T 11062(ISO 6976)計(jì)算發(fā)熱量、壓縮因子、密度等物性參數(shù)[5]。天然氣組成測定結(jié)果的準(zhǔn)確性將直接影響到發(fā)熱量計(jì)算結(jié)果的準(zhǔn)確性。目前現(xiàn)場使用的在線氣相色譜儀普遍采用單點(diǎn)校正方法，采用經(jīng)過原點(diǎn)的一階函數(shù)，由峰面積計(jì)算組成含量，但在實(shí)際工況中多數(shù)色譜儀某些組分的響應(yīng)函數(shù)的截距并不為零，即在沒有樣品進(jìn)入系統(tǒng)的情況下，某些組分仍然有響應(yīng)峰面積。使用單點(diǎn)校正方法不再適合于多氣源能量計(jì)量[6]。根據(jù)GB/T 28766[7]中采用多點(diǎn)校正法由試驗(yàn)得到多個測試點(diǎn)的數(shù)據(jù)，以廣義最小二乘法(GLS)回歸計(jì)算得到色譜儀系統(tǒng)的真實(shí)分析函數(shù)和校準(zhǔn)函數(shù)，并以這些函數(shù)為基礎(chǔ)進(jìn)行各組分含量的定量計(jì)算。

筆者針對氣相色譜儀分析系統(tǒng)性能評價要求，結(jié)合GLS算法擬合天然氣各組分的響應(yīng)函數(shù)，對2臺氣相色譜儀測量的準(zhǔn)確性進(jìn)行評價。基于Visual Basic語言編寫GLS多項(xiàng)式響應(yīng)函數(shù)，在此基礎(chǔ)上建立天然氣組分的分析函數(shù)，研究2臺儀器回歸函數(shù)偏差情況與高位體積發(fā)熱量偏差情況之間的關(guān)系，以此評價2臺氣相色譜儀測量的準(zhǔn)確性。

1 廣義最小二乘法的基本原理

一般化回歸模型中最小二乘法(OLS)的估計(jì)量是無偏估計(jì)，但是其推導(dǎo)過程基于同方差假設(shè)進(jìn)行，針對異方差并不具有有效性[8-10]。

假設(shè)給定的線性回歸模型如下：

Y=Xβ+ε

(1)

E[ε]=0

(2)

Var[ε]=σ2V

(3)

(4)

其中，V是n×n階矩陣，如果V為對角線元素不等的對角矩陣，則y僅存在異方差，無自相關(guān)性，當(dāng)V的對角線外有非零元素時，y則有自相關(guān)性，其公式為：

(5)

(6)

σ2V是異方差，V是方差-協(xié)方差矩陣，必是正定的，其逆矩陣也是正定的，矩陣K滿足V=K′K=KK。

再定義

z=K-1yB=K-1Xg=K-1ε?z=Bβ+g

(7)

則有

E[g]=K-1E[ε]=0

(8)

Var[g]=Var[K-1ε]=K-1Var[ε]K-1=

σ2K-1VK-1=σ2K-1KKK-1=σ2I

(9)

此時滿足最小二乘法的假設(shè)條件，最小二乘法函數(shù)式為：

(10)

(X′V-1X)β=XV-1y

(11)

則廣義最小二乘法(GLS)的估計(jì)量為：

(12)

同時可得：

(X′V-1X)-1XV-1Xβ=β

(13)

σ2(X′V-1X)-1

(14)

廣義最小二乘法(GLS)的理論求解過程如上，說明最小二乘法(OLS)是廣義最小二乘法(GLS)的特殊情況之一。

2 廣義最小二乘法在多點(diǎn)校正曲線回歸中的應(yīng)用

響應(yīng)函數(shù)根據(jù)獨(dú)立變量與非獨(dú)立變量的選擇不同有2種表達(dá)方式：

選用組分含量為獨(dú)立變量的函數(shù)稱為校準(zhǔn)函數(shù)，標(biāo)準(zhǔn)物質(zhì)經(jīng)氣相色譜儀分析后的組分濃度與儀器響應(yīng)之間關(guān)系式為：

yi=F(xi)

(15)

選用儀器響應(yīng)為獨(dú)立變量的函數(shù)稱為分析函數(shù)，標(biāo)準(zhǔn)物質(zhì)經(jīng)氣相色譜儀分析后的儀器響應(yīng)與組分濃度之間關(guān)系式為：

xi=G(yi)

(16)

通過多個WMS經(jīng)氣相色譜儀分析后可以得到xi、yi的樣本矩陣，選用廣義最小二乘法(GLS)對校準(zhǔn)函數(shù)與分析函數(shù)進(jìn)行回歸分析。

2.1 響應(yīng)函數(shù)的線性回歸

觀察推導(dǎo)過程可知，V是假設(shè)值，在實(shí)際求解中無法直接獲取，此時針對響應(yīng)函數(shù)的分析過程中，需要對廣義最小二乘法進(jìn)行一定數(shù)值求解變換，使用可行廣義最小二乘法(GLS)，其主要求解思路如下：

(17)

2.2 響應(yīng)函數(shù)的非線性回歸

當(dāng)響應(yīng)函數(shù)回歸時采用非線性回歸方法，GLS算法應(yīng)用也需要對數(shù)值計(jì)算方法進(jìn)行相應(yīng)修改。主要求解思路為針對誤差平方和引入單位權(quán)重，將誤差平方和進(jìn)行一階泰勒展開后求偏導(dǎo)，此時的正則方程組是非滿秩的，無法求解。這一問題需要對正則方程組進(jìn)行變換，將其變換為y的多次微分方程進(jìn)行迭代求解，以獲得正則方程組中所需參數(shù)；此時可以將誤差平方和轉(zhuǎn)化為線性方程組，再使用GLS即可得到響應(yīng)函數(shù)的參數(shù)估計(jì)值。

3 多點(diǎn)校正分析函數(shù)和現(xiàn)行函數(shù)比較

按照GB/T 28766規(guī)定的方法對色譜儀進(jìn)行性能評價，以多個濃度的氣體標(biāo)準(zhǔn)物質(zhì)多次檢測得到的峰面積，與氣體標(biāo)準(zhǔn)物質(zhì)含量結(jié)合，用回歸算法計(jì)算得到儀器的多點(diǎn)校正真實(shí)分析函數(shù)，并與儀器當(dāng)前使用的單點(diǎn)校正函數(shù)進(jìn)行比較。色譜儀多點(diǎn)校正方法回歸得到的分析函數(shù)均為一次，函數(shù)形式為：

y=b0+b1x

(18)

其中：y為組分的摩爾分?jǐn)?shù)；x為該組分的峰面積，單位與色譜儀型號有關(guān)，因此這里只給出數(shù)據(jù)；b0和b1為回歸分析函數(shù)的系數(shù)。

色譜儀1和色譜儀2的函數(shù)如表1和表2所示。

表1 色譜儀1的多點(diǎn)校正和單點(diǎn)校正分析函數(shù)

表2 色譜儀2的多點(diǎn)校正和單點(diǎn)校正分析函數(shù)

為了直觀地看出不同函數(shù)系數(shù)的差別，繪制出函數(shù)的圖形以便比較2組分析函數(shù)之間的差別情況。2臺色譜儀的分析函數(shù)比較情況如圖1～圖6所示。

圖1 色譜儀1的不同N2分析函數(shù)比較Fig.1 Comparison of different N2 analysis functions of gas chromatography 1

圖2 色譜儀1的C2分析函數(shù)比較Fig.2 Comparison of different N2 analysis functions of gas chromatography 1

圖3 色譜儀1的不同C1分析函數(shù)比較Fig.3 Comparison of different C1 analysis functions of gas chromatography 1

由圖1～圖3中可以看出，色譜儀1的3個主要組分，在N2和C2摩爾分?jǐn)?shù)大于0.6%、C1摩爾分?jǐn)?shù)大于60%時，現(xiàn)行單點(diǎn)校正分析函數(shù)與多點(diǎn)校正法得到的真實(shí)分析函數(shù)之間的偏差隨著摩爾分?jǐn)?shù)的增加顯著增大。

圖4 色譜儀2的N2分析函數(shù)比較Fig.4 Comparison of different N2 analysis functions of gas chromatography 2

圖5 色譜儀2的C2分析函數(shù)比較Fig.5 Comparison of different C2 analysis functions of gas chromatography 2

圖6 色譜儀2的C1分析函數(shù)比較Fig.6 Comparison of different C1 analysis functions of gas chromatography 2

由圖4～圖6中可以看出，色譜儀2的3個主要組分，在N2摩爾分?jǐn)?shù)大于0.6%時，偏差隨著摩爾分?jǐn)?shù)的增加有所增大，C2的偏差情況基本保持在同樣的水平，C1摩爾分?jǐn)?shù)大于70%時，現(xiàn)行單點(diǎn)校正分析函數(shù)與多點(diǎn)校正法得到的真實(shí)分析函數(shù)之間的偏差隨著摩爾分?jǐn)?shù)的增加顯著減小。

通常天然氣中甲烷的體積分?jǐn)?shù)在80%～99%之間，所以甲烷的偏差情況很大程度上決定了發(fā)熱量等物性參數(shù)的偏差情況。

4 高位體積發(fā)熱量偏差情況比較

按照GB/T 28766—2018中規(guī)定的方法，用蒙特-卡洛法對2臺色譜儀在測量范圍內(nèi)的偏差進(jìn)行計(jì)算和比較，結(jié)果如表3所示。

表3 2臺色譜儀單位發(fā)熱量偏差 MJ/m3

由表3中可以看出，色譜儀1的1萬組假設(shè)數(shù)據(jù)的平均發(fā)熱量偏差絕對值比色譜儀2高1個數(shù)量級，色譜儀1的發(fā)熱量最大值和最小值與真實(shí)值的偏差分別為1.17%和1.15%，遠(yuǎn)高于色譜儀2的0.03%和0.08%，已經(jīng)超出GB/T 18603[12]對A級和B級計(jì)量站的發(fā)熱量最大允許誤差要求。

6 結(jié)論與建議

(1)相比單點(diǎn)校正函數(shù)，基于GLS算法的多點(diǎn)校正函數(shù)消除了自誤差項(xiàng)，精度較高。在未來能量計(jì)量實(shí)施的條件下，更適用于氣相色譜儀分析系統(tǒng)的性能評價工作。

(2)由2臺色譜儀的分析函數(shù)和物性參數(shù)計(jì)算結(jié)果的對比情況來看，現(xiàn)場安裝的在線氣相色譜儀目前所采用的方法存在檢測結(jié)果偏差超出標(biāo)準(zhǔn)要求的情況。其中當(dāng)天然氣中占主要成份的組分函數(shù)偏差較大時，相應(yīng)的發(fā)熱量偏差會超出標(biāo)準(zhǔn)規(guī)定的范圍。

(3)由于現(xiàn)場使用的在線組成分析儀器的軟件系統(tǒng)中一般尚不具備更改各組分校準(zhǔn)函數(shù)和響應(yīng)函數(shù)的功能，可能無法滿足能量計(jì)量過程中的數(shù)據(jù)準(zhǔn)確性。在天然氣能量計(jì)量方式正式實(shí)施的情況下，天然氣組成分析儀器應(yīng)定期開展評價工作，確保儀器工作在正常狀態(tài)下能夠給出準(zhǔn)確的分析數(shù)據(jù)。