基于GDAL的遙感圖像變化檢測(cè)技術(shù)

蔣世豪 ，江 洪

1.福州大學(xué) 數(shù)字中國(guó)研究院（福建），福州 350108

2.福州大學(xué) 空間數(shù)據(jù)挖掘與信息共享教育部重點(diǎn)實(shí)驗(yàn)室，衛(wèi)星空間信息技術(shù)綜合應(yīng)用國(guó)家地方聯(lián)合工程研究中心，福州 350108

1 引言

圖像變化檢測(cè)是指通過(guò)對(duì)同一區(qū)域不同時(shí)間觀測(cè)的圖像進(jìn)行處理與分析，來(lái)確定目標(biāo)變化情況的過(guò)程。遙感圖像因其覆蓋地域大的特點(diǎn)被廣泛應(yīng)用于地物目標(biāo)的變化檢測(cè)[1-4]，在各領(lǐng)域均具有廣泛的應(yīng)用，如環(huán)境監(jiān)測(cè)、土地利用、農(nóng)作物生長(zhǎng)狀況監(jiān)測(cè)、森林采伐監(jiān)測(cè)、災(zāi)情估計(jì)等。因此，開(kāi)展遙感圖像變化檢測(cè)方法的研究具有重要的應(yīng)用價(jià)值[5-8]。

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

為了從遙感圖像中提取有用信息，傳統(tǒng)對(duì)遙感圖像的實(shí)驗(yàn)或單機(jī)模式大多是使用專業(yè)的圖像處理軟件（如ENVI、ARCGIS）等。這種方式雖然也能夠?qū)b感圖像進(jìn)行處理和分析，但是更多地是利用軟件封裝好的方法且離不開(kāi)人為重復(fù)性的輸入?yún)?shù)，這樣會(huì)難以避免過(guò)多的人為干預(yù)所造成的誤差[16-18]。同時(shí)，由于遙感圖像包含的信息量巨大，這種方式不能有效地?cái)U(kuò)展遙感圖像的處理算法。本文基于GDAL（Geospatial Data Abstraction Library）讀取遙感圖像的柵格像元數(shù)據(jù)并提取其對(duì)應(yīng)的波段信息、坐標(biāo)投影、仿射變換矩陣等元數(shù)據(jù)，結(jié)合遙感圖像處理技術(shù)和聚類分析技術(shù)等，構(gòu)建了一套遙感圖像變化檢測(cè)的技術(shù)流程，實(shí)現(xiàn)了遙感圖像變化檢測(cè)的自動(dòng)化處理，可以為區(qū)域尺度的植被長(zhǎng)勢(shì)監(jiān)測(cè)提供決策支持。

2 算法設(shè)計(jì)

2.1 技術(shù)路線

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

圖1 技術(shù)流程圖

2.2 GDAL簡(jiǎn)介

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

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

2.3 植被長(zhǎng)勢(shì)參數(shù)反演

在植被光譜特征中，由于植被的葉綠素在紅色波段表現(xiàn)出較強(qiáng)的吸收能力，而植被葉片的結(jié)構(gòu)特征導(dǎo)致在近紅色波段的輻射吸收較少，形成較強(qiáng)的反射。因此，以數(shù)學(xué)組合形式，構(gòu)建以紅波段、近紅波段為主的植被指數(shù)，可以反映出植被的長(zhǎng)勢(shì)情況[5-6，20]。本文采用陰影消除植被指數(shù)（SEVI）進(jìn)行植被長(zhǎng)勢(shì)信息提取，這是因?yàn)樵谀承?fù)雜地形山區(qū)，太陽(yáng)輻射的傳輸易受到地形效應(yīng)的影響。SEVI能夠有效消除崎嶇地形中本影和落影的影響，獲得準(zhǔn)確的植被長(zhǎng)勢(shì)信息[21-22]，計(jì)算公式如下：

式中，Bnir為近紅外波段反射率；Bred為紅色波段反射率；Bnir/Bred為比值型植被指數(shù)RVI；1/Bred為陰影植被指數(shù)SVI；f(Δ)為地形調(diào)節(jié)因子，通過(guò)控制它來(lái)調(diào)節(jié)陰坡、陽(yáng)坡的植被指數(shù)值，調(diào)節(jié)計(jì)算方法如下：

其中，x是SEVI,y1是RVI,y2是SVI,n是子區(qū)域圖像中的像素?cái)?shù)。設(shè)定f(Δ)從0變?yōu)?，固定步長(zhǎng)為0.001。如式（4），當(dāng)兩個(gè)相關(guān)系數(shù)趨近相等時(shí)，確定為最優(yōu)調(diào)節(jié)因子。

2.4 變化檢測(cè)

2.4.1 圖像差值法

圖像差值法是一種簡(jiǎn)單快速的變化檢測(cè)方法[9-12]，它的原理是：圖像中未發(fā)生變化的地類在兩個(gè)時(shí)相的遙感圖像上一般具有相近或相等的灰度值，而當(dāng)?shù)仡惏l(fā)生變化時(shí)，對(duì)應(yīng)位置的灰度值將有較大差別。因此在差值圖像中，沒(méi)有發(fā)生變化區(qū)域的像素值接近0，明顯發(fā)生變化區(qū)域的像素值會(huì)明顯區(qū)別于0，從而使變化信息從圖像中顯現(xiàn)出來(lái)。

因?yàn)椴煌瑫r(shí)期的SEVI 受到的地形陰影影響不同，采用了不同的地形調(diào)節(jié)因子，所以需要進(jìn)行不同時(shí)期圖像的標(biāo)準(zhǔn)化處理。假設(shè)同一自然保護(hù)區(qū)同一季相的不同年份的遙感圖像的植被指數(shù)最大值相近，而在巖石具有相近的最小值。這些最大值和最小值可以用于不同年份SEVI數(shù)據(jù)之間的標(biāo)準(zhǔn)化，公式如下：

式中，SEVI′是標(biāo)準(zhǔn)化的SEVI值；SEVImin是預(yù)先判定為某塊巖石處SEVI的值；SEVImax是預(yù)先判定為某塊植被長(zhǎng)勢(shì)最好處SEVI的值。

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

2.4.2 閾值的選取

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

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

式中，ω1為當(dāng)前滿足閾值的像素點(diǎn)數(shù)占整幅圖像的比例，其平均灰度為μ1；ω2為當(dāng)前不滿足閾值的像素點(diǎn)數(shù)占整幅圖像的比例，其平均灰度為μ2；xi、xj為遙感圖像的像素值。

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

2.5 利用K 均值聚類分割與識(shí)別變化區(qū)域

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

K均值聚類的基本思想是：對(duì)于給定的聚類數(shù)目k，首先隨機(jī)創(chuàng)建一個(gè)初始分類，然后采用迭代方法通過(guò)聚類中心的不斷移動(dòng)來(lái)嘗試著改進(jìn)劃分[26-30]。具體公式如下：

式中，E是所有樣本的聚類誤差；Ci是第i個(gè)簇；x是Ci中的樣本點(diǎn)；μi是Ci中的質(zhì)心。

這種聚類算法簡(jiǎn)單且高效，但其聚類結(jié)果的優(yōu)劣往往取決于預(yù)先設(shè)定的聚類數(shù)目k。為了提高聚類穩(wěn)定性并降低時(shí)間復(fù)雜度，本文選擇平均輪廓系數(shù)來(lái)自動(dòng)確定聚類數(shù)目k，提高自動(dòng)化程度。具體公式如下：

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

結(jié)合以上原理及式（10）、（11）可以找出最優(yōu)的聚類簇?cái)?shù)Kopt，同時(shí)會(huì)產(chǎn)生Kopt個(gè)聚類簇。此時(shí)，為了識(shí)別每一個(gè)簇所在的區(qū)域，可以根據(jù)每一個(gè)簇中所有點(diǎn)位的行列位置來(lái)提取這些簇的外接圖形。其中，外接矩形具有較好的目視效果，且相對(duì)其他圖形比較規(guī)整，更容易看出當(dāng)前簇的形狀大小、內(nèi)部的點(diǎn)位密度等信息。因此，本文選擇繪制這些簇的外接矩形，具體方式如下：

3 實(shí)驗(yàn)與結(jié)果分析

3.1 植被長(zhǎng)勢(shì)參數(shù)反演結(jié)果

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

依據(jù)式（1）～（5）得到了上述兩幅遙感圖像的SEVI標(biāo)準(zhǔn)化結(jié)果（圖3），從目視效果來(lái)看，SEVI 的標(biāo)準(zhǔn)化結(jié)果，整體呈現(xiàn)出平坦的特征，基本沒(méi)有圖2 中山體區(qū)域出現(xiàn)的浮雕情況。這與Jiang 等[21]發(fā)現(xiàn)的SEVI 呈現(xiàn)的目視效果一致。

3.2 差值檢測(cè)結(jié)果

依據(jù)式（6）～（9）得到了兩幅遙感圖像的變化檢測(cè)結(jié)果（圖4 和表2）。已知類別1 是植被長(zhǎng)勢(shì)明顯減少的區(qū)域，在研究區(qū)內(nèi)占120個(gè)像元（0.108 km2）；類別2是植被長(zhǎng)勢(shì)明顯增加的區(qū)域，在研究區(qū)內(nèi)占1 659 個(gè)像元（1.490 km2），是明顯減少區(qū)域的13 倍左右。兩種明顯發(fā)生變化的像元點(diǎn)位共約占所有像元數(shù)的0.10%。

表1 兩幅Landsat8-OLI遙感圖像的相關(guān)參數(shù)

圖2 Landsat8-OLI遙感圖像

圖3 SEVI的標(biāo)準(zhǔn)化結(jié)果

圖4 變化點(diǎn)位

表2 OTSU法閾值檢測(cè)結(jié)果

從圖4可以發(fā)現(xiàn)，植被長(zhǎng)勢(shì)發(fā)生明顯變化的地區(qū)主要集中于研究區(qū)中部偏北區(qū)域，呈現(xiàn)小面積聚集狀態(tài)。

3.3 利用K 均值聚類的識(shí)別結(jié)果

根據(jù)圖像差值以及OTSU 閾值提取的明顯變化點(diǎn)位進(jìn)行聚類分析，結(jié)果如圖5所示。

圖5 變化點(diǎn)位的行列號(hào)

根據(jù)研究區(qū)內(nèi)變化點(diǎn)位的數(shù)量以及分布情況（圖5），設(shè)定聚類簇?cái)?shù)K以1 為步長(zhǎng)，從1 到60，計(jì)算當(dāng)前聚類簇?cái)?shù)K所對(duì)應(yīng)的平均輪廓系數(shù)（式（11）），并繪制兩者的關(guān)系曲線（圖6）?？梢园l(fā)現(xiàn)，當(dāng)聚類簇?cái)?shù)為15 的時(shí)候，類別1 有最高的平均輪廓系數(shù)為0.649 7。而類別2的曲線呈雙峰分布且右邊的峰值比左邊的峰值更大，所以取對(duì)應(yīng)的聚類簇?cái)?shù)為56。

圖6 平均輪廓系數(shù)

依據(jù)確定出的Kopt個(gè)聚類簇以及每簇中所有點(diǎn)位的行列號(hào)，結(jié)合式（12）、（13），就可以繪制出所有聚類簇的外接矩形（圖7）。通過(guò)與圖4的對(duì)比可以發(fā)現(xiàn)，圖7具有更明顯的目視效果，提高了檢測(cè)精度。

圖7 識(shí)別結(jié)果

3.4 不同區(qū)域的測(cè)試結(jié)果

為了充分驗(yàn)證本文技術(shù)的有效性，選擇福建省閩江源國(guó)家級(jí)自然保護(hù)區(qū)作為測(cè)試研究區(qū)（圖8）。該保護(hù)區(qū)內(nèi)的地形同樣崎嶇復(fù)雜，海拔在250～1 858 m之間，平均坡度為14.2°，標(biāo)準(zhǔn)方差為8.45，包含了植被、裸地、建筑物、陰影這些主要地物類別，相比武夷山自然保護(hù)區(qū)，它的地物分布情況較為錯(cuò)綜復(fù)雜。選擇的遙感圖像詳細(xì)信息可見(jiàn)表3。

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

表3 兩幅Landsat8-OLI遙感圖像的相關(guān)參數(shù)

圖9 結(jié)果展示

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

表4 閩江源國(guó)家級(jí)自然保護(hù)區(qū)兩幅遙感圖像的檢測(cè)結(jié)果

4 結(jié)束語(yǔ)

針對(duì)遙感圖像變化檢測(cè)的自動(dòng)化程度較低，現(xiàn)有算法不易擴(kuò)展等問(wèn)題，本文提出一套基于GDAL的遙感圖像變化檢測(cè)技術(shù)，主要包括：利用陰影消除植被指數(shù)（SEVI）反演植被長(zhǎng)勢(shì)；利用圖像差值法及最大類間方差法（OTSU）來(lái)提取植被長(zhǎng)勢(shì)明顯變化點(diǎn)位；利用K均值聚類自動(dòng)分割并識(shí)別變化區(qū)域。實(shí)驗(yàn)結(jié)果表明，本文方法能夠有效地檢測(cè)植被長(zhǎng)勢(shì)變化區(qū)域，可為區(qū)域尺度的植被長(zhǎng)勢(shì)監(jiān)測(cè)提供技術(shù)支持。下一步將研究更多閾值自動(dòng)化提取算法、分割算法在遙感圖像變化檢測(cè)中的應(yīng)用。