水體動態(tài)變化信息遙感監(jiān)測方法研究:以秦皇島桃林口水庫為例

崔名赫，潘軍，蔣立軍，劉祉余，許鶴，麻洪川

1.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026； 2.中國人民解放軍31441部隊，沈陽 110001； 3.長春建筑學(xué)院 交通學(xué)院，長春 130699； 4.中國人民解放軍31009部隊，北京 100088

0 引言

遙感動態(tài)變化監(jiān)測技術(shù)是基于同一地區(qū)的多時相遙感數(shù)據(jù)，對地表地物類型的轉(zhuǎn)化/改變進行識別、屬性變化信息進行獲取以及遙感制圖與分析的一種方法。動態(tài)變化監(jiān)測作為遙感圖像分析中的一項重要應(yīng)用，為環(huán)境監(jiān)測、土地利用和資源勘探等方面提供了有效的技術(shù)手段[1--4]。

為探究有效的動態(tài)監(jiān)測方法，相關(guān)學(xué)者提出了不同的研究方法。例如，李丹等[5]以遙感圖像解譯技術(shù)作為支撐，將遙感影像中目標(biāo)地物的各種影像特征(光譜特征，紋理特征和形狀特征等)為標(biāo)志，運用4個年份的多時相Landsat遙感影像實現(xiàn)對耕地資源空間信息的提取，分析研究區(qū)耕地分布及變化。Li et al.[6]對北京市1984—2013年的城市變化狀況進行了分類與變化檢測，通過目視解譯選取了準(zhǔn)確、穩(wěn)定不變的訓(xùn)練樣本，隨后利用該訓(xùn)練樣本進行了分類，并將分類結(jié)果使用到變化檢測中。喬芳敏等[7]采用最大似然法對影像進行分類，再進行分類后比較，進而得出變化信息。紀(jì)超南[8]采用數(shù)據(jù)分組和直方圖匹配算法對長時序影像進行分類，進而提高了變化檢測的精度。

目前已有的研究方法大體上可分為兩類：直接比較法和分類后比較法。前者是直接利用多時相遙感影像/數(shù)據(jù)對地物類型的光譜、空間變化特征進行增強與分類；后者是在不同時相的遙感影像上單獨進行分類，然后比較圖像同一位置的分類結(jié)果，進而確定變化信息的位置和類型。分類后比較的方法不必考慮不同時相圖像間輻射校正和匹配等問題，但是它的精度主要取決于每個時相圖像的分類精度，存在嚴(yán)重的誤差累積問題。直接比較法大致分為以下幾類：①彩色合成法：將后一時相某一波段賦予R通道、前一時相同一波段同時賦予G、B通道，變化信息一般表現(xiàn)為紅色、青色色彩，易于識別。該方法操作簡單，但精度較低，適于初步分析。②光譜特征變異法：在多時相遙感數(shù)據(jù)中，變化像元的光譜特征發(fā)生了顯著變化，由此可計算兩個時相數(shù)據(jù)中同一像元樣本之間的歐式距離、夾角余弦等參數(shù)來作為光譜特征發(fā)生變化的度量指標(biāo)，但不同地物類型變化可能會有相同或類似的度量值。③差值法、比值法，即：將兩個時相的遙感數(shù)據(jù)對應(yīng)波段進行差值、比值波段運算，則變化信息在差值、比值數(shù)據(jù)中具有顯著區(qū)別于未變化像元的數(shù)據(jù)特征。與上述“光譜特征變異法”同理，不同的變化信息可能會具有相同或相似的“差值光譜特性曲線”或“比值光譜特性曲線”，此種特殊情況會導(dǎo)致信息混淆。

隨著多時相圖像數(shù)據(jù)源可獲取性的增加，人們對遙感圖像變化檢測及趨勢分析的要求也在不斷提高。首先希望高精度獲取變化結(jié)果并盡可能減少預(yù)處理誤差;其次希望了解變化過程,獲得變化區(qū)域更豐富的屬性信息;最后希望充分挖掘多時相圖像的時序信息,并全面結(jié)合變化檢測結(jié)果,對未知時相進行趨勢預(yù)測分析[9]。一個好的變化檢測方法應(yīng)該回答3個方面的問題：變化的區(qū)域和程度、變化類型的分布以及變化的類型[10]。因此，如何有效利用多時相數(shù)據(jù)的光譜信息來識別和提取某種地物類型的動態(tài)變化信息，以及如何減小誤差累積、提高監(jiān)測精度并減少分析的復(fù)雜度，是動態(tài)變化監(jiān)測中有待解決的問題[11--13]。

筆者以秦皇島桃林口水庫為例利用多時相波段序列因子分析法進行水體動態(tài)變化信息遙感監(jiān)測等研究，把前后兩個時相的全部波段組合成一個新的波段序列作為原始變量，對其進行因子分析，從中提取出具有明確指示意義的波段組合作為因子，并利用方差分析方法對各因子進行顯著性檢驗，篩選出適合提取水體變化信息的因子，根據(jù)因子構(gòu)建水體變化指數(shù)，對水體變化信息進行定量識別和提取。

1 多時相數(shù)據(jù)組合法及因子分析原理

1.1 多時相數(shù)據(jù)組合法

多時相數(shù)據(jù)組合法是將兩個時相數(shù)據(jù)的所有波段組合為新的數(shù)據(jù)，所生成的組合數(shù)據(jù)中，各像元的反射率曲線為兩個時相波譜特性曲線的組合，稱為“地物波譜特性時相變化曲線”(圖1)。

圖1 地物波譜特性時相變化曲線示意圖Fig.1 Schematic diagram of spectral characteristic time-phase change curves of ground object

其前半部分(波段1～波段7)為同名像元前一時相(變化前)的地物波譜特性曲線，后半部分(波段8～波段14)為后一時相(變化后)地物波譜特性曲線，整個曲線代表了前后兩時相地物類型的變化信息。

單一時相遙感數(shù)據(jù)中的地物波譜特性曲線代表了某種具體的地物類型，而“地物波譜特性時相變化曲線”代表某種地物類型變化信息。通過地物波譜特性時相變化曲線，可以了解每種地物類型變化的光譜特征，對比不同的地物類型變化的光譜特征，可以發(fā)現(xiàn)每種地物類型的變化都有其各自的獨特性，根據(jù)這種獨特性建立解譯標(biāo)志，從而對此地物類型變化進行識別和提取。

1.2 因子分析

在分析地物波譜特性時相變化曲線的基礎(chǔ)上進行因子分析，確定能夠代表某種地物類型變化的因子，根據(jù)因子識別和提取此類地物類型變化信息，最后得出水體的動態(tài)變化信息。

因子分析是幫助人們對大量觀測數(shù)據(jù)進行綜合分析解釋的一種多元統(tǒng)計方法[14]。其基本思想是從協(xié)方差矩陣或相關(guān)矩陣出發(fā)，根據(jù)相關(guān)性大小把原始變量分組，使同組內(nèi)的變量間相關(guān)性較高，而不同組的變量間相關(guān)性較低。每組變量用一個不可觀測的假想變量表示，這個假想變量稱為公共因子，幾個公共因子能反映原始變量的主要信息[15]。因子分析分為R型與Q型兩種。R型因子分析用于研究變量間的關(guān)系，Q型因子分析用于研究樣本間關(guān)系[16]。本文采用R型因子分析。其矩陣形式為：

X=FAT+EC

(1)

式中：X代表標(biāo)準(zhǔn)化數(shù)據(jù)矩陣，X的(i,j)元為多時相遙感影像中像元i波段j的反射率，X(m×n)為n行m列矩陣，n為像元數(shù)(樣本數(shù))，m為波段數(shù)(變量數(shù)，全部時相的所有波段數(shù)目)。

A為公因子負(fù)載矩陣：

A=(a1,a2,…,ap)=(aij)(m×p)

(2)

式中：aij表示第i個變量在第j個公因子上的負(fù)載，即公因子j所承載的變量i的信息量，反映公因子j與變量i的關(guān)聯(lián)密切程度。

公因子得分矩陣F表示為：

F=(f1,f2,…,fp)=(fij)(m×p)

(3)

式中：fij為第j個公因子在第i個樣本上的取值，即公因子j所承載的樣本i的信息量，反映公因子j與樣本i的關(guān)聯(lián)密切程度。

另外，C和E分別表示單因子負(fù)載矩陣和單因子得分矩陣：

C=diag(c1,c2,…,cm)

(4)

E=(e1,e2,…,em)=(eij)(m×p)

(5)

式中：對角線元素ci為變量i在單因子ei上的負(fù)載(i=1,2,…,m)，即單因子ei所承載的變量i的信息量，反映單因子ei與變量i的關(guān)聯(lián)密切程度。而元素eij表示單因子ej在樣本i中的取值，即單因子ej所承載的樣本i的信息量，反映單因子ej與樣本i的關(guān)聯(lián)密切程度。

根據(jù)因子負(fù)載矩陣A的求解方式，確定最佳取值。為了合理解釋公因子數(shù)，在實際應(yīng)用中一般將多個原始變量歸結(jié)為少數(shù)公因子。其中，公因子數(shù)的選取原則如下：①信息損失無幾原則。方差貢獻代表了因子的貢獻程度，一般選定方差貢獻排序靠前的少數(shù)因子，以使得在信息損失很少的前提下解決實際問題，通常將方差貢獻累計百分比閾值設(shè)為80%，85%，90%等。②因子具有實際意義。公因子負(fù)載矩陣反映公因子與變量的關(guān)聯(lián)密切程度，可用以因子的專業(yè)意義進行理解，通常選擇能解決那些具有明確專業(yè)意義的少數(shù)公因子來對實際問題進行合理解釋。

2 研究區(qū)概況及數(shù)據(jù)預(yù)處理

研究區(qū)為秦皇島桃林口水庫，該水庫位于河北省秦皇島市西北部青龍縣境內(nèi)，主要為秦皇島、唐山兩市提供農(nóng)業(yè)生產(chǎn)和城市用水，其水位變化直接影響兩市的生活和經(jīng)濟發(fā)展。

由于近年來該流域及周邊地區(qū)受人類活動影響較為顯著，局部氣候條件發(fā)生改變，流域內(nèi)降水量明顯下降，干旱現(xiàn)象較為嚴(yán)重，造成水位下降。因此，對秦皇島桃林口水庫的水體變化信息進行遙感動態(tài)變化監(jiān)測，是該流域及水庫水資源開發(fā)、利用、優(yōu)化配置和運行管理的有力理論依據(jù)。本文的研究數(shù)據(jù)選用該地區(qū)2013年5月30日及2015年9月25日的兩時相Landsat8衛(wèi)星遙感影像，對原始數(shù)據(jù)進行了輻射定標(biāo)、FLAASH大氣校正等預(yù)處理，研究區(qū)如圖2所示。

3 水體動態(tài)變化信息提取及精度檢驗

3.1 地物類型動態(tài)變化分析

通過對兩時相的遙感影像進行解譯得出，研究區(qū)有水體、植被、裸地和工礦用地共4種地物類型，2013—2015年主要變化包括：水體變?yōu)橹脖?、水體變?yōu)樗w、水體變?yōu)槁愕?、植被變?yōu)橹脖弧⒅脖蛔優(yōu)槁愕?、裸地變?yōu)橹脖弧⒙愕刈優(yōu)橹脖缓凸さV用地變?yōu)橹脖坏裙?種地物類型變化(表1)。

表1 土地利用類型變化矩陣

將兩幅遙感影像按照時相的先后順序進行合并，選取這8種地物變化類型的訓(xùn)練樣本，并檢驗這8類訓(xùn)練樣本的J--M距離和散度，結(jié)果表明，訓(xùn)練樣本的可分離性滿足要求[17]。經(jīng)統(tǒng)計分析得出研究區(qū)內(nèi)各類訓(xùn)練樣本反射率的均值，繪制地物波譜特性時相變化曲線(圖3)，根據(jù)曲線特征分析地物波譜特性。

a. 2013；b. 2015。圖2 秦皇島桃林口水庫Landsat8衛(wèi)星遙感影像Fig.2 Landsat8 satellite remote sensing image of Taolinkou reservoir in Qinhuangdao City

圖3 地物波譜特性時相變化曲線圖Fig.3 Spectral characteristic time-phase change curves of ground object

根據(jù)該曲線可以初步了解每種地物類型變化的光譜特征，水體在前一時相(2013年)中均為水體的光譜特征，其反射率在5、6、7波段明顯低于其他地物；在后一時相(2015年)中，水體地類曲線表現(xiàn)為變化趨勢，其反射率在波段12的大小依次為：水體<裸地<植被。但也存在某種地物類型變化不顯著或與其他地物類型變化光譜特征相似的情況。從圖3中可以看出：水體變化與其他地物變化類型變化在波段1和2的光譜特征相似，因此認(rèn)為波段1和2就是無關(guān)波段，在監(jiān)測中起到干擾作用。

地物波譜特性時相變化曲線雖然能夠反映每種地物類型變化的光譜特征，但在監(jiān)測中有干擾作用的無關(guān)波段，且對各波段所含有的目標(biāo)地物動態(tài)變化信息量缺少定量評判標(biāo)準(zhǔn)，只憑借原始的地物波譜特性時相變化曲線很難準(zhǔn)確地識別和提取具體的水體地物類型變化信息。

3.2 因子分析及顯著性檢驗

為準(zhǔn)確提取變化信息，在分析地物波譜特性時相變化曲線的基礎(chǔ)上進行因子分析，通過因子分析可以對波段進行分類和歸組，定義各波段的權(quán)重系數(shù)，確定能夠代表水體地物類型變化的因子，根據(jù)水體變化因子識別和提取水體地物類型變化信息。

水體變化信息量隱含在各波段中的權(quán)重有所不同，根據(jù)各波段表達水體變化信息量的多少確定系數(shù)，構(gòu)建水體變化因子。在對原始數(shù)據(jù)進行標(biāo)準(zhǔn)化處理的基礎(chǔ)上，對研究區(qū)內(nèi)的8種地物類型變化的訓(xùn)練樣本在14個波段上的反射率進行R型因子分析，利用最大方差法提取旋轉(zhuǎn)因子解。經(jīng)過計算公因子方差，統(tǒng)計特征值、方差貢獻以及累積方差貢獻提取了方差貢獻1的前3個公因子，3個因子一共表達了91.82%的信息。為直觀分析3個因子的類別區(qū)分能力，將8種地物類型變化的訓(xùn)練樣本的反射率在3個因子上的取值兩兩組合作為X軸和Y軸制成散點圖(圖4)。

從散點圖中可以初步發(fā)現(xiàn)，因子1(F1)和因子2(F2)不能區(qū)分任何地物類型變化，而因子3(F3)可以明顯區(qū)分出水體變化信息，且可以將3種具體的水體變化類型區(qū)分開。為定量檢驗各因子對水體變化信息的動態(tài)識別能力，采用方差分析方法進行顯著性檢驗，引入可分性度量評價指標(biāo)F，判斷組間離差平方和是否顯著大于組內(nèi)離差平方和，以此衡量因子變化動態(tài)識別水平。F值計算如下：

(6)

式中：QA表示組間離差平方和；Qe表示組內(nèi)離差平方和；fA為組間自由度；fe為組內(nèi)自由度。

對水體整體變化與其他地物類型變化的整體兩組數(shù)據(jù)(記為Ⅰ)、水體變化為水體與水體變化為植被的整體兩組數(shù)據(jù)(記為Ⅱ)，對水體變化為水體與水體變化為裸地整體兩組數(shù)據(jù)(記為Ⅲ)以及水體變化為植被與水體變化為裸地的整體兩組數(shù)據(jù)(記為Ⅳ)分別進行方差分析，結(jié)果如表2所示。

表2 方差分析結(jié)果

由表2可知，在信度α=0.05下，因子F2在區(qū)分水體變?yōu)橹脖慌c水體變?yōu)槁愕貢r的FFα，所以在F1和F3上水體變化信息與其他地物變化信息差異顯著，但在3組檢驗中，F(xiàn)3的F值都遠(yuǎn)遠(yuǎn)大于F1的F值，表明F3對水體變化信息的動態(tài)識別效果最好。

3.3 水體變化因子選定及水體變化信息提取

根據(jù)因子分析及顯著性檢驗的結(jié)果，確定F3為水體變化因子。并通過計算旋轉(zhuǎn)成分矩陣得出F3的表達結(jié)果：

F3≈0.87×B5+0.71×B6+0.74×B12

(7)

式中：B5、B6、B12分別表示組合的多時相遙感數(shù)據(jù)的5、6、12波段。

基于公式(7)計算了研究區(qū)內(nèi)8種地物類型(A～H)變化樣本的因子取值(圖5)。

a. F1-F2；b. F1-F3；c. F2-F3。圖4 因子組合變化信息分析圖Fig.4 Analysis diagrams of factor combination change information

圖5 因子取值正態(tài)分布曲線Fig.5 Normal distribution curves of factor values

a. 詳圖；b. 全景圖。圖6 水體動態(tài)變化信息圖Fig.6 Water dynamic change information

由圖5可直觀看出，水體變化因子(F3)對于3種水體變化類型區(qū)分效果較好。水體地物類型變化的值均較小，與其他地物類型變化的值有較大差異，且3種具體的水體地物類型變化也能相互明顯區(qū)分，因子值大小為：水體變?yōu)樗w<水體變?yōu)槁愕?水體變?yōu)橹脖弧?/p>

基于F3水體變化因子，在計算出水體變化信息識別閾值(均值加減2倍的標(biāo)準(zhǔn)差)的基礎(chǔ)上，對研究區(qū)整幅圖像進行水體變化信息提取，包括水體變?yōu)樗w(藍色)、水體變?yōu)橹脖?紅色)、水體變?yōu)槁愕?黃色)(圖6)。

3.4 精度分析

利用選取的217個水體變?yōu)樗w樣本，162個水體變?yōu)橹脖粯颖竞?7個水體變?yōu)槁愕貥颖具M行提取結(jié)果的精度評價，基于提取結(jié)果并參考實際變化樣本，統(tǒng)計提取精度(表3)。

表3 水體變化信息提取混淆矩陣

由混淆矩陣計算得到總體精度為96.8%，Kappa系數(shù)為0.949。Kappa系數(shù)在0.81～1范圍內(nèi)表明所分結(jié)果與實際水體變化類型具有高度的一致性。由以上結(jié)果說明本文方法的水體提取精度較高，可用于水體變化信息的識別與提取。

4 結(jié)論

(1)通過多時相數(shù)據(jù)組合法將變化前后兩地物的光譜特性同時體現(xiàn)在“地物波譜特性時相變化曲線”上，可以初步了解每種地物類型變化的光譜特征。

(2)利用因子分析的理論手段，通過分析多時相遙感數(shù)據(jù)中波段序列對水體動態(tài)變化的響應(yīng)程度，得出監(jiān)測水體變化信息的3種波段組合方式(因子)，3個因子一共表達了91.82%的信息。利用方差分析方法對各因子提取水體變化信息的動態(tài)識別能力進行顯著性檢驗，得出因子F3的F值為3 373.44，遠(yuǎn)遠(yuǎn)大于其他兩個因子的F值，表明F3對水體變化信息的動態(tài)識別效果最好，可以作為水體變化識別指數(shù)。

(3)經(jīng)精度驗證分析，水體變化識別指數(shù)F3對水體變化信息的識別精度達到96.8 %，Kappa系數(shù)為0.949。Kappa系數(shù)在0.81～1范圍內(nèi),表明所分結(jié)果與實際水體變化類型具有高度的一致性。由以上結(jié)果說明本文方法的水體提取精度較高，可用于水體變化信息的識別與提取。