邢學(xué)文， 劉松， 許德剛， 錢(qián)凱俊

(1.中國(guó)石油勘探開(kāi)發(fā)研究院，北京 100083； 2.中國(guó)石油集團(tuán)安全環(huán)保技術(shù)研究院，北京 102206)

0 引言

隨著水上油氣勘探、油氣運(yùn)輸活動(dòng)的日益活躍，溢油事故時(shí)有發(fā)生。國(guó)外比較有名的墨西哥灣溢油事故，國(guó)內(nèi)影響較大的渤海灣溢油事故，都嚴(yán)重污染了附近海域的海水，對(duì)海洋生態(tài)環(huán)境造成了巨大的破壞。溢油事故發(fā)生后，溢油量估算是后期事故評(píng)估、處理的重要指標(biāo)，其中水面油膜厚度是溢油量估算的關(guān)鍵參數(shù)。

目前，采用最多的方法是通過(guò)遙感等監(jiān)測(cè)技術(shù)獲得水面油膜顏色，根據(jù)波恩協(xié)議油膜顏色與厚度的對(duì)應(yīng)關(guān)系表獲得油膜厚度，結(jié)合監(jiān)測(cè)獲取的面積估算溢油量[1-2]。國(guó)內(nèi)近年來(lái)也開(kāi)展了大量油膜厚度的遙感光譜特征分析和反演建模研究，趙冬至等[3]應(yīng)用安徽光機(jī)所生產(chǎn)的VF921-B地物光譜儀分別對(duì)遼河原油、輕柴油和潤(rùn)滑油進(jìn)行了光譜測(cè)量和分析； 劉旭攏等[4]開(kāi)展了水面浮油光譜測(cè)量及光譜特征分析； 臧影[5]開(kāi)展了高光譜溢油圖像波段選擇在油膜厚度估算中的應(yīng)用； 蘭國(guó)新[6]開(kāi)展了海上溢油遙感光譜信息挖掘與應(yīng)用研究； 肖劍偉等[7]開(kāi)展了基于生物光學(xué)模型的水面薄油膜厚度高光譜反演實(shí)驗(yàn)研究； 孫鵬等[8]應(yīng)用AvaSpec光譜儀開(kāi)展了高光譜油膜厚度估計(jì)模型分析，利用曲線擬合、神經(jīng)網(wǎng)絡(luò)和基于奇異值分解的迭代方法構(gòu)造響應(yīng)函數(shù)，建立了多個(gè)基于高光譜指標(biāo)的厚度模型； 劉丙新等[9]開(kāi)展了不同厚度海上油膜高光譜遙感波段敏感性研究。以上研究大多針對(duì)350～1 000 nm譜段，利用與油膜厚度相關(guān)性較大的光譜特征指標(biāo)及特征指標(biāo)組合選擇不同建模方法進(jìn)行油膜厚度-特征指標(biāo)(組合)建模，而忽略了大量與油膜厚度相關(guān)性較小的譜段所攜帶的油膜厚度信息。特別是1 000～2 500 nm譜段，由于其在水面的反射率非常小，以往研究大多不予關(guān)注，但原油的烴類(lèi)組分和官能團(tuán)的吸收特征卻基本都位于該譜段范圍[10]。

本文以水面原油油膜高光譜測(cè)量實(shí)驗(yàn)為基礎(chǔ)，獲取不同厚度油膜的全譜段(350～2 500 nm)反射率光譜，針對(duì)實(shí)驗(yàn)數(shù)據(jù)自變量(光譜數(shù)據(jù))多的特點(diǎn)，選擇采用偏最小二乘方法(partial least squares，PLS)進(jìn)行油膜厚度-遙感光譜反射率建模，為水面油膜厚度估算提供一種新的途徑。

1 實(shí)驗(yàn)

1.1 實(shí)驗(yàn)裝置及材料

油膜厚度高光譜遙感檢測(cè)實(shí)驗(yàn)裝置和材料包括： ASD FieldSpec3光譜儀、手柄支架、白板、石英鹵素?zé)簟? 000 ml燒杯、黑色消光布和油量注射器，實(shí)驗(yàn)油品為江蘇油田原油，實(shí)驗(yàn)用水為自來(lái)水，圖1為實(shí)驗(yàn)裝置的示意圖。

圖1 油膜遙感檢測(cè)實(shí)驗(yàn)示意圖

其中ASD FieldSpec3光譜儀采集背景水和水面不同厚度油膜的反射率光譜，手柄支架固定光譜儀的光纖探頭，保持探頭垂直水面，2個(gè)石英鹵素?zé)裟M太陽(yáng)光源，在1 000 ml燒杯中裝入一定量實(shí)驗(yàn)背景水，黑色消光布包裹燒杯底部和外圍，消除外界光線干擾，注射器用于實(shí)驗(yàn)用油的定量化。由于油膜厚度直接測(cè)量難度較大，本次實(shí)驗(yàn)選擇體積法進(jìn)行油膜厚度估算，基本原理如圖2所示。

圖2 油膜厚度計(jì)算方法示意圖

具體表達(dá)式為

(1)

式中：V代表注射到燒杯中的油量；r為燒杯的半徑；h為油膜厚度。

1.2 實(shí)驗(yàn)數(shù)據(jù)

1.2.1 原油水面光譜特征

原油油膜模擬的厚度范圍為0.05～0.6 mm，F(xiàn)ieldSpec3光譜儀測(cè)量波譜范圍為350～2 500 nm。不同厚度油膜的光譜曲線如圖3所示，yy.001(黑色曲線)為背景水光譜，其他分別為不同厚度油膜的光譜曲線，隨著油膜厚度增加，反射率不斷降低，在可見(jiàn)光—近紅外譜段范圍比較直觀，短波紅外譜段由于水的強(qiáng)烈吸收，反射率非常小，變化不明顯。

圖3 不同厚度油膜反射率光譜曲線

1.2.2 光譜數(shù)據(jù)的線性相關(guān)分析

光譜數(shù)據(jù)是作為自變量進(jìn)行油膜厚度回歸建模，當(dāng)自變量之間存在高度線性相關(guān)時(shí)，對(duì)回歸系數(shù)統(tǒng)計(jì)檢驗(yàn)造成困難，回歸系數(shù)估計(jì)值的穩(wěn)定性也會(huì)降低，不能很好地解釋回歸系數(shù)的物理含義，最終會(huì)對(duì)預(yù)測(cè)結(jié)果造成影響。

圖4為油膜在350～2 500 nm光譜范圍反射率變量的線性相關(guān)矩陣，其中350～950 nm譜段范圍變量之間的相關(guān)性非常大，相關(guān)系數(shù)大多超過(guò)0.8，也就是至少600多個(gè)變量之間具有高度線性相關(guān)。

圖4 線性相關(guān)性矩陣

產(chǎn)生高度線性相關(guān)的原因?yàn)椋?①變量之間的物理含義決定它們之間的多重相關(guān)性； ②測(cè)量樣本點(diǎn)個(gè)數(shù)較少[11]。

2 PLS模型建立

2.1 PLS

PLS是一種多元統(tǒng)計(jì)數(shù)據(jù)分析方法，集多元線性回歸分析、典型相關(guān)分析和主成分分析的基本功能于一體。與傳統(tǒng)多元線性回歸模型相比，PLS回歸的特點(diǎn)為： ①能夠在自變量存在嚴(yán)重多重相關(guān)性的情況下進(jìn)行回歸建模； ②允許在樣本點(diǎn)個(gè)數(shù)少于變量個(gè)數(shù)的情況下進(jìn)行回歸建模； ③PLS在最終模型中將包含原有的所有自變量； ④PLS回歸模型中，每一個(gè)自變量的回歸系數(shù)將更容易解釋[12]。

具體建模的原理為： 設(shè)有p個(gè)自變量X={x1，x2，…，xp}、q個(gè)因變量Y={y1，y2，…，yq}和n個(gè)樣本點(diǎn)，分別在X和Y中提取主成分分量t1和u1，要求t1和u1應(yīng)盡可能大地?cái)y帶各自數(shù)據(jù)表中的變異信息，同時(shí)t1和u1的相關(guān)程度能夠達(dá)到最大，使得t1和u1應(yīng)盡可能好地代表數(shù)據(jù)表X與Y，且自變量成分t1對(duì)因變量成分u1又具有最強(qiáng)的解釋能力。在第一個(gè)主成分分量t1和u1被提取后，分別實(shí)施X和Y對(duì)t1的回歸，如果回歸方程已經(jīng)達(dá)到滿意的程度，則算法終止，否則，分別利用X和Y被t1解釋后的殘余信息進(jìn)行第二輪的成分提取，如此往復(fù)，直到達(dá)到一個(gè)較為滿意的精度為止。若對(duì)X提取了m個(gè)主成分分量t1，t2，… ，tm，PLS將通過(guò)實(shí)施yk對(duì)t1，t2，…，tm的回歸，然后再表達(dá)成yk對(duì)原變量x1，x2，…，xm的回歸方程(k=1，2，…，q)[13-14]。

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

2.2.1 樣本點(diǎn)分布結(jié)構(gòu)觀察與特異點(diǎn)的發(fā)現(xiàn)

通過(guò)將高維數(shù)據(jù)系統(tǒng)降維至二維平面上，就可以對(duì)樣本點(diǎn)的分布結(jié)構(gòu)進(jìn)行直接考察。高維數(shù)據(jù)系統(tǒng)提取成分t1和t2后，繪制以t1為橫坐標(biāo)、t2為縱坐標(biāo)的t1-t2平面圖，如圖5所示。

圖5 t1-t2平面圖和T2橢圓

圖5中(t1(i)，t2(i))代表了每一個(gè)樣本點(diǎn)的位置，當(dāng)2個(gè)樣本點(diǎn)位置很接近時(shí)，它們?cè)谠宰兞靠臻g的高維性質(zhì)就可能很近似，因此在t1-t2平面圖上可以觀察樣本點(diǎn)的分布情況和相似性結(jié)構(gòu)。此外，在t1-t2平面圖上繪制T2橢圓還可以判斷樣本是否為特異點(diǎn)[14-15]。如果所有的樣本點(diǎn)均落在橢圓區(qū)內(nèi)，則認(rèn)為所有樣本點(diǎn)的分布是均勻的，否則，認(rèn)為落在橢圓區(qū)外的樣本點(diǎn)為奇異值，它們的取值遠(yuǎn)離所有樣本點(diǎn)的平均水平，建模時(shí)應(yīng)將其剔除。圖中樣本點(diǎn)的取值分布基本上是均勻的，絕大多數(shù)樣本點(diǎn)都落在橢圓區(qū)內(nèi)，唯有一個(gè)樣本點(diǎn)落在了T2橢圓區(qū)外，對(duì)應(yīng)樣點(diǎn)編號(hào)為yy.002，是注射0.05 ml油量時(shí)測(cè)量的數(shù)據(jù)組，有可能是油量太少，油膜在燒杯中擴(kuò)散不完全所致，這一樣本在后期建模中需要剔除。

2.2.2 最佳主成分?jǐn)?shù)確定及其解釋能力

PLS回歸建模中，并不需要將所有主成分分量進(jìn)行回歸建模，究竟應(yīng)該選擇多少個(gè)主成分分量，可以通過(guò)考察增加一個(gè)新的主成分分量后，對(duì)模型的預(yù)測(cè)功能是否有明顯的改進(jìn)來(lái)確定[12]。

圖6 模型最佳主成分個(gè)數(shù)的確定

通過(guò)評(píng)價(jià)主成分分量t1，t2，...，tm，對(duì)X和Y的解釋能力可以判斷PLS模型的精度[14]。表1給出了前5個(gè)主成分分量對(duì)X和Y的解釋能力，對(duì)自變量X來(lái)說(shuō)，第一主成分分量的解釋能力最強(qiáng)，反映了32.6%的自變量信息和49.0%的因變量信息，從累積解釋能力來(lái)看，前5個(gè)主成分分量累積解釋了74.0%的自變量信息和99.8%的因變量信息，達(dá)到了較高的解釋水平，說(shuō)明利用PSL擬合的回歸模型能夠概括原始數(shù)據(jù)所攜帶的大部分信息，所建模型的精度較高。

表1 PLS模型各主成分分量對(duì)X和Y的解釋能力

2.2.3X與Y之間相關(guān)關(guān)系判斷

自變量集合X與因變量集合Y之間是否存在相關(guān)關(guān)系，是檢驗(yàn)是否可以建立Y對(duì)X的回歸方程的基本條件[14]。主成分分量t1和u1分別攜帶了X和Y主要的典型成分特征，通過(guò)繪制t1-u1平面圖，就可以從整體出發(fā)判斷X與Y的相關(guān)關(guān)系，如圖7所示。如果X和Y之間存在潛在的線性關(guān)系，模型的主成分分量數(shù)等于模型的維數(shù)，樣本點(diǎn)在t1-u1圖中的排列會(huì)近似一條直線，如果X和Y之間存在潛在的非線性關(guān)系，那就需要額外的主成分分量來(lái)描述非線性特征[16]。

(a) t1-u1相關(guān)關(guān)系 (b) t2-u2相關(guān)關(guān)系 (c) t3-u3相關(guān)關(guān)系

(d) t4-u4相關(guān)關(guān)系 (e) t5-u5相關(guān)關(guān)系

從表1中可以看出，PLS模型確定了前5個(gè)主成分分量，其中t1主成分分量攜帶了最大的變異信息，t1-u1平面圖中，樣本點(diǎn)呈現(xiàn)出比較明顯的非線性特征，說(shuō)明因變量(油膜厚度)與自變量(反射率)之間主要是存在非線性關(guān)系，還需要其他4個(gè)主成分分量來(lái)描述非線性特征(圖7)，直到回歸模型達(dá)到滿意的程度。

2.2.4 PLS模型建立

實(shí)測(cè)油膜厚度與模型預(yù)測(cè)油膜厚度對(duì)比如圖8所示。

(a) 建模數(shù)據(jù) (b) 驗(yàn)證數(shù)據(jù)

根據(jù)最佳主成分分量進(jìn)行模型建立，結(jié)果表明所建立的模型具有較高的精度，累積預(yù)測(cè)能力(Q2)達(dá)到92.8%(表1)。20個(gè)建模樣本實(shí)際厚度與模型預(yù)測(cè)厚度最大相對(duì)誤差為9.49%，最小相對(duì)誤差為0.16%。20個(gè)樣點(diǎn)的均方根誤差(root mean squared error，RMSE)為0.01。實(shí)測(cè)值與預(yù)測(cè)值具有很好的對(duì)應(yīng)關(guān)系，基本上均勻分布在1∶ 1線的兩側(cè)(圖8)，這也說(shuō)明模型具有相對(duì)較高的精度。

模型精度高并不意味著模型具有較好的預(yù)測(cè)能力，為了驗(yàn)證模型的預(yù)測(cè)能力及評(píng)價(jià)模型的穩(wěn)定性，將未參與建模的6組樣本數(shù)據(jù)代入上述建立的PLS模型進(jìn)行驗(yàn)證。結(jié)果顯示，實(shí)測(cè)油膜厚度與預(yù)測(cè)油膜厚度的相關(guān)系數(shù)(R)達(dá)到0.91，RMSE為0.04，總體上來(lái)說(shuō)，所建立的PLS模型可以較好地反演油膜厚度，具有較好的預(yù)測(cè)能力和穩(wěn)定性，可以用于水面油膜厚度估算。

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

自變量因子對(duì)因變量的解釋能力越強(qiáng)，其重要性越大。從PLS建模過(guò)程可知，若所提取的主成分分量th對(duì)Y的解釋能力很強(qiáng)，而自變量因子xj在構(gòu)造th時(shí)又起到了相當(dāng)重要的作用，則xj對(duì)Y的解釋能力就很大。xj的重要性程度可以用變量投影重要性指標(biāo)(variable importance in projection，VIP)來(lái)測(cè)定，如果每個(gè)xj在解釋Y時(shí)的作用都相同，則所有VIPj均等于1； 否則，對(duì)于VIPj很大(>1)的xj，在解釋Y時(shí)候就有更加重要的作用(圖9)[14-15]。

圖9 PLS模型自變量因子的重要性指標(biāo)

從350～2 500 nm譜段范圍反射率的重要性指標(biāo)可以看出，反映油膜厚度最顯著的譜段為近紅外譜段，這些譜段在油膜厚度建模過(guò)程中起到了重要的作用。

此外，自變量X與因變量Y的關(guān)系密切程度還可以通過(guò)w*c圖來(lái)展示，其中w*為構(gòu)造不同主成分分量t時(shí)各自變量的權(quán)重值，這些權(quán)值使構(gòu)造的t和u相關(guān)性最大，間接指示了X和Y的相關(guān)性大小，c為構(gòu)造主成分分量u時(shí)各變量的權(quán)重值[14]，針對(duì)主成分分量t1和t2，制作w*c1-w*c2平面圖(圖10)，可以看出油膜厚度與近紅外譜段距離最近，也就是相關(guān)性最好，這可能與原油烴類(lèi)組分和官能團(tuán)的吸收特征位于近紅外譜段有關(guān)[10]。

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

3 模型對(duì)比與評(píng)價(jià)

除了PLS模型，本次研究還選擇常用的曲線擬合方法對(duì)油膜厚度和反射率進(jìn)行建模。根據(jù)各個(gè)波段反射率和油膜厚度的相關(guān)系數(shù)，發(fā)現(xiàn)1 086 nm處的相關(guān)系數(shù)(R)最大(0.79)，以1 086 nm處的反射率為自變量，油膜厚度為因變量，分別進(jìn)行指數(shù)、對(duì)數(shù)、乘冪和多項(xiàng)式擬合，其中指數(shù)模型精度最高，決定系數(shù)(R2)達(dá)到0.918 9(圖11)。

圖11 油膜厚度-反射率曲線擬合

PLS模型和指數(shù)模型的建模精度都比較高，R2分別為0.998 0和0.918 9，但PLS模型的RMSE要明顯小于指數(shù)模型。6個(gè)驗(yàn)證樣本分別帶入建立的2個(gè)模型，PLS模型的RMSE同樣明顯小于指數(shù)模型(表2)，對(duì)比后認(rèn)為PLS模型具有更高的建模精度，而且模型的穩(wěn)定性也相對(duì)更好。

表2 不同建模方法比較

4 結(jié)論

1)油膜遙感檢測(cè)實(shí)驗(yàn)發(fā)現(xiàn)，0.05～0.60 mm厚度范圍的新鮮原油油膜，隨著油膜厚度的增加，350～1 000 nm譜段范圍的反射率不斷變小，光譜上就可以直觀反映油膜厚度變化； 而1 000～2 500 nm譜段范圍，由于反射率過(guò)小，光譜特征不明顯。

2)偏最小二乘法(PLS)適合于油膜厚度-全譜段反射率光譜數(shù)據(jù)組的回歸建模，通過(guò)最佳主成分分量確定和提取，最大程度利用了所有油膜光譜中攜帶的油膜厚度信息進(jìn)行建模。對(duì)于建模結(jié)果，還可以直觀地分析各個(gè)自變量在建模中的重要性程度，來(lái)判斷模型的合理性。研究發(fā)現(xiàn)，傳統(tǒng)油膜厚度建模時(shí)經(jīng)常被剔除的1 000～2 500 nm譜段，雖然其反射率很小，但對(duì)油膜厚度的PLS模型貢獻(xiàn)比較大，而且由于烴類(lèi)組分和官能團(tuán)的吸收特征全部位于這一譜段范圍，模型解釋更為合理，因而選擇其他方法進(jìn)行建模時(shí)，這一譜段也應(yīng)該給予關(guān)注。

3)相對(duì)于傳統(tǒng)的曲線擬合建模方法，PLS模型無(wú)論是建模樣本還是驗(yàn)證樣本的誤差均優(yōu)于傳統(tǒng)的經(jīng)驗(yàn)?zāi)Ｐ?，適合于高光譜數(shù)據(jù)的水面油膜厚度估算。