融合光譜-紋理-顏色信息識(shí)別地下天然氣微泄漏點(diǎn)

吳自勇，蔣金豹，郭鑒威，王鑫達(dá)，季 楊

中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083

引 言

天然氣是目前應(yīng)用最為廣泛的清潔能源之一。中國(guó)的天然氣消費(fèi)量由2001年的274億m3增長(zhǎng)至2017年的2 373億m3，年均增長(zhǎng)14.4%，增幅約為同期世界平均水平的7倍[1]。在此背景下，加快天然氣地下儲(chǔ)氣庫(kù)的建設(shè)成為必然選擇，但是儲(chǔ)氣庫(kù)在運(yùn)行過(guò)程中因受到地質(zhì)構(gòu)造、腐蝕以及施工等多種因素的影響會(huì)發(fā)生泄漏，泄漏的天然氣會(huì)通過(guò)土壤緩慢滲透到大氣中，不僅會(huì)對(duì)地表土壤以及植被產(chǎn)生影響，而且當(dāng)泄漏的天然氣濃度達(dá)到一定程度后還會(huì)引發(fā)爆炸、人員凍傷、窒息等事故，對(duì)生命財(cái)產(chǎn)安全造成巨大損失。此外，天然氣雖作為一種清潔能源，在未來(lái)仍有可能會(huì)成為溫室效應(yīng)的重要來(lái)源[2]。目前檢測(cè)地下儲(chǔ)氣庫(kù)泄漏的方法存在效率低、技術(shù)要求高、效果差的缺陷[3]，因此，迫切需要一種高效準(zhǔn)確的檢測(cè)方法。由于泄漏的天然氣會(huì)致使地表相應(yīng)位置的植被產(chǎn)生脅迫反應(yīng)，因此，可通過(guò)探究植被的脅迫反應(yīng)間接地識(shí)別地下儲(chǔ)氣庫(kù)泄漏點(diǎn)，提高檢測(cè)效率，節(jié)約檢測(cè)成本，避免安全事故的發(fā)生。

目前，高光譜技術(shù)在天然氣脅迫下植被的識(shí)別研究方面多集中于應(yīng)用單一光譜信息進(jìn)行目標(biāo)識(shí)別[4-5]。Smith等[6]研究發(fā)現(xiàn)小麥、草地在天然氣泄漏脅迫下一階微分光譜在紅邊內(nèi)逐漸減小，并伴有多峰與藍(lán)移現(xiàn)象，并且利用725與702 nm一階微分比值指數(shù)可以識(shí)別天然氣脅迫下的草地、小麥和大豆； Noomen等[7]發(fā)現(xiàn)指數(shù)R440/R740能探測(cè)僅有的4個(gè)天然氣泄漏點(diǎn)； Lassalle等[4]發(fā)現(xiàn)受油氣脅迫61 d后的植被光譜曲線在可見(jiàn)光、近紅外區(qū)域均有所升高。這些研究中僅利用植被的光譜信息進(jìn)行識(shí)別，雖取得一定的效果，卻忽略了空間信息。

成像光譜儀能夠同時(shí)獲取待檢測(cè)對(duì)象的光譜和空間信息，這在天然氣微泄漏點(diǎn)檢測(cè)中具有重要意義，能夠高效精確地檢測(cè)出天然氣泄漏點(diǎn)的空間位置。Noomen等[7]在美國(guó)Ojai地區(qū)發(fā)現(xiàn)天然氣長(zhǎng)期泄漏導(dǎo)致地表草地出現(xiàn)一個(gè)半徑為20 m的環(huán)形裸土區(qū)，裸土區(qū)外則是茂密的草地，在考慮天然氣泄漏脅迫植被癥狀空間特征的基礎(chǔ)上，融合光譜信息成功探測(cè)天然氣泄漏點(diǎn)。Werff等[8]發(fā)現(xiàn)地下天然氣泄漏導(dǎo)致地表直徑約6～14 m的環(huán)狀區(qū)域內(nèi)植被出現(xiàn)了明顯的脅迫癥狀，并嘗試?yán)肏ough變換提取空間信息，融合光譜分類信息探測(cè)天然氣泄漏脅迫植被。本文運(yùn)用高光譜成像技術(shù)基于野外模擬實(shí)驗(yàn)檢測(cè)天然氣泄漏點(diǎn)。

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

1.1 實(shí)驗(yàn)設(shè)計(jì)

實(shí)驗(yàn)區(qū)位于北京市大興區(qū)(地理位置： 39°39′2.56″ N，116°34′33.10″ E)。實(shí)驗(yàn)區(qū)長(zhǎng)約40 m，寬約20 m，區(qū)域內(nèi)共設(shè)置24個(gè)邊長(zhǎng)2.5 m的方形地塊，其中12個(gè)地塊為通氣脅迫組，剩余12個(gè)地塊為控制組，地塊之間間隔0.5 m，地塊下埋藏有天然氣管道，管道泄漏口在地塊中心以下0.6 m處。自2018年10月1日起，通過(guò)天然氣控制系統(tǒng)以1 L·min-1的速度連續(xù)不斷通氣，截止到2018年11月18日。實(shí)驗(yàn)期內(nèi)每天進(jìn)行巡視，保證植被不受其他因素(如缺肥、缺水)的影響。

1.2 數(shù)據(jù)采集

使用SOC710VP成像光譜儀在天氣狀況優(yōu)良，當(dāng)?shù)貢r(shí)間10:00—14:00時(shí)間段內(nèi)，獲取各地塊的高光譜影像。冬小麥影像數(shù)據(jù)采集時(shí)間分別為10月11日、10月24日、11月1日、11月9日和11月18日。模型驗(yàn)證所選的數(shù)據(jù)分別為大豆、玉米和草地，所選取的大豆影像為8月24日，玉米影像為8月7日，草地影像為9月17日。

1.3 影像處理

1.3.1 預(yù)處理

(1)標(biāo)準(zhǔn)反射率校正

SOC710VP采集的影像數(shù)據(jù)并非地物反射率數(shù)據(jù)，而是原始DN值，需經(jīng)過(guò)輻射校正處理，將原始DN值轉(zhuǎn)換為反射率數(shù)據(jù)，本文采用SOC710VP的配套軟件SRAnal710對(duì)影像數(shù)據(jù)校正。

(2)光譜平滑

對(duì)校正后的影像使用五點(diǎn)移動(dòng)均值平滑進(jìn)行平滑處理，平滑公式為[9]

(1)

式(1)中，Rs為每個(gè)像元平滑后的光譜，Ri第i波段的反射率值，且i∈[3，k-2]k為高光譜圖像的波段數(shù)。其他影像操作還包括影像剪裁。

1.3.2 后處理

(1)波段篩選

高光譜圖像能夠提供眾多精細(xì)波段，但是信息冗余的問(wèn)題也會(huì)凸顯，因此有必要篩選出包含有用信息的特征波段[10]。蔣金豹等對(duì)高光譜數(shù)據(jù)進(jìn)行ANOVA處理，選擇敏感波段構(gòu)建光譜指數(shù)來(lái)識(shí)別不同品質(zhì)的建筑涂料。本文運(yùn)用單因素方差分析(one-way ANOVA)評(píng)價(jià)脅迫與非脅迫小麥的波段可分性，選取p<0.05置信水平下的有效波段，去除可分性較差的波段[11]。

(2)紋理及顏色特征提取

天然氣泄漏脅迫下作物長(zhǎng)勢(shì)及生理特征都會(huì)發(fā)生變化，這些變化不僅體現(xiàn)在光譜特征上還體現(xiàn)在紋理特征上，因此，可利用紋理特征輔助識(shí)別天然氣脅迫下的植被[12]?；诮y(tǒng)計(jì)學(xué)的灰度共生矩陣(GLCM)在紋理提取中應(yīng)用較為廣泛，本文采用GLCM計(jì)算紋理特征，灰度共生矩陣的距離參數(shù)設(shè)置為1，分別計(jì)算0°，45°，90°和135°方向上的同質(zhì)性、對(duì)比度、非相似性、熵值、能量和相關(guān)性。此外，考慮到受天然氣脅迫的植被會(huì)變黃，本文利用顏色矩計(jì)算RGB合成影像的前三階顏色矩[13]，分別為一階矩(Mean，μi)、二階矩(Variance，σi)和三階矩(Skewness，si)。

(2)

1.4 建模方法

1.4.1 歸一化植被指數(shù)(NDVI)

NDVI是目前應(yīng)用最廣泛的植被指數(shù)之一，在已有研究中，該指數(shù)常被用來(lái)區(qū)分植被和裸土區(qū)域[14]。本文應(yīng)用NDVI設(shè)置閾值提取脅迫下的裸露土壤，為防止閾值設(shè)置不合理影響后期數(shù)據(jù)處理，本實(shí)驗(yàn)通過(guò)多期影像進(jìn)行結(jié)果對(duì)比，并最終設(shè)置閾值為0.4。

1.4.2 最小二乘支持向量機(jī)

支持向量機(jī)(SVM)在高光譜影像分類中應(yīng)用廣泛，最小二乘支持向量機(jī)(least squares support vector machine，LSSVM)作為一種新型支持向量機(jī)方法，用等式約束代替了傳統(tǒng)支持向量機(jī)中的不等式約束，用誤差平方替換了SVM中的松弛變量，通過(guò)求解一組等式方程得出最優(yōu)分類超平面，從而避開(kāi)了求解計(jì)算量相對(duì)繁重的二次規(guī)劃問(wèn)題，有效提高了求解速度[15]，LSSVM模型的輸出可表示為

(3)

式中，ai對(duì)應(yīng)xi的Lagrange乘子，對(duì)于核函數(shù)k(x)，本文選用高斯徑向基(RBF)。

1.4.3 模型建立

本文首先基于NDVI對(duì)影像閾值分割，對(duì)于植被部分，提取光譜、紋理以及顏色特征，運(yùn)用LSSVM識(shí)別脅迫下的植被區(qū)域，并與裸露土壤部分融合，得到識(shí)別結(jié)果，模型可簡(jiǎn)單表示如式(4)

(4)

其中，G(x)為脅迫區(qū)域，X為光譜、紋理以及顏色特征向量所構(gòu)成的矩陣，SNDVI<0.4為裸露土壤區(qū)域。

為更好提取出天然氣泄漏脅迫區(qū)域，本文運(yùn)用圖像形態(tài)學(xué)處理中腐蝕、膨脹和填充等算法對(duì)提取結(jié)果進(jìn)行處理，并利用最小二乘擬合方法對(duì)識(shí)別的脅迫區(qū)域邊界進(jìn)行擬合。

2 結(jié)果與討論

2.1 波段篩選與特征提取結(jié)果

對(duì)預(yù)處理后的影像進(jìn)行ANOVA，在p<0.05水平下，選取F較大的波段，并剔除300 nm附近噪聲較大的波段，優(yōu)選出510，520，570，625，645，680和690 nm七個(gè)敏感波段，如圖1所示。

由圖2可知，脅迫第11天，受脅迫的小麥與為控制組小麥差異不顯著。脅迫第24天之后，脅迫組小麥變黃，長(zhǎng)勢(shì)稀疏，脅迫組小麥光譜反射率在可見(jiàn)光范圍內(nèi)(370～760 nm)增大，而近紅外范圍(760～1 100 nm)降低； 脅迫組影像的同質(zhì)性、能量和相關(guān)性三個(gè)紋理特征均低于控制組，對(duì)比度、非相似性和熵值均高于控制組； 脅迫組的影像的一階矩和三階矩明顯高于控制組。隨著脅迫的持續(xù)進(jìn)行，受脅迫的小麥干枯、凋落，脅迫組小麥近紅外區(qū)域光譜反射率仍低于控制組，但是與控制組差距減小。

圖1 波段篩選Fig.1 Band Selection

2.2 模型識(shí)別

2.2.1 識(shí)別結(jié)果

針對(duì)植被的識(shí)別，將上述所提取的光譜、紋理以及顏色特征向量歸一化作為L(zhǎng)SSVM的輸入，運(yùn)用MATLAB編程，LSSVM中懲罰系數(shù)C=100，RBF參數(shù)γ=20，系數(shù)采用10-fold交叉驗(yàn)證得到； 對(duì)于裸露土壤，計(jì)算影像的NDVI，將NDVI值小于0.4的區(qū)域作為裸露土壤區(qū)域。將植被的識(shí)別結(jié)果與裸露土壤的識(shí)別結(jié)果融合，并對(duì)識(shí)別結(jié)果分別運(yùn)用2×2模板進(jìn)行腐蝕、膨脹處理，提取脅迫邊界，得到如圖3的結(jié)果。

圖3 (a) 原始影像； (b) LSSVM識(shí)別的天然氣脅迫下的小麥； (c) 裸露土壤； (d) 融合結(jié)果腐蝕膨脹操作后的結(jié)果； (e) 提取的脅迫邊界和脅迫區(qū)域圓擬合的邊界

2.2.2 結(jié)果分析

地塊持續(xù)通氣11 d后，脅迫組小麥僅有少量葉片邊緣有微黃跡象，脅迫癥狀不明顯，并且小麥生長(zhǎng)初期易受到干旱、養(yǎng)分缺失等其他條件的相互作用，致使脅迫區(qū)與非脅迫區(qū)小麥在光譜、紋理以及顏色特征方面差異不顯著，故模型未能較好地識(shí)別出脅迫第11天的區(qū)域。隨脅迫的持續(xù)進(jìn)行，脅迫中心區(qū)域小麥逐漸變黃，并且多數(shù)并在脅迫后期枯萎，與對(duì)照區(qū)相比，脅迫中心區(qū)小麥葉片萎縮，植株矮小、稀疏，模型均可較好地識(shí)別出脅迫區(qū)域，由圖4的識(shí)別結(jié)果可以看出，脅迫發(fā)生第11天至第40天，脅迫區(qū)域逐漸變大，而第40天之后，脅迫區(qū)域趨于穩(wěn)定，通過(guò)表1中對(duì)4個(gè)地塊的脅迫區(qū)域半徑統(tǒng)計(jì)發(fā)現(xiàn)，脅迫第40天至第49天，脅迫區(qū)域有縮小的趨勢(shì)，通過(guò)進(jìn)一步觀測(cè)分析可以看出，在脅迫與非脅迫的交界處，隨著小麥的生長(zhǎng)，未受到脅迫的小麥或受脅迫較輕微的小麥葉片會(huì)逐漸伸長(zhǎng)至脅迫區(qū)域，這種現(xiàn)象映射到冠層的影像中，則會(huì)致使脅迫區(qū)域的縮小。

2.3 討論

2.3.1 模型驗(yàn)證

為驗(yàn)證本文所構(gòu)建的模型能否適用于其他類型植被，本文運(yùn)用構(gòu)建的模型對(duì)同一試驗(yàn)區(qū)的大豆、玉米和草地進(jìn)行分析，均得到較好地識(shí)別結(jié)果，如圖5所示。

表1 脅迫天數(shù)與脅迫半徑關(guān)系Table 1 Relationship between stress days and stress radius

2.3.2 綠暈現(xiàn)象

通過(guò)觀察影像可以發(fā)現(xiàn)，在距脅迫區(qū)0～100像素(0.25 m)的范圍內(nèi)，小麥長(zhǎng)勢(shì)較邊緣對(duì)照區(qū)域旺盛，這一現(xiàn)象與Noomen所發(fā)現(xiàn)的“green halo”現(xiàn)象較為一致，Noomen結(jié)合前人的研究推斷油氣泄漏會(huì)將水汽帶到地表并且在暈圈附近聚集，同時(shí)泄漏的氣體會(huì)使周圍的氣體分布發(fā)生變化，在兩者的作用下，將會(huì)出現(xiàn)綠暈現(xiàn)象[7]。隨著時(shí)間的推移，該區(qū)域的小麥葉片會(huì)向脅迫區(qū)伸進(jìn)，這也進(jìn)一步解釋了脅迫后期脅迫區(qū)域有縮小趨勢(shì)。

圖5 模型驗(yàn)證Fig.5 Model Validation

圖6 綠暈現(xiàn)象Fig.6 Green halo phenomenon

3 結(jié) 論

通過(guò)野外模擬天然氣微泄漏對(duì)地表植被的脅迫實(shí)驗(yàn)，分析不同時(shí)間序列天然氣微泄漏脅迫下冬小麥的高光譜影像數(shù)據(jù)，基于LSSVM并結(jié)合光譜、紋理以及顏色信息，構(gòu)建天然氣微泄漏識(shí)別模型，得到如下結(jié)論：

(1) 天然氣泄漏脅迫下的冬小麥與對(duì)照區(qū)小麥相比，可見(jiàn)光范圍光譜反射率上升，近紅外區(qū)域下降； 通過(guò)方差分析選擇510，520，570，625，645，680和690 nm七個(gè)可分性較好的特征波段，對(duì)特征波段分析可知，受脅迫的小麥影像同質(zhì)性、能量和相關(guān)性三個(gè)紋理特征較對(duì)照區(qū)降低，而對(duì)比度、非相似性和熵值三個(gè)紋理特征較對(duì)照區(qū)提高； 脅迫下的冬小麥影像前三階顏色矩總體上高于控制組。

(2) 本文的模型融合了光譜、紋理和顏色信息，并且在脅迫第24天有效識(shí)別出天然氣泄漏區(qū)域。模型同樣適用于大豆、玉米和草地，具有一定的適用性。

(3) 通過(guò)對(duì)脅迫區(qū)域圓擬合可知，冬小麥?zhǔn)芴烊粴饷{迫的范圍呈現(xiàn)先擴(kuò)大后縮小的時(shí)序特征。

本文構(gòu)建的模型通過(guò)識(shí)別植被的變化能夠較好的識(shí)別天然氣微泄漏，但模型仍有不足之處，比如混合像元、陰影以及脅迫區(qū)域周圍土壤的存在都會(huì)對(duì)模型的識(shí)別效果產(chǎn)生一定的影響。在接下來(lái)的研究中，將會(huì)結(jié)合天然氣疑似泄漏區(qū)的多期影像、天然氣管道埋藏圖、地質(zhì)圖等資料解決這些問(wèn)題，及時(shí)高效地識(shí)別出天然氣泄漏點(diǎn)，為以后的工程應(yīng)用提供理論基礎(chǔ)和技術(shù)支持。