基于GDAL的遙感圖像變化檢測技術

蔣世豪 ，江 洪

1.福州大學 數字中國研究院（福建），福州 350108

2.福州大學 空間數據挖掘與信息共享教育部重點實驗室，衛星空間信息技術綜合應用國家地方聯合工程研究中心，福州 350108

1 引言

圖像變化檢測是指通過對同一區域不同時間觀測的圖像進行處理與分析，來確定目標變化情況的過程。遙感圖像因其覆蓋地域大的特點被廣泛應用于地物目標的變化檢測[1-4]，在各領域均具有廣泛的應用，如環境監測、土地利用、農作物生長狀況監測、森林采伐監測、災情估計等。因此，開展遙感圖像變化檢測方法的研究具有重要的應用價值[5-8]。

近年來，國內外的學者們針對不同的遙感數據源以及不同的待檢測地物特征，提出了許多快捷、有效的變化檢測方法。其中，由于較高的計算效率以及較好的適用度，直接比較法是最為常用的一種方法，其原理是：將兩幅同位置不同時相的遙感圖像進行逐像元點位的分析比對，構造出差異圖像后，按照一定的閾值，確定出發生變化的點位[9-11]。但是，構造差異圖像及閾值的判定一直是直接比較法的兩個難點。構造差異圖像需要依據待檢測地物的特征，選取敏感波段或進行波段組合來有效反映該種地物的變化狀態。利用閾值可以將差異圖像的像元分為兩類，一類是明顯發生變化的區域，另一類是未明顯發生變化的區域[6，12]。其難點是如何自動確定閾值，使變化閾值具有更好的普適性。目前，一般采用閾值分割的方法來確定變化閾值，并得到明顯發生變化的點位[13-15]。

為了從遙感圖像中提取有用信息，傳統對遙感圖像的實驗或單機模式大多是使用專業的圖像處理軟件（如ENVI、ARCGIS）等。這種方式雖然也能夠對遙感圖像進行處理和分析，但是更多地是利用軟件封裝好的方法且離不開人為重復性的輸入參數，這樣會難以避免過多的人為干預所造成的誤差[16-18]。同時，由于遙感圖像包含的信息量巨大，這種方式不能有效地擴展遙感圖像的處理算法。本文基于GDAL（Geospatial Data Abstraction Library）讀取遙感圖像的柵格像元數據并提取其對應的波段信息、坐標投影、仿射變換矩陣等元數據，結合遙感圖像處理技術和聚類分析技術等，構建了一套遙感圖像變化檢測的技術流程，實現了遙感圖像變化檢測的自動化處理，可以為區域尺度的植被長勢監測提供決策支持。

2 算法設計

2.1 技術路線

本文算法的技術流程如圖1。首先，利用GDAL 讀取遙感圖像的像元數據及空間元數據；第二步，通過波段組合，計算出能夠有效消除地形陰影影響（包括本影與落影）且能夠有效反演植被長勢信息的植被指數SEVI（Shadow Elimination Vegetation Index）；第三步，對不同時相的SEVI 結果進行標準化，以及利用圖像差值法計算出差異圖像；第四步，利用OTSU 法自動提取明顯發生變化的像元點位，該點位包括植被長勢增加的點位及植被長勢減少的點位；第五步，利用K均值聚類法提取出變化點位聚集區域的核心，并根據每一核心的所有點位的相對行列號，自動識別出變化區域；第六步，依據前幾步獲取到的檢測結果及第一步提取的空間元數據，輸出檢測結果。

圖1 技術流程圖

2.2 GDAL簡介

GDAL 是使用C/C++語言編寫的一套用于讀取空間數據的跨平臺開源庫，它利用抽象數據模型來表達所支持的各種文件格式，還有一系列命令行工具來進行數據轉換與處理。它提供了對多種柵格數據的支持，包括Arc/Info ASCII Grid（asc）、GeoTiff（tiff）、Erdas Imagine-Images（img）、ASCII DEM（dem）等格式[19]。值得一提的是，GDAL 使用抽象數據模型（Abstract Data Model）來解析它所支持的數據格式。抽象數據模型包括數據集（Dataset）、坐標系統、仿射地理坐標轉換（Affine Geo-Transform）、大地控制點（GCPs）、元數據（Metadata）、柵格波段（Raster Band）、顏色表（ColorTable）、子數據集域（Subdatasets Domain）、圖像結構域（Image_Structure Domain）、XML域（XML Domains）。

基于GDAL 開發遙感圖像算法具有明顯的優勢。這是因為遙感圖像是一種包含空間位置特征的地理空間數據，這使它區別于一般的計算機數據。普通的計算機圖像處理庫不能有效地根據空間數據結構解析出空間元數據，而使用GDAL不僅能夠快速讀取多種格式遙感圖像上的像素值，還能夠根據原有數據結構有效地提取遙感圖像的投影坐標系、仿射變換等信息，這有助于進行空間分析以及輸出包含地理坐標位置的結果圖像?；贕DAL的諸多功能及優點，利用它開發遙感圖像處理的算法，能極大地提升計算效率，簡化人為利用軟件處理圖像的步驟，也能根據不同的具體應用需求，開發不同的遙感圖像處理算法。

2.3 植被長勢參數反演

在植被光譜特征中，由于植被的葉綠素在紅色波段表現出較強的吸收能力，而植被葉片的結構特征導致在近紅色波段的輻射吸收較少，形成較強的反射。因此，以數學組合形式，構建以紅波段、近紅波段為主的植被指數，可以反映出植被的長勢情況[5-6，20]。本文采用陰影消除植被指數（SEVI）進行植被長勢信息提取，這是因為在某些復雜地形山區，太陽輻射的傳輸易受到地形效應的影響。SEVI能夠有效消除崎嶇地形中本影和落影的影響，獲得準確的植被長勢信息[21-22]，計算公式如下：

式中，Bnir為近紅外波段反射率；Bred為紅色波段反射率；Bnir/Bred為比值型植被指數RVI；1/Bred為陰影植被指數SVI；f(Δ)為地形調節因子，通過控制它來調節陰坡、陽坡的植被指數值，調節計算方法如下：

其中，x是SEVI,y1是RVI,y2是SVI,n是子區域圖像中的像素數。設定f(Δ)從0變為1，固定步長為0.001。如式（4），當兩個相關系數趨近相等時，確定為最優調節因子。

2.4 變化檢測

2.4.1 圖像差值法

圖像差值法是一種簡單快速的變化檢測方法[9-12]，它的原理是：圖像中未發生變化的地類在兩個時相的遙感圖像上一般具有相近或相等的灰度值，而當地類發生變化時，對應位置的灰度值將有較大差別。因此在差值圖像中，沒有發生變化區域的像素值接近0，明顯發生變化區域的像素值會明顯區別于0，從而使變化信息從圖像中顯現出來。

因為不同時期的SEVI 受到的地形陰影影響不同，采用了不同的地形調節因子，所以需要進行不同時期圖像的標準化處理。假設同一自然保護區同一季相的不同年份的遙感圖像的植被指數最大值相近，而在巖石具有相近的最小值。這些最大值和最小值可以用于不同年份SEVI數據之間的標準化，公式如下：

式中，SEVI′是標準化的SEVI值；SEVImin是預先判定為某塊巖石處SEVI的值；SEVImax是預先判定為某塊植被長勢最好處SEVI的值。

式中，E(t1-t2)代表從t1到t2時間段內SEVI的變化量；SEVIt1代表t1時相的植被指數；SEVIt2代表t2時相的植被指數。

2.4.2 閾值的選取

本文采用最大類間方差法（OTSU）來自動獲取閾值。這是一種被公認的閾值自動選取方法的最優方法之一，具有算法簡單、分割速度快等優點[14，23-25]。閾值的自動獲取方式較人機交互方式更加智能化，且避免了因人的主觀判斷所帶來的誤差。最大類間方差法（OTSU）的基本思想是如果一幅圖像由一物體和背景構成，物體與背景有不同的灰度值，以差異圖像中的某一灰度為閾值將圖像分為目標和背景兩組并計算兩組間的方差，當被分成的兩組之間的方差最大時，就以這個灰度值作為閾值分割圖像[25]。

假設圖像的大小為M×N，圖像中像素的灰度值滿足（大于或小于）閾值T的像素個數記作N1，像素值不滿足（大于或小于）閾值T的像素個數記作N2，則有：

式中，ω1為當前滿足閾值的像素點數占整幅圖像的比例，其平均灰度為μ1；ω2為當前不滿足閾值的像素點數占整幅圖像的比例，其平均灰度為μ2；xi、xj為遙感圖像的像素值。

式中，g為類間方差，按照圖像的灰度特征，采用遍歷的方法得到使類間方差最大的閾值T。

2.5 利用K 均值聚類分割與識別變化區域

大部分的地理現象都具有空間相關特性，即距離越近的兩事物越相似。因此，植被長勢的變化點位在一定的區域范圍內往往符合聚集規律，但是整體仍會呈現出不均勻、不規則的分布狀態，甚至在某些區域還會出現異常點。本文根據所有已檢測到點位的行列位置，利用K均值聚類算法將這些點位分為不同區域的簇類，然后再根據每個簇類中所有點位的行列位置，繪制出更具直觀效果的形狀，來輔助地物判斷與分析。

K均值聚類的基本思想是：對于給定的聚類數目k，首先隨機創建一個初始分類，然后采用迭代方法通過聚類中心的不斷移動來嘗試著改進劃分[26-30]。具體公式如下：

式中，E是所有樣本的聚類誤差；Ci是第i個簇；x是Ci中的樣本點；μi是Ci中的質心。

這種聚類算法簡單且高效，但其聚類結果的優劣往往取決于預先設定的聚類數目k。為了提高聚類穩定性并降低時間復雜度，本文選擇平均輪廓系數來自動確定聚類數目k，提高自動化程度。具體公式如下：

式中，ai是當前第i個樣本與同簇其他樣本的平均距離，稱為凝聚度；bi是當前第i個樣本與最近簇中所有樣本的平均距離，稱為分離度。Si是當前第i個樣本的個體輪廓系數，當求出所有樣本的輪廓系數后再求平均值就得到了平均輪廓系數[31]。平均輪廓系數的取值范圍為[-1,1]，且簇內樣本的距離越近，簇間樣本距離越遠，平均輪廓系數越大，聚類效果越好[32-34]。

結合以上原理及式（10）、（11）可以找出最優的聚類簇數Kopt，同時會產生Kopt個聚類簇。此時，為了識別每一個簇所在的區域，可以根據每一個簇中所有點位的行列位置來提取這些簇的外接圖形。其中，外接矩形具有較好的目視效果，且相對其他圖形比較規整，更容易看出當前簇的形狀大小、內部的點位密度等信息。因此，本文選擇繪制這些簇的外接矩形，具體方式如下：

3 實驗與結果分析

3.1 植被長勢參數反演結果

為驗證遙感圖像變化檢測技術的效果，本文選擇了Landsat8-OLI 的兩幅遙感圖像作為數據源（表1），裁剪邊界為福建省武夷山自然保護區（圖2），該保護區內地形崎嶇復雜，最小坡度為0°，平均坡度為27.3°，最大坡度大于60°，標準方差為10.04。地物類型主要是植被，同時包含了一些水體、裸地和建筑物。

依據式（1）～（5）得到了上述兩幅遙感圖像的SEVI標準化結果（圖3），從目視效果來看，SEVI 的標準化結果，整體呈現出平坦的特征，基本沒有圖2 中山體區域出現的浮雕情況。這與Jiang 等[21]發現的SEVI 呈現的目視效果一致。

3.2 差值檢測結果

依據式（6）～（9）得到了兩幅遙感圖像的變化檢測結果（圖4 和表2）。已知類別1 是植被長勢明顯減少的區域，在研究區內占120個像元（0.108 km2）；類別2是植被長勢明顯增加的區域，在研究區內占1 659 個像元（1.490 km2），是明顯減少區域的13 倍左右。兩種明顯發生變化的像元點位共約占所有像元數的0.10%。

表1 兩幅Landsat8-OLI遙感圖像的相關參數

圖2 Landsat8-OLI遙感圖像

圖3 SEVI的標準化結果

圖4 變化點位

表2 OTSU法閾值檢測結果

從圖4可以發現，植被長勢發生明顯變化的地區主要集中于研究區中部偏北區域，呈現小面積聚集狀態。

3.3 利用K 均值聚類的識別結果

根據圖像差值以及OTSU 閾值提取的明顯變化點位進行聚類分析，結果如圖5所示。

圖5 變化點位的行列號

根據研究區內變化點位的數量以及分布情況（圖5），設定聚類簇數K以1 為步長，從1 到60，計算當前聚類簇數K所對應的平均輪廓系數（式（11）），并繪制兩者的關系曲線（圖6）?？梢园l現，當聚類簇數為15 的時候，類別1 有最高的平均輪廓系數為0.649 7。而類別2的曲線呈雙峰分布且右邊的峰值比左邊的峰值更大，所以取對應的聚類簇數為56。

圖6 平均輪廓系數

依據確定出的Kopt個聚類簇以及每簇中所有點位的行列號，結合式（12）、（13），就可以繪制出所有聚類簇的外接矩形（圖7）。通過與圖4的對比可以發現，圖7具有更明顯的目視效果，提高了檢測精度。

圖7 識別結果

3.4 不同區域的測試結果

為了充分驗證本文技術的有效性，選擇福建省閩江源國家級自然保護區作為測試研究區（圖8）。該保護區內的地形同樣崎嶇復雜，海拔在250～1 858 m之間，平均坡度為14.2°，標準方差為8.45，包含了植被、裸地、建筑物、陰影這些主要地物類別，相比武夷山自然保護區，它的地物分布情況較為錯綜復雜。選擇的遙感圖像詳細信息可見表3。

圖8 Landsat8-OLI遙感圖像（圖像是753波段的假彩色合成）

表3 兩幅Landsat8-OLI遙感圖像的相關參數

圖9 結果展示

通過本文技術檢測和識別的結果（圖9 和表4），可以得出以下結論：（1）植被長勢明顯增加的點位主要分布在兩個區域，分別是閩江源的最北與最南部地區，面積較大，共占2 260個像元，約2.034 0 km2。（2）植被長勢明顯減少的點位比較少，僅占25 個像元，分布在6 個不同的區域，共約0.022 5 km2。（3）圖9（b）的目視效果要優于圖9（a），更能直觀地看出植被長勢增加的有兩塊區域，且南部區域的密度要整體大于北部區域，而且突出了減少點位區域的目視效果（如需要更加突出，則僅需繪制時加粗線條）。

表4 閩江源國家級自然保護區兩幅遙感圖像的檢測結果

4 結束語

針對遙感圖像變化檢測的自動化程度較低，現有算法不易擴展等問題，本文提出一套基于GDAL的遙感圖像變化檢測技術，主要包括：利用陰影消除植被指數（SEVI）反演植被長勢；利用圖像差值法及最大類間方差法（OTSU）來提取植被長勢明顯變化點位；利用K均值聚類自動分割并識別變化區域。實驗結果表明，本文方法能夠有效地檢測植被長勢變化區域，可為區域尺度的植被長勢監測提供技術支持。下一步將研究更多閾值自動化提取算法、分割算法在遙感圖像變化檢測中的應用。