基于冬小麥高光譜圖像的天然氣微泄漏脅迫區(qū)域提取

李 輝, 劉姁升, 蔣金豹, 陳緒慧, 張 帥, 唐 珂, 趙新偉, 杜興強(qiáng), 玉龍飛雪

1. 中國(guó)四維測(cè)繪技術(shù)有限公司, 北京 100086

2. 鄂爾多斯應(yīng)用技術(shù)學(xué)院, 內(nèi)蒙古 鄂爾多斯 017000

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

4. 生態(tài)環(huán)境部衛(wèi)星環(huán)境應(yīng)用中心, 北京 100094

5. 中國(guó)資源衛(wèi)星應(yīng)用中心, 北京 100094

引 言

天然氣在能源結(jié)構(gòu)中逐漸占據(jù)重要地位, 近年來(lái)我國(guó)天然氣消費(fèi)量逐年增加, 2020年天然氣消費(fèi)量已超越3 600億m3, 在一次能源消費(fèi)比例中占據(jù)8.3%～10%的高比[1]。 天然氣儲(chǔ)氣庫(kù)儲(chǔ)存量大、 經(jīng)濟(jì)合理、 安全系數(shù)高, 在優(yōu)化供氣系統(tǒng)和保障供氣安全上建設(shè)儲(chǔ)氣庫(kù)頗為重要。 2015年我國(guó)天然氣管道建設(shè)的長(zhǎng)度約7.2萬(wàn)km[2], 已建成的地下儲(chǔ)氣庫(kù)據(jù)不完全統(tǒng)計(jì)為25座[3]。 由于地下管道和地下儲(chǔ)氣庫(kù)常年埋藏于地下, 無(wú)氧腐蝕、 自然災(zāi)害、 注入井和管道口松懈等因素會(huì)導(dǎo)致氣體發(fā)生泄漏[4]。 天然氣管線(xiàn)或儲(chǔ)氣庫(kù)泄漏引起的爆炸事故, 帶來(lái)的災(zāi)害、 傷亡嚴(yán)重威脅人類(lèi)生命、 財(cái)產(chǎn)、 環(huán)境安全。

地下儲(chǔ)存的油氣等烴類(lèi)化合物輕微泄漏會(huì)導(dǎo)致土地表面的巖石、 土壤發(fā)生裂變和腐蝕, 研究表明利用高光譜遙感研究巖石、 土壤的裂變腐蝕來(lái)探尋油氣等烴類(lèi)化合物的泄漏點(diǎn)可行[5]。 而埋藏在地下的天然氣管道或是儲(chǔ)氣庫(kù)輕微泄漏初期, 腐蝕能力較弱, 難以使地面的巖石、 土壤產(chǎn)生裂變和腐蝕。 但天然氣微泄漏會(huì)占據(jù)土壤中O2含量限制地表植被光合和呼吸作用, 進(jìn)而導(dǎo)致地面植被生長(zhǎng)異常[6]。 高光譜遙感技術(shù)以植被為媒介可間接探測(cè)地下天然氣管道泄漏點(diǎn)。

目前國(guó)內(nèi)外已有許多學(xué)者致力于利用高光譜遙感技術(shù)分析植被脅迫癥狀并間接識(shí)別油氣微泄漏點(diǎn)[7-13]。 崔鑫利用航空高光譜數(shù)據(jù), 結(jié)合野外實(shí)測(cè)光譜, 提取準(zhǔn)噶爾盆地東北緣油氣微滲漏引起的烴類(lèi)及相關(guān)蝕變礦物信息, 結(jié)果顯示, 研究區(qū)的油氣滲漏異常區(qū)具有環(huán)帶狀分布特征[11]。 Noomen等設(shè)計(jì)盆栽玉米試驗(yàn), 分別通入天然氣、 CH4、 C2H6進(jìn)行脅迫, 并利用波段深度(band depth, BD)指數(shù)來(lái)識(shí)別不同氣體泄漏脅迫下的玉米; 在實(shí)驗(yàn)中發(fā)現(xiàn), 三種氣體的泄漏會(huì)導(dǎo)致土壤中氧氣含量的減少[12]。 趙欣梅利用EO-1 Hyprion成像光譜數(shù)據(jù)以植被為媒介, 觀(guān)察植被光譜變化進(jìn)而識(shí)別油氣滲漏點(diǎn)。 但這些研究都只局限于室內(nèi), 未在野外開(kāi)展實(shí)驗(yàn)。 van der Werff等利用機(jī)載高光譜數(shù)據(jù), 通過(guò)多次霍夫變換(Hough transform)結(jié)合天然氣微泄漏的行為規(guī)律, 檢測(cè)出近似圓形的天然氣微泄漏點(diǎn)[8]; van der Werff等依據(jù)植被光譜反射率一階微分最大值在受到外界脅迫后會(huì)“紅邊藍(lán)移”, 利用Hymap成像高光譜數(shù)據(jù), 對(duì)野外管道微泄漏點(diǎn)進(jìn)行了成功探尋[13]; 李輝[7]等模擬天然氣地下管道和儲(chǔ)氣庫(kù)微泄漏對(duì)地表冬小麥的脅迫實(shí)驗(yàn), 采集冬小麥冠層光譜信息, 利用連續(xù)小波變換構(gòu)建CWTmexh指數(shù)模型, 可以較好識(shí)別受脅迫和健康冬小麥。 但現(xiàn)有研究大多利用單一的植被冠層光譜或成像光譜數(shù)據(jù), 將冠層尺度構(gòu)建的指數(shù)模型應(yīng)用到成像尺度的案例不多。

本工作擬將冬小麥冠層光譜構(gòu)建的CWTmexh指數(shù)模型應(yīng)用到成像高光譜數(shù)據(jù), 分析CWTmexh指數(shù)模型的可識(shí)別性、 穩(wěn)定性, 探尋天然氣微泄漏脅迫下的空間特征。

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

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

實(shí)驗(yàn)區(qū)位于北京市大興區(qū)長(zhǎng)子營(yíng)鎮(zhèn), 實(shí)驗(yàn)場(chǎng)地長(zhǎng)40 m, 寬20 m, 設(shè)計(jì)8個(gè)冬小麥實(shí)驗(yàn)小區(qū), 分別為4個(gè)天然氣脅迫實(shí)驗(yàn)區(qū)和4個(gè)對(duì)照區(qū), 實(shí)驗(yàn)區(qū)與對(duì)照區(qū)間隔排列, 每個(gè)區(qū)大小為2.5 m×2.5 m, 區(qū)之間的間隔均為0.5 m。 天然氣泄漏點(diǎn)位于實(shí)驗(yàn)區(qū)的中心下方60 cm處(圖1中紅色圓點(diǎn)即微泄漏位置)。 實(shí)驗(yàn)時(shí)間為2016.10—2017.06。 天然氣泄漏速率為1 L·min-1, 于2017年4月11日在小麥返青期前開(kāi)始持續(xù)通氣至6月1日實(shí)驗(yàn)結(jié)束。 種植的冬小麥品種為京冬14號(hào), 該品種成穗率較高, 抗倒伏、 抗病性較好。

圖1 實(shí)驗(yàn)區(qū)空間分布圖

圖2 (a) 實(shí)驗(yàn)田概況; (b)高光譜成像數(shù)據(jù)采集

1.2 高光譜圖像采集與預(yù)處理

實(shí)驗(yàn)采用SOC710-VP高光譜成像儀, 其波長(zhǎng)范圍為400～1 000 nm, 鏡頭焦距為8 mm, 視場(chǎng)角為21°。 采集高光譜成像數(shù)據(jù)時(shí), 為完全覆蓋實(shí)驗(yàn)田2.5 m×2.5 m的范圍, 將成像光譜儀搭載在5 m高的升降平臺(tái)上, 并保持成像光譜儀的鏡頭垂直向下。 選擇天氣晴朗, 風(fēng)力小于三級(jí), 無(wú)云或少云的時(shí)段進(jìn)行采集, 觀(guān)測(cè)時(shí)間為10：00—14：00之間, 此時(shí)太陽(yáng)高度角足夠大。 共選取4期圖像數(shù)據(jù), 分別為5月1日, 5月11日, 5月21日, 5月31日, 其中5月1日的數(shù)據(jù)是在天然氣脅迫第21天之后采集。 數(shù)據(jù)采集后用SOC710-VP配套的軟件SRAnal710進(jìn)行反射率轉(zhuǎn)換, 得到反射率圖像。 反射率圖像數(shù)據(jù)選用MNF正逆變換和5點(diǎn)平滑相結(jié)合的方法進(jìn)行降噪平滑處理。

1.3 圖像處理及分析方法

處理及分析方法有連續(xù)統(tǒng)去除、 連續(xù)小波變換、 Kapur最大熵閾值分割、 PCA及SVM分類(lèi)等。 具體處理流程如圖3所示。

圖3 基于冬小麥高光譜圖像天然氣微泄漏脅迫區(qū)域提取流程

(1)CWTmexh指數(shù)模型

利用CWTmexh指數(shù)模型提取天然氣微泄漏脅迫區(qū)域。CWTmexh指數(shù)是基于墨西哥帽(Mexihat)母小波(尺度參數(shù)32)對(duì)連續(xù)統(tǒng)去除后的冠層光譜進(jìn)行連續(xù)小波變換的特征, 具體請(qǐng)參照文獻(xiàn)[7]。

(1)

式(1)中,CW為連續(xù)小波能量系數(shù)值, 其下標(biāo)為波長(zhǎng)。

(2) Kapur最大熵閾值分割算法

利用Kapur最大熵閾值分割區(qū)分圖像的目標(biāo)和背景。 Kapur最大熵閾值熵越大, 分布越均勻。 具體原理為： 使選擇的閾值分割圖像目標(biāo)、 背景兩部分灰度統(tǒng)計(jì)的信息量為最大。

(3) 主成分分析(PCA)

利用主成分分析(PCA)對(duì)高光譜數(shù)據(jù)有效降維并減小數(shù)據(jù)冗余。 PCA前幾個(gè)主成分基本能夠概括所有波段的95%以上信息, 且每個(gè)主成分能反應(yīng)高光譜圖像的不同信息和特征[10]。

(4) 支持向量機(jī)(SVM)

利用SVM對(duì)圖像數(shù)據(jù)的PCA特征進(jìn)行監(jiān)督分類(lèi), 采用的SVM核函數(shù)為徑向基(RBF)核函數(shù)。 其核心思想是在線(xiàn)性可分的數(shù)據(jù)當(dāng)中, 尋找最優(yōu)分隔面來(lái)對(duì)數(shù)據(jù)進(jìn)行分隔, 能夠最大程度地將待分樣本分離, 且保證分隔距離最大[10]。

(5) 圖像形態(tài)學(xué)分析

利用數(shù)學(xué)形態(tài)學(xué)開(kāi)運(yùn)算達(dá)到圖像增強(qiáng)。 數(shù)學(xué)形態(tài)學(xué)分析的基本方法有腐蝕(erosion)、 膨脹(dilation)、 開(kāi)運(yùn)算(open)、 閉運(yùn)算(close)。 開(kāi)運(yùn)算在操作上具有優(yōu)勢(shì), 可以做到不明顯改變目標(biāo)面積的同時(shí)平滑目標(biāo)的邊緣。

(6) 最小二乘擬合圓

最小二乘在回歸問(wèn)題解算中用于估計(jì)和預(yù)測(cè)輸入的量和輸出變量?jī)烧咧g存在的關(guān)聯(lián)。 如圖4所示為最小二乘擬合圓的示意圖, 點(diǎn)(Xi,Yi)到圓心的距離平方和半徑的平方差為σi, 則

圖4 最小二乘擬合圓示意圖

(2)

令Q為σi的平方和, 則Q為關(guān)于A(yíng),B,R的函數(shù),A,B,R使得Q值最小時(shí)的最優(yōu)解即為所求的圓心坐標(biāo)和半徑。

(3)

在最小二乘擬合前, 先對(duì)形態(tài)學(xué)分析開(kāi)運(yùn)算處理后的冬小麥二值化圖像進(jìn)行連通區(qū)域標(biāo)記, 利用8鄰接連通區(qū)域分析, 獲取冬小麥在天然氣脅迫區(qū)域下的輪廓和邊界。

2 結(jié)果與討論

2.1 基于CWTmexh指數(shù)的天然氣微泄漏脅迫區(qū)域提取

CWTmexh指數(shù)對(duì)高光譜圖像進(jìn)行單一波段計(jì)算, 運(yùn)算得到的結(jié)果如圖5所示, 經(jīng)過(guò)CWTmexh指數(shù)運(yùn)算后, 土壤等出現(xiàn)異常值, 從灰度圖上明顯看到脅迫區(qū)域小麥較亮, 隨著天然氣脅迫的持續(xù), 脅迫范圍的暈圈在不斷擴(kuò)大。

圖5 冬小麥高光譜指數(shù)運(yùn)算結(jié)果

Kapur最大熵閾值分割對(duì)指數(shù)運(yùn)算結(jié)果進(jìn)行處理, 得到二值化圖像, 在Kapur最大熵閾值分割算法前先對(duì)灰度圖像進(jìn)行中值濾波和歸一化處理。 定義白色區(qū)域(1值)為天然氣泄漏脅迫區(qū)域, 黑色區(qū)域(0值)為其他, 如圖6所示, 可有效區(qū)分出脅迫和健康小麥區(qū)域。 整體上脅迫區(qū)域都呈現(xiàn)出暈圈狀空間特征, 隨天然氣泄漏脅迫持續(xù), 暈圈狀逐漸變大, 且逐漸趨于圓形。

圖6 Kapur最大熵閾值分割后的二值化圖像

2.2 基于SVM分類(lèi)器的天然氣微泄漏脅迫區(qū)域提取

對(duì)四期圖像數(shù)據(jù)進(jìn)行PCA降維處理, 并保留其前四個(gè)主成分PC1, PC2, PC3, PC4, 圖7所示為5月1日?qǐng)D像數(shù)據(jù)的前四個(gè)主成分的結(jié)果。 前四個(gè)主成分基本上能占據(jù)圖像99.82%的信息量, 能很好的反映原始高光譜圖像。

圖7 5月1日數(shù)據(jù)的第四個(gè)主成分： (a)PC1, (b)PC2, (c)PC3, (d)PC4

對(duì)圖像數(shù)據(jù)進(jìn)行SVM分類(lèi)。 先對(duì)每期圖像進(jìn)行訓(xùn)練樣本選取, 每期圖像分布均勻的選取受脅迫、 健康小麥訓(xùn)練樣本各40個(gè), 驗(yàn)證樣本各20個(gè), 訓(xùn)練樣本和驗(yàn)證樣本不重疊, 且驗(yàn)證樣本盡量只選擇純凈像元。 在選取訓(xùn)練樣本前, 為防止土壤對(duì)分類(lèi)結(jié)果產(chǎn)生影響, 利用掩膜濾除部分土壤。 基于樣本選取的結(jié)果, 計(jì)算兩類(lèi)樣本的J-M距離, 區(qū)分兩類(lèi)樣本的可分離性, J-M踞離可分離性原理請(qǐng)參照參考文獻(xiàn)[7]。 最終選取的訓(xùn)練樣本特征分離度如表1所示。

表1 選取的訓(xùn)練樣本特征分離度

對(duì)PCA降維后的圖像數(shù)據(jù)進(jìn)行SVM分類(lèi), 生成分類(lèi)圖像, 紅色代表受脅迫區(qū)域的小麥, 如圖8所示。 從分類(lèi)結(jié)果來(lái)看, 隨著天然氣的持續(xù)通入, 冬小麥?zhǔn)芴烊粴饷{迫區(qū)域逐漸增大, 且在空間特征上逐漸接近圓形, 與指數(shù)模型的分類(lèi)結(jié)果保持一致。 使用驗(yàn)證樣本建立混淆矩陣, 得到每期圖像的分類(lèi)精度和Kappa系數(shù), 如表2所示。 4期圖像數(shù)據(jù)分類(lèi)精度都優(yōu)于93%, Kappa系數(shù)大于0.83, 5.31日數(shù)據(jù)分類(lèi)精度達(dá)到99.25%, Kappa系數(shù)為0.97。 由表2也可以看到, 冬小麥隨著天然氣脅迫持續(xù), 受脅迫小麥的分類(lèi)精度在增大, Kappa系數(shù)在增大, 表明可分性逐漸增加。

表2 4期數(shù)據(jù)的分類(lèi)精度

圖8 SVM分類(lèi)結(jié)果

2.3 數(shù)學(xué)形態(tài)學(xué)圖像處理和最小二乘擬合圓分析

天然氣小孔泄漏脅迫區(qū)域呈暈圈狀空間特征, 隨著脅迫時(shí)間的持續(xù), 黑色暈圈在擴(kuò)大, 逐漸呈圓形分布且具有連續(xù)性, 可以通過(guò)數(shù)學(xué)形態(tài)學(xué)和連通區(qū)域分析來(lái)提取天然氣泄漏脅迫區(qū)域。 對(duì)二值化圖像進(jìn)行腐蝕、 膨脹開(kāi)運(yùn)算操作, 之后對(duì)其進(jìn)行填充、 連續(xù)邊界提取, 處理結(jié)果如圖9所示。 利用CWT和SVM提取的冬小麥?zhǔn)芴烊粴饷{迫的天數(shù)與半徑關(guān)系如表3所示。 結(jié)合天然氣泄漏行為及其空間信息, 天然氣脅迫天數(shù)和脅迫區(qū)域半徑呈一元線(xiàn)性回歸關(guān)系, 結(jié)果如圖10所示。CWTmexh指數(shù)模型下冬小麥?zhǔn)芴烊粴庑孤┟{迫天數(shù)與脅迫半徑的關(guān)系為：y=0.013x+0.492(x≥21),R2=0.97, 基于SVM分類(lèi)下冬小麥?zhǔn)芴烊粴庑孤┟{迫天數(shù)與半徑的關(guān)系為：y=0.015x+0.439(x≥21),R2=0.94, 其中x均表示脅迫天數(shù)(單位： d),y均表示天然氣脅迫半徑(單位： m)。 該結(jié)果表明冬小麥在天然氣泄漏脅迫下的半徑與泄漏時(shí)間成近似線(xiàn)性關(guān)系。

表3 不同模型下冬小麥?zhǔn)芴烊粴饷{迫的天數(shù)與半徑關(guān)系

圖9 冬小麥?zhǔn)芴烊粴饷{迫范圍擬合結(jié)果

圖10 CWTmexh指數(shù)、 SVM分類(lèi)提取出的冬小麥?zhǔn)苊{迫范圍和脅迫時(shí)間的一元回歸分析

3 結(jié) 論

通過(guò)建立野外試驗(yàn)場(chǎng), 模擬天然氣地下管道和儲(chǔ)氣庫(kù)微泄露對(duì)地表冬小麥的脅迫實(shí)驗(yàn), 基于時(shí)間序列的冬小麥高光譜圖像數(shù)據(jù), 融合其圖像、 空間、 時(shí)相特征, 探尋被天然氣脅迫的冬小麥脅迫范圍半徑和脅迫時(shí)長(zhǎng)之間的關(guān)系, 以提高并驗(yàn)證通過(guò)植被光譜信息間接檢測(cè)天然氣泄漏點(diǎn)的可識(shí)別性、 穩(wěn)定性。 主要結(jié)論如下：

基于冬小麥實(shí)驗(yàn)數(shù)據(jù),CWTmexh指數(shù)模型在圖像尺度識(shí)別天然氣泄漏特征上同樣適用;

基于SVM分類(lèi)器可以把受脅迫的冬小麥區(qū)域提取出來(lái), 分類(lèi)精度較好, 最大分類(lèi)精度可以達(dá)到99.25%, Kappa系數(shù)為0.97。 冬小麥隨著天然氣脅迫持續(xù), 受脅迫小麥的分類(lèi)精度在增大;

小麥?zhǔn)苊{迫區(qū)域半徑和通氣時(shí)間呈現(xiàn)強(qiáng)烈的線(xiàn)性相關(guān), 可以預(yù)測(cè)地下天然氣微泄漏隨著時(shí)間變化引起的脅迫區(qū)域變化。

將CWTmexh指數(shù)應(yīng)用到同期成像高光譜數(shù)據(jù), 表現(xiàn)出較好的識(shí)別性能。 雖然該實(shí)驗(yàn)研究在探索利用高光譜識(shí)別天然氣微泄漏上取得了一些進(jìn)展, 但本實(shí)驗(yàn)的高光譜圖像數(shù)據(jù)僅采用的是5 m高的地面升降平臺(tái)系統(tǒng), 而隨著衛(wèi)星數(shù)據(jù)源的不斷豐富, 星載高光譜遙感大范圍對(duì)天然氣管道和儲(chǔ)氣庫(kù)的監(jiān)測(cè)會(huì)逐漸普及, 如何利用星載高光譜圖像數(shù)據(jù)來(lái)識(shí)別天然氣微泄漏點(diǎn)還有待探究。