王曉東，何 浩，侯 東，孫冠楠，朱文泉

(1.北京師范大學資源學院，北京 100875;2.北京師范大學地表過程與資源生態(tài)國家重點實驗室，北京 100875)

基于TM圖像的農(nóng)業(yè)區(qū)域植被覆蓋變化檢測

王曉東1，2，何 浩1，2，侯 東1，2，孫冠楠1，2，朱文泉1，2

(1.北京師范大學資源學院，北京 100875;2.北京師范大學地表過程與資源生態(tài)國家重點實驗室，北京 100875)

以交叉相關(guān)光譜匹配(cross correlogram spectral matching，CCSM)為基礎構(gòu)建土地覆蓋變化強度指標，利用華北農(nóng)業(yè)植被覆蓋區(qū)2期不同時相的TM圖像計算該地區(qū)土地覆蓋變化強度圖像。認為變化強度圖像任意二階鄰域中像素的變化強度服從隱馬爾可夫模型，用馬爾可夫隨機場－最大后驗估計(maxium a posteriori estimation of markov random field，MRF－MAP)的方法從變化強度圖像中提取植被變化區(qū)域。實驗證明:該方法能夠有效識別各種外源噪聲造成的農(nóng)業(yè)植被覆蓋區(qū)域同物異譜的現(xiàn)象，可準確提取植被變化區(qū)域;但對于水體區(qū)域存在誤判現(xiàn)象。

植被覆蓋;變化檢測;交叉相關(guān)光譜匹配(CCSM);交叉相關(guān)系數(shù);隱馬爾可夫隨機場模型

0 引言

農(nóng)業(yè)植被覆蓋區(qū)域變化檢測對于及時掌握農(nóng)作物種植狀況、了解區(qū)域土地利用變化具有重要意義。衛(wèi)星傳感器獲取的遙感圖像包含有大量的地物特征信息。利用同一地區(qū)、不同時期的遙感圖像來確定地物狀態(tài)變化的過程稱為變化檢測。目前已經(jīng)提出的遙感影像變化檢測算法主要包括分類后比較法和直接比較法［1－4］。直接比較法是通過定量評價不同時相、同一像素的光譜曲線獲得不同時期地表土地覆蓋變化的強度，進而分割變化強度圖像提取變化區(qū)域。由于能夠綜合利用變化像元的強度信息和空間信息，直接比較法成為變化檢測的一個重要研究方向。

直接比較法包括檢測指標和圖像分割兩個關(guān)鍵技術(shù)［5］。當前的檢測指標主要從不同時期地物光譜在數(shù)值、曲線形狀、相似性及歸屬概率等特征上的差異出發(fā)進行設計，包括歐式距離［6］、夾角余弦［7］、相關(guān)系數(shù)［8］、交叉相關(guān)系數(shù)［9］、變化向量分析［6］、變化概率［10］等。目前用于圖像分割［11］和描述遙感圖像數(shù)據(jù)相關(guān)關(guān)系的統(tǒng)計模型包括有限高斯混合模型［12］、馬爾可夫模型［13］、隱馬爾可夫模型［14］、模糊馬爾可夫模型［15］等。隱馬爾可夫模型將有限高斯混合模型和馬爾可夫隨機場模型兩者的優(yōu)點相結(jié)合，是目前比較理想的差異圖像先驗模型。

華北農(nóng)作物種植區(qū)田塊面積較大、種植規(guī)模成片、種類相對單一，而TM遙感數(shù)據(jù)具有分辨率適中、單景覆蓋面積大的特點，基本能夠滿足這種農(nóng)作物種植區(qū)變化區(qū)域提取的需要。農(nóng)作物種植區(qū)土地覆蓋的光譜受到較大的人為干擾，農(nóng)作物輪作造成的種類差異、灌溉造成的土壤背景影響以及農(nóng)作物本身的物候特征，都導致了植被覆蓋區(qū)較多的同物異譜現(xiàn)象，而TM遙感數(shù)據(jù)較高的光譜分辨率能夠較好地反映地表植被覆蓋的連續(xù)的光譜曲線，為準確提取變化區(qū)域提供了光譜信息保障。

本研究以2期TM遙感圖像為數(shù)據(jù)源，由于華北農(nóng)作物種植區(qū)的植被光譜受到各種干擾，導致各波段光譜曲線形狀相同而光譜值統(tǒng)一偏高或偏低。針對這一特點，本文選用可比較光譜形狀差異的交叉相關(guān)系數(shù)為地表變化強度的指標;利用隱馬爾可夫模型來描述八鄰域像素中各像素之間的影響;按照貝葉斯最小規(guī)則，根據(jù)區(qū)域變化強度灰度圖像來估計類別場最大后驗概率，然后自動分割圖像，提取植被變化區(qū)域。

1 研究方法

1.1 變化檢測指標設計

檢測指標的性能是變化檢測的關(guān)鍵。就遙感圖像來說，檢測指標要能夠準確反映地表覆蓋的微小變化，同時要抑制遙感圖像獲取過程中大氣和傳感器的噪聲以及農(nóng)田灌溉作業(yè)等地表農(nóng)事活動對土地覆蓋光譜造成的干擾。

本研究以Van der Meer［9］提出的交叉相關(guān)系數(shù)作為檢測指標的設計基礎。交叉相關(guān)系數(shù)反映兩條光譜曲線在不同匹配位置上的相似度，其定義為

式中:a={ai|i∈K}，b={bi|i∈K}，a 和 b 為需要比較的兩條光譜曲線數(shù)據(jù)，K為波段集合數(shù);Rm為b光譜曲線移動m位置后和a光譜曲線的交叉相關(guān)系數(shù)，對b光譜曲線進行移動變換，其中m為b光譜曲線的整體移動程度，m=0表示沒有移動，m=1表示b曲線整體向右移動一個波段，并將最后一個波段補在第一波段的前面;i為a光譜和移動變換后的b光譜中的波段位置;n為移動后兩曲線的重疊波段數(shù)。

所有匹配位置的交叉相關(guān)系數(shù)可以組成一條差異曲線，即R={Rm|m∈K}。當兩條比較光譜采用同一圖像的光譜時，差異曲線反映土地覆蓋完全不變情況下的標準曲線;當兩條比較光譜采用不同時相的光譜時，則差異曲線反映兩個時相之間的土地覆蓋現(xiàn)狀變化的實際曲線。通過比較標準曲線和實際曲線的差別，可以刻畫土地覆蓋變化的強度。越接近標準曲線，土地覆蓋變化情況越小，反之，土地覆蓋的變化程度越大。土地覆蓋變化的檢測指標可以通過實際曲線與標準曲線之間的均方根差(RMS)進行表達［16］，即

式中Rm，R'm分別代表標準曲線和實際曲線。顯然，土地覆蓋的變化必然導致地物光譜曲線的變化，進一步引起不同時相實際曲線和標準曲線的差異，并最終反映到RMS的數(shù)值上。因此，RMS綜合代表了不同時相土地覆蓋變化的強度。

1.2 隱馬爾可夫隨機場模型

差異灰度圖像可以表示為一個變化強度的集合Y={y1，…，yi，…，yN|yi∈D，i∈S}，是可以觀測到的隨機場。其中，S={1，2，…，N}，表示所有像素的位置;D表示所有檢測指標值。令L代表類別集，則L={0，1}，其中0代表沒有變化，1代表有變化。圖像分割就是對圖像集合Y中每個像素i進行判斷并標記為L中的一個類別(l)的過程，所有像素的標記結(jié)果就是類別圖像?？梢员硎緸轭悇e的集合X={x1，…，xi，…，xN|xi∈L，i∈S}，是無法直接觀測到的隱隨機場，也是變化檢測的結(jié)果［17］。

根據(jù)有限高斯模型，圖像的概率可以用幾個分量組成的混合模型來表示，即

式中:xi和yi分別代表類別圖像集合X和強度圖像集合Y中的元素，i∈S;θ為類別參數(shù);p(xi=l)為混合比例參數(shù)，由類別數(shù)量決定，由灰度圖像統(tǒng)計特征決定，沒有考慮不同像素的空間屬性。

馬爾可夫隨機場(markov random field，MRF)很好地描述了當前像素與其鄰域中像素之間的相互關(guān)系［18］。當前點的鄰域定義為 N={Ni，i∈S}，其中鄰域不包括當前點本身i?Ni，并且鄰域關(guān)系相互對稱，即i∈Nj?j∈Ni。MRF表達了當前像素點的灰度僅與其鄰域像素的灰度有關(guān)，與其他位置像素的灰度無關(guān)，即

有限高斯模型通過馬爾可夫隨機場加入空間信息就可以得到隱馬爾可夫隨機場模型，即

式中P(xi=l|xNi)為混合比例參數(shù)，由像元在類別場X中的鄰域xNi來決定。隱馬爾可夫模型不僅能夠描述灰度圖像統(tǒng)計特征，還可以描述圖像的空間結(jié)構(gòu)，產(chǎn)生連續(xù)的空間聚類，因而可以去除單獨的噪聲，得到更準確的分類結(jié)果。

1.3 MRF－MAP分類器

基于灰度場的條件獨立性，可以得到

式中θ為類別參數(shù)，且 θ=(μxi，σxi)在遞歸過程中獲得并更新。根據(jù)Hammersley－Clifford定理，MRF分布可以等效地刻畫為Gibbs隨機場分布［19］，即

式中:U(X)是能量函數(shù);Z為配分函數(shù)，即

將式(8)代入式(7)后取對數(shù)并推導，可以得到

關(guān)系式。其中，令

則類別場X的后驗概率最大式(6)等價于后驗勢能的最小，即

1.4 檢測流程

本流程的實現(xiàn)分為4步:①分別對兩期TM圖像進行輻射校正和幾何糾正等預處理，獲取位置準確的地表覆蓋表觀反射率;②逐像素計算交叉相關(guān)系數(shù)的均方根差，得到圖像范圍內(nèi)的土地覆蓋變化強度圖像;③假設土地覆蓋變化強度和空間關(guān)系符合隱馬爾可夫模型，采用MRF－MAP分類器逐像素判斷土地類型是否變化，其中變化強度圖像的初始分類采用K均值算法，后驗勢能最小化采用條件迭代模式(iterated conditional mode，ICM)［20］，隱馬爾可夫模型的隱含參數(shù)模擬采用期望最大值算法(expection－maxium，EM)［17］;④計算整體后驗勢能，與 MRF－MAP分類前的整體后驗勢能進行比較，如果2次分類前后整體后驗勢能未達到穩(wěn)定，重新進行逐像素MRF－MAP分類，如果2次分類前后整體后驗勢能穩(wěn)定，則停止計算，輸出變化區(qū)域提取結(jié)果。整個變化檢測流程如圖1所示。

圖1 植被覆蓋區(qū)域變化檢測流程Fig.1 Flow chart of vegetation cover change detection

2 應用實例

2.1 研究區(qū)概況與數(shù)據(jù)源

研究區(qū)范圍在 E 115.95°～116.06°，N 39.67°～39.59°之間，位于北京市大興區(qū)內(nèi)。土地覆蓋類型包括耕地、林地、草地、裸地和城鎮(zhèn)用地。為驗證本文方法的實際效果，采用2006年4月7日和2006年12月3日TM圖像進行試驗。圖像時相選擇在冬小麥覆蓋區(qū)域光譜與土壤背景對比較大的抽穗和分蘗階段，以消除冬小麥季節(jié)變化對檢測結(jié)果的影響。

2.2 數(shù)據(jù)預處理和變化區(qū)域提取

對遙感圖像進行輻射校正和幾何糾正，去除與本研究無關(guān)的熱紅外波段，利用FLAASH大氣校正模塊，去除氣溶膠和視場角的影響。將其他6個波段的DN值轉(zhuǎn)化為地表反射率［21］，然后選擇二次多項式進行幾何糾正，幾何誤差控制在一個像素之內(nèi)。兩期TM地表反射率圖像如圖2所示。

圖2 研究區(qū)TM地表反射率圖像Fig.2 TM surface reflection images in study area

采用1.1節(jié)所述的變化檢測指標逐像素進行計算，獲得研究區(qū)土地覆蓋變化強度圖，如圖3(a)所示。以此為基礎，用1.3節(jié)所述的MRF－MAP分類算法進行變化區(qū)域提取，結(jié)果如圖3(b)所示。

圖3 變化強度圖像及MRF－MAP提取結(jié)果Fig.3 Change intensity images and MRF －MAP extract results

2.3 提取結(jié)果定性對比

為了更有效地對比分析提取結(jié)果，在實驗區(qū)選擇了有代表性的4個典型區(qū)域，對比了典型區(qū)域的原始地表反射率、變化強度圖像和變化區(qū)域MRFMAP提取結(jié)果，如表1所示。

表1 典型區(qū)域?qū)Ρ萒ab.1 Comparison of typical region

2.3.1 檢測指標性能分析

區(qū)域A前后兩期土地覆蓋為裸地和植被，在圖像光譜上體現(xiàn)出較大的差異，在變化強度圖上灰度值較高，說明針對土地類型不同、光譜曲線差異的土地覆蓋變化，本文采用的檢測指標能夠準確地指示出該變化的強度。

區(qū)域C前后兩期土地覆蓋主要為植被，實際并未發(fā)生變化，但后時相由于進行了灌溉作業(yè)，土壤背景發(fā)生變化，從而導致了地物在圖像上的色調(diào)也出現(xiàn)差異。在變化強度圖上區(qū)域C的灰度值較區(qū)域A灰度低，與城鎮(zhèn)區(qū)域B的灰度值相似，說明本文檢測指標對于同類地物中由于噪聲導致的光譜曲線差異并不敏感，能有效反映地表覆蓋的變化情況。

區(qū)域D前后兩期土地覆蓋為水體，實際未發(fā)生變化。在變化強度圖上區(qū)域D的灰度值較高，說明檢測指標對于水體的覆蓋變化存在較大誤差。原因是本文檢測指標主要對各波段按指定順序形成的光譜曲線形狀差異程度進行評價，而水體區(qū)域在各波段上反射率都很低，大氣輻射的影響體現(xiàn)得尤為明顯，況且微量的大氣輻射就能導致水體光譜曲線形狀發(fā)生變化，因此區(qū)域D水體形成了較高的變化強度。

2.3.2 MRF－MAP提取結(jié)果分析

區(qū)域A在后時相冬小麥種植區(qū)域內(nèi)存在線狀的田間道路，與前時相土地覆蓋類型一致，變化強度圖中線狀道路灰度值低，準確體現(xiàn)了地表覆蓋未變化的情況。MRF－MAP提取方法在判斷像素變化與否的時候考慮了相鄰像素的情況，由于田間道路周圍都是變化的像素，所以MRF－MAP將田間道路誤認為是變化區(qū)域中的噪聲，并判斷為變化區(qū)域。說明MRF－MAP方法在去除點狀噪聲的同時，沒有很好地保留呈現(xiàn)線狀特征的不變區(qū)域。

區(qū)域B是土地覆蓋類型基本未發(fā)生變化的城鎮(zhèn)區(qū)域，但該區(qū)域中存在植被綠化等現(xiàn)象，因此零星覆蓋類型發(fā)生的變化導致了變化強度圖中該區(qū)域存在較多的點狀變化區(qū)。由于區(qū)域B中的變化區(qū)域周圍為未變化的房屋或道路，綜合鄰域像素的信息后，MRF－MAP法將點狀變化區(qū)域判斷為未變化，與實際情況一致。

2.4 提取結(jié)果精度評價

為了進一步定量評價本文方法的有效性，在全圖范圍內(nèi)隨機選擇600個點作為評價總體。通過目視解譯獲取評價總體的真實變化情況，并根據(jù)本方法的提取結(jié)果分別計算總體精度、總體Kappa系數(shù)、用戶精度和生產(chǎn)者精度4個評價指標，結(jié)果如表2所示。

表2 提取結(jié)果精度評價Tab.2 Accuracy evaluation of the extract results

本文方法總體Kappa系數(shù)為0.40，說明本文方法的提取結(jié)果和真實情況接近中等吻合的程度?？傮w精度達到0.81，說明在以像素為基本單位的精度評價體系下，本方法的精度還存在較大的提升空間，主要錯分出現(xiàn)在線狀區(qū)域和水體區(qū)域。

3 結(jié)論

本文從比較不同時期土地覆蓋變化導致的地物光譜曲線形狀差異出發(fā)，以交叉相關(guān)系數(shù)為基礎構(gòu)建檢測指標，對土地覆蓋的變化強度進行表達，生成土地覆蓋變化強度圖像;再以隱馬爾可夫模型為基礎，利用變化強度和變化區(qū)域空間分布信息對變化強度圖像進行圖像分割，提取變化區(qū)域。本方法初步實現(xiàn)了變化區(qū)域的自動提取，能夠避免傳統(tǒng)方法的閾值選取和偽像素人工去除兩大問題，提高了變化區(qū)域提取的效率，為快速、準確獲取農(nóng)作物種植區(qū)植被覆蓋變化區(qū)域提供了借鑒。通過定性分析和精度評價，可以得出如下結(jié)論:

1)檢測指標能夠避免地面灌溉和不同時期農(nóng)區(qū)植被長勢差異造成的干擾，可準確反映農(nóng)區(qū)土地覆蓋的變化強度;

2)能夠有效利用土地覆蓋變化的空間信息，抑制城鎮(zhèn)和農(nóng)田中的點狀噪聲對提取結(jié)果的干擾;

3)整個過程無需先驗知識和人工交互，實現(xiàn)了自動化;

4)由于本方法是以地物光譜曲線形狀的差異來衡量土地覆蓋變化強度，而水體本身的反射率接近于零，很小的干擾就能引起光譜曲線形狀的變化，因此本文檢測指標對于水體區(qū)域的檢測結(jié)果誤差較大。今后需在提高檢測指標在水體區(qū)域的適應性或者對未變化水體區(qū)域進行掩模等方面開展進一步研究。

［1］Lu D，Mausel P，Brondizio E，et al.Change Detection Techniques［J］.International Journal of Remote Sensing，2004，25(12):2365 －2401.

［2］李淑坤，李培軍，程 濤.加入多時相紋理的遙感變化檢測［J］.國土資源遙感，2009(3):35－40.

［3］謝仁偉，牛 錚，孫 睿，等.基于多波段統(tǒng)計檢驗的土地利用變化檢測［J］.國土資源遙感，2009(2):66－70.

［4］Foody G M.Monitoring the Magnitude of Land － cover Change Around the Southern Limits of the Sahara［J］.Photogrammetric Engineering and Remote Sensing，2001，67(7):841 －847.

［5］李月臣，陳 晉，宮 鵬，等.基于NDVI時間序列數(shù)據(jù)的土地覆蓋變化檢測指標設計［J］.應用基礎與工程科學學報，2005，13(3):261－275.

［6］陳 晉，何春陽，史培軍，等.基于變化向量分析的土地利用/覆蓋變化動態(tài)監(jiān)測(I)——變化閾值的確定方法［J］.遙感學報，2001，5(4):259 －266.

［7］許衛(wèi)東，尹 球，匡定波.地物光譜匹配模型比較研究［J］.紅外與毫米波學報，2005，24(4):296 －300.

［8］申邵洪，萬幼川，龔 浩，等.遙感影像變化檢測自適應閾值分割的Kriging方法［J］.武漢大學學報:信息科學版，2009，34(8):902－905.

［9］Van der Meer.Spectral Curve Shape Matching with a Continuum Removed CCSM Algorithm［J］.International Journal of Remote Sensing，2000，21(16):3179 －3185.

［10］Chen X H，Jin C，Shen M G，et al.Land － use/Land － cover Change Detection Using Change－vector Analysis in Posterior Probability Space［C］//Proc.SPIE 7144，714405，2008:10.1117/12.812671.

［11］劉永學，李滿春，毛 亮.基于邊緣的多光譜遙感圖像分割方法［J］.遙感學報，2006，10(3):350 －356.

［12］Bilmes J A.A Gentle Tutorial of the EM Algorithm and Its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models［R］.Technical Report:ICSI.TR － 97 －02，International Computer Science Institute，Berkeley CA，USA，1998.

［13］劉偉強，陳 鴻，夏德深.基于馬爾可夫隨機場的遙感圖像分割和描述［J］.東南大學學報:自然科學版，1999，29(11):11－15.

［14］鐘家強，王潤生.基于自適應參數(shù)估計的多時相遙感圖像變化檢測［J］.測繪學報，2005，34(4):331 －336.

［15］鄭 瑋，康戈文，陳武凡，等.基于模糊馬爾可夫隨機場的無監(jiān)督遙感圖像分割算法［J］.遙感學報，2008，12(2):246 －252.

［16］Wang L，Chen J，Gong P，et al.Land Cover Change Detection with a Cross－ correlogram Spectral Matching Algorithm［J］.International Journal of Remote Sensing，2009，30(12):3259 －3273.

［17］Zhang Y，Brady M，Smith S.Segmentation of Brain MR Images Through a Hidden Markov Random Field Model and the Expectation－ maximization Algorithm［J］.IEEE Transactions on Medical Image，2001，20(1):45 －57.

［18］Liu D S，Kelly M，Gong P.A Spatial－temporal Approach to Monitoring Forest Disease Spread Using Multi－temporal High Spatial Resolution Imagery［J］.Remote Sensing of Environment，2006，101(2):167－180.

［19］馮衍秋，梁 斌，陳 明，等.基于gibbs隨機場的有限混合模型改進與腦部MR圖像的穩(wěn)健分割［J］.中國生物醫(yī)學工程學報，2003，22(3):193 －198.

［20］Besag J.On the Statistical Analysis of Dirty Pictures［J］.Journal of the Royal Statistical Society:Series B(Methodological)，1986，48(3):259－302.

［21］Matthew M W，Adler Gdden S M，Berk A，et al.Atmospheric Correction ofSpectralImagery:Evaluation oftheFLAASH Algorithm with AVIRIS Data［C］//Presented at SPIE Proceeding，Algorithm and Technologies for Multispectral，Hyperspectral，and Ultraspectral Imagery IX，2003.

Vegetation Cover Change Detection in the Cropping Area Based on TM Image

WANG Xiao － dong1，2，HE Hao1，2，HOU Dong1，2，SUN Guan － nan1，2，ZHU Wen － quan1，2
(1.College of Resource Sciences＆ Technology，Beijing Normal University，Beijing 100875，China;2.State Key Laboratory of Earth Surface Processes and Resource Ecology，Beijing Normal University，Beijing 100875，China)

In this paper，a change intensity indicator of land cover based on cross correlogram spectral matching(CCSM)technique was employed to generate the change intensity image of the cropping vegetation cover area in North China between two TM images in different periods.It was first considered that the change intensity of image pixel of the two－order neighbor in the change intensity image obeyed the hidden markov random field model，and then the vegetation cover change area was extracted from the change intensity image using maximum a posteriori estimation of markov random field(MRF － MAP)model.The experiment has proved that the proposed method could precisely extract vegetation cover change and inhibit effectively the same object with different spectra due to exogenous noises in the cropping vegetation cover area.However，this method seems to perform unsatisfactorily over the water area.

vegetation cover;change detection;cross correlogram spectral matching(CCSM);cross－correlation coefficient;hidden markov random field model

TP 79

A

1001－070X(2012)02－0092－06

王曉東(1983－)，男，博士，主要從事統(tǒng)計遙感及GIS軟件工程方面研究。E－mail:wangxd22@gmail.com。

何 浩(1980－)，男，博士研究生，主要從事資源環(huán)境遙感及3S應用技術(shù)研究。E－mail:hehaowhu@163.com。

(責任編輯:邢 宇)

10.6046/gtzyyg.2012.02.17

2011－07－28;

2011－09－09

國家高技術(shù)研究發(fā)展計劃項目(編號:2006AA120101)和國家自然科學基金項目(編號:40871191)共同資助。