王欣玮，赵 欣，宋 悦

(北京市燃气集团有限责任公司， 北京 100035)

2020年，中国天然气消费量为3 281×108m3，较2019年增长7.8%[1]，国内天然气市场上并存着页岩气、煤层气、液化天然气和管道进口天然气等多种气源[2]。这些气源的热值不尽相同，目前，在天然气交接计量方面，国际上主要以能量计量为主，而我国通常采用体积计量，由于体积计量会将天然气中无效的组分也考虑在内，不能反映天然气作为能源的最终价值，也不便与其他类型的能源进行定量比较[3-4]。在能量计量中发热量的获取至关重要，然而直接测量法不易在现场实施，国内天然气计量站现场天然气单位发热量的获取方式主要是用气相色谱仪检测天然气的组成，然后通过GB/T 11062(ISO 6976)计算发热量、压缩因子、密度等物性参数[5]。天然气组成测定结果的准确性将直接影响到发热量计算结果的准确性。目前现场使用的在线气相色谱仪普遍采用单点校正方法，采用经过原点的一阶函数，由峰面积计算组成含量，但在实际工况中多数色谱仪某些组分的响应函数的截距并不为零，即在没有样品进入系统的情况下，某些组分仍然有响应峰面积。使用单点校正方法不再适合于多气源能量计量[6]。根据GB/T 28766[7]中采用多点校正法由试验得到多个测试点的数据，以广义最小二乘法(GLS)回归计算得到色谱仪系统的真实分析函数和校准函数，并以这些函数为基础进行各组分含量的定量计算。

笔者针对气相色谱仪分析系统性能评价要求，结合GLS算法拟合天然气各组分的响应函数，对2台气相色谱仪测量的准确性进行评价。基于Visual Basic语言编写GLS多项式响应函数，在此基础上建立天然气组分的分析函数，研究2台仪器回归函数偏差情况与高位体积发热量偏差情况之间的关系，以此评价2台气相色谱仪测量的准确性。

1 广义最小二乘法的基本原理

一般化回归模型中最小二乘法(OLS)的估计量是无偏估计，但是其推导过程基于同方差假设进行，针对异方差并不具有有效性[8-10]。

假设给定的线性回归模型如下：

Y=Xβ+ε

(1)

E[ε]=0

(2)

Var[ε]=σ2V

(3)

(4)

其中，V是n×n阶矩阵，如果V为对角线元素不等的对角矩阵，则y仅存在异方差，无自相关性，当V的对角线外有非零元素时，y则有自相关性，其公式为：

(5)

(6)

σ2V是异方差，V是方差-协方差矩阵，必是正定的，其逆矩阵也是正定的，矩阵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)

此时满足最小二乘法的假设条件，最小二乘法函数式为：

(10)

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

(11)

则广义最小二乘法(GLS)的估计量为：

(12)

同时可得：

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

(13)

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

(14)

广义最小二乘法(GLS)的理论求解过程如上，说明最小二乘法(OLS)是广义最小二乘法(GLS)的特殊情况之一。

2 广义最小二乘法在多点校正曲线回归中的应用

响应函数根据独立变量与非独立变量的选择不同有2种表达方式：

选用组分含量为独立变量的函数称为校准函数，标准物质经气相色谱仪分析后的组分浓度与仪器响应之间关系式为：

yi=F(xi)

(15)

选用仪器响应为独立变量的函数称为分析函数，标准物质经气相色谱仪分析后的仪器响应与组分浓度之间关系式为：

xi=G(yi)

(16)

通过多个WMS经气相色谱仪分析后可以得到xi、yi的样本矩阵，选用广义最小二乘法(GLS)对校准函数与分析函数进行回归分析。

2.1 响应函数的线性回归

观察推导过程可知，V是假设值，在实际求解中无法直接获取，此时针对响应函数的分析过程中，需要对广义最小二乘法进行一定数值求解变换，使用可行广义最小二乘法(GLS)，其主要求解思路如下：

(17)

2.2 响应函数的非线性回归

当响应函数回归时采用非线性回归方法，GLS算法应用也需要对数值计算方法进行相应修改。主要求解思路为针对误差平方和引入单位权重，将误差平方和进行一阶泰勒展开后求偏导，此时的正则方程组是非满秩的，无法求解。这一问题需要对正则方程组进行变换，将其变换为y的多次微分方程进行迭代求解，以获得正则方程组中所需参数；此时可以将误差平方和转化为线性方程组，再使用GLS即可得到响应函数的参数估计值。

3 多点校正分析函数和现行函数比较

按照GB/T 28766规定的方法对色谱仪进行性能评价，以多个浓度的气体标准物质多次检测得到的峰面积，与气体标准物质含量结合，用回归算法计算得到仪器的多点校正真实分析函数，并与仪器当前使用的单点校正函数进行比较。色谱仪多点校正方法回归得到的分析函数均为一次，函数形式为：

y=b0+b1x

(18)

其中：y为组分的摩尔分数；x为该组分的峰面积，单位与色谱仪型号有关，因此这里只给出数据；b0和b1为回归分析函数的系数。

色谱仪1和色谱仪2的函数如表1和表2所示。

表1 色谱仪1的多点校正和单点校正分析函数

表2 色谱仪2的多点校正和单点校正分析函数

为了直观地看出不同函数系数的差别，绘制出函数的图形以便比较2组分析函数之间的差别情况。2台色谱仪的分析函数比较情况如图1～图6所示。

图1 色谱仪1的不同N2分析函数比较Fig.1 Comparison of different N2 analysis functions of gas chromatography 1

图2 色谱仪1的C2分析函数比较Fig.2 Comparison of different N2 analysis functions of gas chromatography 1

图3 色谱仪1的不同C1分析函数比较Fig.3 Comparison of different C1 analysis functions of gas chromatography 1

由图1～图3中可以看出，色谱仪1的3个主要组分，在N2和C2摩尔分数大于0.6%、C1摩尔分数大于60%时，现行单点校正分析函数与多点校正法得到的真实分析函数之间的偏差随着摩尔分数的增加显著增大。

图4 色谱仪2的N2分析函数比较Fig.4 Comparison of different N2 analysis functions of gas chromatography 2

图5 色谱仪2的C2分析函数比较Fig.5 Comparison of different C2 analysis functions of gas chromatography 2

图6 色谱仪2的C1分析函数比较Fig.6 Comparison of different C1 analysis functions of gas chromatography 2

由图4～图6中可以看出，色谱仪2的3个主要组分，在N2摩尔分数大于0.6%时，偏差随着摩尔分数的增加有所增大，C2的偏差情况基本保持在同样的水平，C1摩尔分数大于70%时，现行单点校正分析函数与多点校正法得到的真实分析函数之间的偏差随着摩尔分数的增加显著减小。

通常天然气中甲烷的体积分数在80%～99%之间，所以甲烷的偏差情况很大程度上决定了发热量等物性参数的偏差情况。

4 高位体积发热量偏差情况比较

按照GB/T 28766—2018中规定的方法，用蒙特-卡洛法对2台色谱仪在测量范围内的偏差进行计算和比较，结果如表3所示。

表3 2台色谱仪单位发热量偏差 MJ/m3

由表3中可以看出，色谱仪1的1万组假设数据的平均发热量偏差绝对值比色谱仪2高1个数量级，色谱仪1的发热量最大值和最小值与真实值的偏差分别为1.17%和1.15%，远高于色谱仪2的0.03%和0.08%，已经超出GB/T 18603[12]对A级和B级计量站的发热量最大允许误差要求。

6 结论与建议

(1)相比单点校正函数，基于GLS算法的多点校正函数消除了自误差项，精度较高。在未来能量计量实施的条件下，更适用于气相色谱仪分析系统的性能评价工作。

(2)由2台色谱仪的分析函数和物性参数计算结果的对比情况来看，现场安装的在线气相色谱仪目前所采用的方法存在检测结果偏差超出标准要求的情况。其中当天然气中占主要成份的组分函数偏差较大时，相应的发热量偏差会超出标准规定的范围。

(3)由于现场使用的在线组成分析仪器的软件系统中一般尚不具备更改各组分校准函数和响应函数的功能，可能无法满足能量计量过程中的数据准确性。在天然气能量计量方式正式实施的情况下，天然气组成分析仪器应定期开展评价工作，确保仪器工作在正常状态下能够给出准确的分析数据。