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

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

1. 中國四維測繪技術有限公司, 北京 100086

2. 鄂爾多斯應用技術學院, 內蒙古 鄂爾多斯 017000

3. 中國礦業大學(北京)地球科學與測繪工程學院, 北京 100083

4. 生態環境部衛星環境應用中心, 北京 100094

5. 中國資源衛星應用中心, 北京 100094

引 言

天然氣在能源結構中逐漸占據重要地位, 近年來我國天然氣消費量逐年增加, 2020年天然氣消費量已超越3 600億m3, 在一次能源消費比例中占據8.3%～10%的高比[1]。 天然氣儲氣庫儲存量大、 經濟合理、 安全系數高, 在優化供氣系統和保障供氣安全上建設儲氣庫頗為重要。 2015年我國天然氣管道建設的長度約7.2萬km[2], 已建成的地下儲氣庫據不完全統計為25座[3]。 由于地下管道和地下儲氣庫常年埋藏于地下, 無氧腐蝕、 自然災害、 注入井和管道口松懈等因素會導致氣體發生泄漏[4]。 天然氣管線或儲氣庫泄漏引起的爆炸事故, 帶來的災害、 傷亡嚴重威脅人類生命、 財產、 環境安全。

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

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

本工作擬將冬小麥冠層光譜構建的CWTmexh指數模型應用到成像高光譜數據, 分析CWTmexh指數模型的可識別性、 穩定性, 探尋天然氣微泄漏脅迫下的空間特征。

1 實驗部分

1.1 實驗設計

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

圖1 實驗區空間分布圖

圖2 (a) 實驗田概況; (b)高光譜成像數據采集

1.2 高光譜圖像采集與預處理

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

1.3 圖像處理及分析方法

處理及分析方法有連續統去除、 連續小波變換、 Kapur最大熵閾值分割、 PCA及SVM分類等。 具體處理流程如圖3所示。

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

(1)CWTmexh指數模型

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

(1)

式(1)中,CW為連續小波能量系數值, 其下標為波長。

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

利用Kapur最大熵閾值分割區分圖像的目標和背景。 Kapur最大熵閾值熵越大, 分布越均勻。 具體原理為： 使選擇的閾值分割圖像目標、 背景兩部分灰度統計的信息量為最大。

(3) 主成分分析(PCA)

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

(4) 支持向量機(SVM)

利用SVM對圖像數據的PCA特征進行監督分類, 采用的SVM核函數為徑向基(RBF)核函數。 其核心思想是在線性可分的數據當中, 尋找最優分隔面來對數據進行分隔, 能夠最大程度地將待分樣本分離, 且保證分隔距離最大[10]。

(5) 圖像形態學分析

利用數學形態學開運算達到圖像增強。 數學形態學分析的基本方法有腐蝕(erosion)、 膨脹(dilation)、 開運算(open)、 閉運算(close)。 開運算在操作上具有優勢, 可以做到不明顯改變目標面積的同時平滑目標的邊緣。

(6) 最小二乘擬合圓

最小二乘在回歸問題解算中用于估計和預測輸入的量和輸出變量兩者之間存在的關聯。 如圖4所示為最小二乘擬合圓的示意圖, 點(Xi,Yi)到圓心的距離平方和半徑的平方差為σi, 則

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

(2)

令Q為σi的平方和, 則Q為關于A,B,R的函數,A,B,R使得Q值最小時的最優解即為所求的圓心坐標和半徑。

(3)

在最小二乘擬合前, 先對形態學分析開運算處理后的冬小麥二值化圖像進行連通區域標記, 利用8鄰接連通區域分析, 獲取冬小麥在天然氣脅迫區域下的輪廓和邊界。

2 結果與討論

2.1 基于CWTmexh指數的天然氣微泄漏脅迫區域提取

CWTmexh指數對高光譜圖像進行單一波段計算, 運算得到的結果如圖5所示, 經過CWTmexh指數運算后, 土壤等出現異常值, 從灰度圖上明顯看到脅迫區域小麥較亮, 隨著天然氣脅迫的持續, 脅迫范圍的暈圈在不斷擴大。

圖5 冬小麥高光譜指數運算結果

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

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

2.2 基于SVM分類器的天然氣微泄漏脅迫區域提取

對四期圖像數據進行PCA降維處理, 并保留其前四個主成分PC1, PC2, PC3, PC4, 圖7所示為5月1日圖像數據的前四個主成分的結果。 前四個主成分基本上能占據圖像99.82%的信息量, 能很好的反映原始高光譜圖像。

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

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

表1 選取的訓練樣本特征分離度

對PCA降維后的圖像數據進行SVM分類, 生成分類圖像, 紅色代表受脅迫區域的小麥, 如圖8所示。 從分類結果來看, 隨著天然氣的持續通入, 冬小麥受天然氣脅迫區域逐漸增大, 且在空間特征上逐漸接近圓形, 與指數模型的分類結果保持一致。 使用驗證樣本建立混淆矩陣, 得到每期圖像的分類精度和Kappa系數, 如表2所示。 4期圖像數據分類精度都優于93%, Kappa系數大于0.83, 5.31日數據分類精度達到99.25%, Kappa系數為0.97。 由表2也可以看到, 冬小麥隨著天然氣脅迫持續, 受脅迫小麥的分類精度在增大, Kappa系數在增大, 表明可分性逐漸增加。

表2 4期數據的分類精度

圖8 SVM分類結果

2.3 數學形態學圖像處理和最小二乘擬合圓分析

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

表3 不同模型下冬小麥受天然氣脅迫的天數與半徑關系

圖9 冬小麥受天然氣脅迫范圍擬合結果

圖10 CWTmexh指數、 SVM分類提取出的冬小麥受脅迫范圍和脅迫時間的一元回歸分析

3 結 論

通過建立野外試驗場, 模擬天然氣地下管道和儲氣庫微泄露對地表冬小麥的脅迫實驗, 基于時間序列的冬小麥高光譜圖像數據, 融合其圖像、 空間、 時相特征, 探尋被天然氣脅迫的冬小麥脅迫范圍半徑和脅迫時長之間的關系, 以提高并驗證通過植被光譜信息間接檢測天然氣泄漏點的可識別性、 穩定性。 主要結論如下：

基于冬小麥實驗數據,CWTmexh指數模型在圖像尺度識別天然氣泄漏特征上同樣適用;

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

小麥受脅迫區域半徑和通氣時間呈現強烈的線性相關, 可以預測地下天然氣微泄漏隨著時間變化引起的脅迫區域變化。

將CWTmexh指數應用到同期成像高光譜數據, 表現出較好的識別性能。 雖然該實驗研究在探索利用高光譜識別天然氣微泄漏上取得了一些進展, 但本實驗的高光譜圖像數據僅采用的是5 m高的地面升降平臺系統, 而隨著衛星數據源的不斷豐富, 星載高光譜遙感大范圍對天然氣管道和儲氣庫的監測會逐漸普及, 如何利用星載高光譜圖像數據來識別天然氣微泄漏點還有待探究。