基于偏最小二乘法的高光谱水面油膜厚度估算

邢学文， 刘松， 许德刚， 钱凯俊

(1.中国石油勘探开发研究院，北京 100083； 2.中国石油集团安全环保技术研究院，北京 102206)

0 引言

随着水上油气勘探、油气运输活动的日益活跃，溢油事故时有发生。国外比较有名的墨西哥湾溢油事故，国内影响较大的渤海湾溢油事故，都严重污染了附近海域的海水，对海洋生态环境造成了巨大的破坏。溢油事故发生后，溢油量估算是后期事故评估、处理的重要指标，其中水面油膜厚度是溢油量估算的关键参数。

目前，采用最多的方法是通过遥感等监测技术获得水面油膜颜色，根据波恩协议油膜颜色与厚度的对应关系表获得油膜厚度，结合监测获取的面积估算溢油量[1-2]。国内近年来也开展了大量油膜厚度的遥感光谱特征分析和反演建模研究，赵冬至等[3]应用安徽光机所生产的VF921-B地物光谱仪分别对辽河原油、轻柴油和润滑油进行了光谱测量和分析； 刘旭拢等[4]开展了水面浮油光谱测量及光谱特征分析； 臧影[5]开展了高光谱溢油图像波段选择在油膜厚度估算中的应用； 兰国新[6]开展了海上溢油遥感光谱信息挖掘与应用研究； 肖剑伟等[7]开展了基于生物光学模型的水面薄油膜厚度高光谱反演实验研究； 孙鹏等[8]应用AvaSpec光谱仪开展了高光谱油膜厚度估计模型分析，利用曲线拟合、神经网络和基于奇异值分解的迭代方法构造响应函数，建立了多个基于高光谱指标的厚度模型； 刘丙新等[9]开展了不同厚度海上油膜高光谱遥感波段敏感性研究。以上研究大多针对350～1 000 nm谱段，利用与油膜厚度相关性较大的光谱特征指标及特征指标组合选择不同建模方法进行油膜厚度-特征指标(组合)建模，而忽略了大量与油膜厚度相关性较小的谱段所携带的油膜厚度信息。特别是1 000～2 500 nm谱段，由于其在水面的反射率非常小，以往研究大多不予关注，但原油的烃类组分和官能团的吸收特征却基本都位于该谱段范围[10]。

本文以水面原油油膜高光谱测量实验为基础，获取不同厚度油膜的全谱段(350～2 500 nm)反射率光谱，针对实验数据自变量(光谱数据)多的特点，选择采用偏最小二乘方法(partial least squares，PLS)进行油膜厚度-遥感光谱反射率建模，为水面油膜厚度估算提供一种新的途径。

1 实验

1.1 实验装置及材料

油膜厚度高光谱遥感检测实验装置和材料包括： ASD FieldSpec3光谱仪、手柄支架、白板、石英卤素灯、1 000 ml烧杯、黑色消光布和油量注射器，实验油品为江苏油田原油，实验用水为自来水，图1为实验装置的示意图。

图1 油膜遥感检测实验示意图

其中ASD FieldSpec3光谱仪采集背景水和水面不同厚度油膜的反射率光谱，手柄支架固定光谱仪的光纤探头，保持探头垂直水面，2个石英卤素灯模拟太阳光源，在1 000 ml烧杯中装入一定量实验背景水，黑色消光布包裹烧杯底部和外围，消除外界光线干扰，注射器用于实验用油的定量化。由于油膜厚度直接测量难度较大，本次实验选择体积法进行油膜厚度估算，基本原理如图2所示。

图2 油膜厚度计算方法示意图

具体表达式为

(1)

式中：V代表注射到烧杯中的油量；r为烧杯的半径；h为油膜厚度。

1.2 实验数据

1.2.1 原油水面光谱特征

原油油膜模拟的厚度范围为0.05～0.6 mm，FieldSpec3光谱仪测量波谱范围为350～2 500 nm。不同厚度油膜的光谱曲线如图3所示，yy.001(黑色曲线)为背景水光谱，其他分别为不同厚度油膜的光谱曲线，随着油膜厚度增加，反射率不断降低，在可见光—近红外谱段范围比较直观，短波红外谱段由于水的强烈吸收，反射率非常小，变化不明显。

图3 不同厚度油膜反射率光谱曲线

1.2.2 光谱数据的线性相关分析

光谱数据是作为自变量进行油膜厚度回归建模，当自变量之间存在高度线性相关时，对回归系数统计检验造成困难，回归系数估计值的稳定性也会降低，不能很好地解释回归系数的物理含义，最终会对预测结果造成影响。

图4为油膜在350～2 500 nm光谱范围反射率变量的线性相关矩阵，其中350～950 nm谱段范围变量之间的相关性非常大，相关系数大多超过0.8，也就是至少600多个变量之间具有高度线性相关。

图4 线性相关性矩阵

产生高度线性相关的原因为： ①变量之间的物理含义决定它们之间的多重相关性； ②测量样本点个数较少[11]。

2 PLS模型建立

2.1 PLS

PLS是一种多元统计数据分析方法，集多元线性回归分析、典型相关分析和主成分分析的基本功能于一体。与传统多元线性回归模型相比，PLS回归的特点为： ①能够在自变量存在严重多重相关性的情况下进行回归建模； ②允许在样本点个数少于变量个数的情况下进行回归建模； ③PLS在最终模型中将包含原有的所有自变量； ④PLS回归模型中，每一个自变量的回归系数将更容易解释[12]。

具体建模的原理为： 设有p个自变量X={x1，x2，…，xp}、q个因变量Y={y1，y2，…，yq}和n个样本点，分别在X和Y中提取主成分分量t1和u1，要求t1和u1应尽可能大地携带各自数据表中的变异信息，同时t1和u1的相关程度能够达到最大，使得t1和u1应尽可能好地代表数据表X与Y，且自变量成分t1对因变量成分u1又具有最强的解释能力。在第一个主成分分量t1和u1被提取后，分别实施X和Y对t1的回归，如果回归方程已经达到满意的程度，则算法终止，否则，分别利用X和Y被t1解释后的残余信息进行第二轮的成分提取，如此往复，直到达到一个较为满意的精度为止。若对X提取了m个主成分分量t1，t2，… ，tm，PLS将通过实施yk对t1，t2，…，tm的回归，然后再表达成yk对原变量x1，x2，…，xm的回归方程(k=1，2，…，q)[13-14]。

2.2 油膜厚度-反射率回归建模

2.2.1 样本点分布结构观察与特异点的发现

通过将高维数据系统降维至二维平面上，就可以对样本点的分布结构进行直接考察。高维数据系统提取成分t1和t2后，绘制以t1为横坐标、t2为纵坐标的t1-t2平面图，如图5所示。

图5 t1-t2平面图和T2椭圆

图5中(t1(i)，t2(i))代表了每一个样本点的位置，当2个样本点位置很接近时，它们在原自变量空间的高维性质就可能很近似，因此在t1-t2平面图上可以观察样本点的分布情况和相似性结构。此外，在t1-t2平面图上绘制T2椭圆还可以判断样本是否为特异点[14-15]。如果所有的样本点均落在椭圆区内，则认为所有样本点的分布是均匀的，否则，认为落在椭圆区外的样本点为奇异值，它们的取值远离所有样本点的平均水平，建模时应将其剔除。图中样本点的取值分布基本上是均匀的，绝大多数样本点都落在椭圆区内，唯有一个样本点落在了T2椭圆区外，对应样点编号为yy.002，是注射0.05 ml油量时测量的数据组，有可能是油量太少，油膜在烧杯中扩散不完全所致，这一样本在后期建模中需要剔除。

2.2.2 最佳主成分数确定及其解释能力

PLS回归建模中，并不需要将所有主成分分量进行回归建模，究竟应该选择多少个主成分分量，可以通过考察增加一个新的主成分分量后，对模型的预测功能是否有明显的改进来确定[12]。

图6 模型最佳主成分个数的确定

通过评价主成分分量t1，t2，...，tm，对X和Y的解释能力可以判断PLS模型的精度[14]。表1给出了前5个主成分分量对X和Y的解释能力，对自变量X来说，第一主成分分量的解释能力最强，反映了32.6%的自变量信息和49.0%的因变量信息，从累积解释能力来看，前5个主成分分量累积解释了74.0%的自变量信息和99.8%的因变量信息，达到了较高的解释水平，说明利用PSL拟合的回归模型能够概括原始数据所携带的大部分信息，所建模型的精度较高。

表1 PLS模型各主成分分量对X和Y的解释能力

2.2.3X与Y之间相关关系判断

自变量集合X与因变量集合Y之间是否存在相关关系，是检验是否可以建立Y对X的回归方程的基本条件[14]。主成分分量t1和u1分别携带了X和Y主要的典型成分特征，通过绘制t1-u1平面图，就可以从整体出发判断X与Y的相关关系，如图7所示。如果X和Y之间存在潜在的线性关系，模型的主成分分量数等于模型的维数，样本点在t1-u1图中的排列会近似一条直线，如果X和Y之间存在潜在的非线性关系，那就需要额外的主成分分量来描述非线性特征[16]。

(a) t1-u1相关关系 (b) t2-u2相关关系 (c) t3-u3相关关系

(d) t4-u4相关关系 (e) t5-u5相关关系

从表1中可以看出，PLS模型确定了前5个主成分分量，其中t1主成分分量携带了最大的变异信息，t1-u1平面图中，样本点呈现出比较明显的非线性特征，说明因变量(油膜厚度)与自变量(反射率)之间主要是存在非线性关系，还需要其他4个主成分分量来描述非线性特征(图7)，直到回归模型达到满意的程度。

2.2.4 PLS模型建立

实测油膜厚度与模型预测油膜厚度对比如图8所示。

(a) 建模数据 (b) 验证数据

根据最佳主成分分量进行模型建立，结果表明所建立的模型具有较高的精度，累积预测能力(Q2)达到92.8%(表1)。20个建模样本实际厚度与模型预测厚度最大相对误差为9.49%，最小相对误差为0.16%。20个样点的均方根误差(root mean squared error，RMSE)为0.01。实测值与预测值具有很好的对应关系，基本上均匀分布在1∶ 1线的两侧(图8)，这也说明模型具有相对较高的精度。

模型精度高并不意味着模型具有较好的预测能力，为了验证模型的预测能力及评价模型的稳定性，将未参与建模的6组样本数据代入上述建立的PLS模型进行验证。结果显示，实测油膜厚度与预测油膜厚度的相关系数(R)达到0.91，RMSE为0.04，总体上来说，所建立的PLS模型可以较好地反演油膜厚度，具有较好的预测能力和稳定性，可以用于水面油膜厚度估算。

2.2.5 自变量因子重要性程度分析

自变量因子对因变量的解释能力越强，其重要性越大。从PLS建模过程可知，若所提取的主成分分量th对Y的解释能力很强，而自变量因子xj在构造th时又起到了相当重要的作用，则xj对Y的解释能力就很大。xj的重要性程度可以用变量投影重要性指标(variable importance in projection，VIP)来测定，如果每个xj在解释Y时的作用都相同，则所有VIPj均等于1； 否则，对于VIPj很大(>1)的xj，在解释Y时候就有更加重要的作用(图9)[14-15]。

图9 PLS模型自变量因子的重要性指标

从350～2 500 nm谱段范围反射率的重要性指标可以看出，反映油膜厚度最显著的谱段为近红外谱段，这些谱段在油膜厚度建模过程中起到了重要的作用。

此外，自变量X与因变量Y的关系密切程度还可以通过w*c图来展示，其中w*为构造不同主成分分量t时各自变量的权重值，这些权值使构造的t和u相关性最大，间接指示了X和Y的相关性大小，c为构造主成分分量u时各变量的权重值[14]，针对主成分分量t1和t2，制作w*c1-w*c2平面图(图10)，可以看出油膜厚度与近红外谱段距离最近，也就是相关性最好，这可能与原油烃类组分和官能团的吸收特征位于近红外谱段有关[10]。

图10 w*c1-w*c2平面图

3 模型对比与评价

除了PLS模型，本次研究还选择常用的曲线拟合方法对油膜厚度和反射率进行建模。根据各个波段反射率和油膜厚度的相关系数，发现1 086 nm处的相关系数(R)最大(0.79)，以1 086 nm处的反射率为自变量，油膜厚度为因变量，分别进行指数、对数、乘幂和多项式拟合，其中指数模型精度最高，决定系数(R2)达到0.918 9(图11)。

图11 油膜厚度-反射率曲线拟合

PLS模型和指数模型的建模精度都比较高，R2分别为0.998 0和0.918 9，但PLS模型的RMSE要明显小于指数模型。6个验证样本分别带入建立的2个模型，PLS模型的RMSE同样明显小于指数模型(表2)，对比后认为PLS模型具有更高的建模精度，而且模型的稳定性也相对更好。

表2 不同建模方法比较

4 结论

1)油膜遥感检测实验发现，0.05～0.60 mm厚度范围的新鲜原油油膜，随着油膜厚度的增加，350～1 000 nm谱段范围的反射率不断变小，光谱上就可以直观反映油膜厚度变化； 而1 000～2 500 nm谱段范围，由于反射率过小，光谱特征不明显。

2)偏最小二乘法(PLS)适合于油膜厚度-全谱段反射率光谱数据组的回归建模，通过最佳主成分分量确定和提取，最大程度利用了所有油膜光谱中携带的油膜厚度信息进行建模。对于建模结果，还可以直观地分析各个自变量在建模中的重要性程度，来判断模型的合理性。研究发现，传统油膜厚度建模时经常被剔除的1 000～2 500 nm谱段，虽然其反射率很小，但对油膜厚度的PLS模型贡献比较大，而且由于烃类组分和官能团的吸收特征全部位于这一谱段范围，模型解释更为合理，因而选择其他方法进行建模时，这一谱段也应该给予关注。

3)相对于传统的曲线拟合建模方法，PLS模型无论是建模样本还是验证样本的误差均优于传统的经验模型，适合于高光谱数据的水面油膜厚度估算。