本發(fā)明涉及洪澇災(zāi)害監(jiān)測(cè)領(lǐng)域，特別是指一種基于時(shí)序水范圍遙感影像的洪澇災(zāi)前預(yù)警方法。

背景技術(shù)：

衛(wèi)星遙感觀測(cè)技術(shù)由于觀測(cè)范圍廣、周期性重訪等獨(dú)特優(yōu)勢(shì)，已經(jīng)應(yīng)用到包括水利、環(huán)境保護(hù)、防災(zāi)減災(zāi)等多個(gè)領(lǐng)域。應(yīng)用遙感技術(shù)進(jìn)行洪澇災(zāi)害的監(jiān)測(cè)，已經(jīng)成為當(dāng)前防災(zāi)減災(zāi)的研究重點(diǎn)。

目前，遙感技術(shù)在洪澇災(zāi)害監(jiān)測(cè)領(lǐng)域的應(yīng)用主要包括光譜遙感監(jiān)測(cè)和微波遙感監(jiān)測(cè)兩種監(jiān)測(cè)方法。對(duì)于光譜遙感監(jiān)測(cè)，光學(xué)遙感影像的光譜信息量大，易于從中有效地提取水體信息。目前應(yīng)用比較廣泛的光學(xué)數(shù)據(jù)主要有TM、ETM、MSS、SPOT、CBERS、AVHRR、MODIS、CCD、HRV、HRG、MERSI等。EOS/MODIS數(shù)據(jù)具有時(shí)間分辨率高，重訪周期短，可以獲得泛洪區(qū)域優(yōu)質(zhì)遙感影像，并為收集多時(shí)相歷史數(shù)據(jù)提供條件，同時(shí)該數(shù)據(jù)覆蓋范圍廣，適合大面積洪水監(jiān)測(cè)。MSS、TM、CBERS、SPOT等光學(xué)遙感數(shù)據(jù)中的空間分辨率比較高，一般被廣泛用于洪災(zāi)發(fā)生前土地利用信息的提取，以進(jìn)行洪澇監(jiān)測(cè)的背景數(shù)據(jù)分析，以及洪澇發(fā)生前后水體變化信息的監(jiān)測(cè)。

洪澇遙感監(jiān)測(cè)評(píng)估離不開對(duì)災(zāi)前、災(zāi)中與災(zāi)后的連續(xù)性監(jiān)測(cè)。洪澇災(zāi)害發(fā)生前主要進(jìn)行常態(tài)化的動(dòng)態(tài)監(jiān)測(cè)，收集多時(shí)相的易洪澇多發(fā)區(qū)域的水體覆蓋歷史數(shù)據(jù)，對(duì)洪峰、降水過程進(jìn)行預(yù)測(cè)預(yù)警；在洪澇災(zāi)害發(fā)生過程中，由于光譜遙感監(jiān)測(cè)技術(shù)受到觀測(cè)條件的限制，常常處于非常被動(dòng)的地位，而微波遙感監(jiān)測(cè)技術(shù)中，微波傳感器能夠有效穿透云霧，獲取洪澇災(zāi)害期間地面的水情信息，微波遙感數(shù)據(jù)主要有SAR、ASAR、SSM/I、MRI等，同時(shí)由于水體對(duì)雷達(dá)波束的鏡面反射，使得水體能夠較好地從雷達(dá)影像中提取出來，因此，在災(zāi)中主要利用微波遙感技術(shù)對(duì)洪峰經(jīng)過區(qū)域進(jìn)行監(jiān)控；洪澇災(zāi)害發(fā)生過后，主要是通過對(duì)災(zāi)前、災(zāi)中和災(zāi)后多時(shí)相遙感影像的水體覆蓋范圍的變化檢測(cè)并結(jié)合社會(huì)經(jīng)濟(jì)數(shù)據(jù)進(jìn)行災(zāi)情評(píng)估。

利用光譜傳感器對(duì)洪澇易發(fā)區(qū)域進(jìn)行長(zhǎng)期動(dòng)態(tài)監(jiān)測(cè)，并結(jié)合遙感技術(shù)建立災(zāi)前預(yù)測(cè)模型，對(duì)洪水或河道的變化進(jìn)行預(yù)測(cè)，已經(jīng)成為衛(wèi)星遙感災(zāi)前預(yù)測(cè)預(yù)警的重要組成部分。雖然不少學(xué)者在這方面仍然在進(jìn)行不懈的努力，然而目前建立災(zāi)前監(jiān)測(cè)預(yù)測(cè)模型，對(duì)洪澇災(zāi)害進(jìn)行較為準(zhǔn)確的災(zāi)前監(jiān)測(cè)仍然存在較大的困難。因此在洪澇災(zāi)前預(yù)警應(yīng)用方面還存在一定的難度，目前還沒有形成較好的方法加以解決，尚處于探索階段。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于提供一種能夠有效進(jìn)行洪澇災(zāi)前監(jiān)測(cè)預(yù)警的基于時(shí)序水范圍遙感影像的洪澇災(zāi)前預(yù)警方法。

為實(shí)現(xiàn)上述目的，本發(fā)明所提供的基于時(shí)序水范圍遙感影像的洪澇災(zāi)前預(yù)警方法，包括如下步驟：

1)輸入多時(shí)相高分辨率地表反射率遙感影像，結(jié)合歷史洪澇災(zāi)害記錄，剔除發(fā)生洪澇災(zāi)害的遙感影像，通過配準(zhǔn)和相對(duì)輻射校正對(duì)輸入的遙感影像進(jìn)行數(shù)據(jù)預(yù)處理；多時(shí)相指的是研究區(qū)6年或6年以上多年同期或者年內(nèi)某時(shí)段內(nèi)的多景影像，高分辨率指的是遙感影像的空間分辨率要優(yōu)于30m；

2)根據(jù)步驟1)中經(jīng)過預(yù)處理的遙感影像，計(jì)算水體指數(shù)，并根據(jù)水體指數(shù)提取水體范圍二值影像；所述水體指數(shù)為歸一化水體指數(shù)NDWI、改進(jìn)的歸一化差異水體指數(shù)MNDWI或新型水體指數(shù)NWI中的一種；

3)根據(jù)步驟2)中的水體范圍二值影像，提取待檢測(cè)時(shí)期水體特征和正常期水體特征；其中，待檢測(cè)時(shí)期水體特征包括待檢測(cè)時(shí)期的水體最大覆蓋范圍、水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)和水體覆蓋面積變化趨勢(shì)規(guī)律；正常期水體特征包括正常期的水體最大覆蓋范圍、水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)和水體覆蓋面積變化趨勢(shì)規(guī)律；說明：本專利中，各項(xiàng)正常期水體特征的時(shí)間段應(yīng)與待檢測(cè)時(shí)期的時(shí)間段是對(duì)應(yīng)的；各正常期水體特征在求解過程所涉及影像均不包括洪澇災(zāi)害時(shí)期的影像；

4)分別將待檢測(cè)時(shí)期的水體覆蓋范圍、水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)、水體面積變化趨勢(shì)與正常期的水體最大覆蓋范圍、水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)、水體覆蓋面積變化趨勢(shì)規(guī)律進(jìn)行疊加對(duì)比分析(一一對(duì)比)，按下述方法檢測(cè)異常信息：

4.1)將待檢測(cè)時(shí)期的水體覆蓋范圍與正常期的水體最大覆蓋范圍進(jìn)行疊置變化檢測(cè)，若某像元點(diǎn)在待檢測(cè)時(shí)期為有水，而在正常期無水，則該像元點(diǎn)作為異常點(diǎn)檢出；

4.2)將研究區(qū)待檢測(cè)時(shí)期內(nèi)的水體覆蓋時(shí)長(zhǎng)、水體覆蓋頻率分別與正常期的水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)進(jìn)行疊置統(tǒng)計(jì)分析，若待檢測(cè)時(shí)期的值背離正常期的值且超出設(shè)定閾值(不同水體特征有不同閾值)，則該像元點(diǎn)作為異常點(diǎn)檢出；

4.3)將待檢測(cè)時(shí)期內(nèi)的水體覆蓋面積變化趨勢(shì)與正常期的水體覆蓋面積變化趨勢(shì)規(guī)律進(jìn)行對(duì)比分析，若待檢測(cè)時(shí)期內(nèi)的水體覆蓋面積變化規(guī)律與正常期的變化規(guī)律相背離，則研究區(qū)的水體覆蓋面積變化趨勢(shì)異常；

5)分析跟蹤步驟4)中檢測(cè)的異常信息，若觀測(cè)區(qū)域原先本不應(yīng)該有水體覆蓋的地方被水體連續(xù)覆蓋，說明該觀測(cè)區(qū)發(fā)生洪澇災(zāi)害的風(fēng)險(xiǎn)增加，應(yīng)予以警示。

優(yōu)選地，步驟2)中，采用歸一化水體指數(shù)NDWI提取水體范圍二值影像，包括如下步驟為：

A、根據(jù)步驟1)中經(jīng)過預(yù)處理的遙感影像，利用NDWI計(jì)算公式，計(jì)算每個(gè)像素對(duì)應(yīng)的NDWI值，獲得由每個(gè)像素NDWI值組成的指數(shù)影像，所述NDWI計(jì)算公式為：

式中，p(Green)為綠波段反射率值，p(NIR)為近紅外反射率值；

B、對(duì)于生成的NDWI指數(shù)影像，選取適當(dāng)?shù)拈撝郸吝M(jìn)行二值分割，閾值α的取值區(qū)間為(-1,1)，進(jìn)一步優(yōu)選為[-0.5,0.5]，其具體合適的值還需要反復(fù)試驗(yàn)并結(jié)合實(shí)際的影像經(jīng)過目視解譯得出，以盡量準(zhǔn)確地區(qū)分出水體像元為準(zhǔn)；一般情況下，α的經(jīng)驗(yàn)值是0。指數(shù)影像中像元值大于等于閾值α的像元被提取為水體像元，其他像元為非水體像元，水體像元取值為1，非水體像元取值為0，得到的二值影像即為水體范圍二值影像。

優(yōu)選地，步驟3)中，所述正常期的水體最大覆蓋范圍NA(Normal Area)按如下方式確定：

A、按下式對(duì)研究區(qū)多年同期的水體范圍二值影像進(jìn)行疊加分析：

式中，A_k(i,j)為研究區(qū)多年同期的水體范圍二值影像，n₁為年數(shù)，i，j為對(duì)應(yīng)圖像的行列號(hào)；

B、設(shè)置適當(dāng)?shù)拈撝郸拢碌娜≈捣秶鸀?～20％，提取出正常期水體最大覆蓋范圍影像NA(i,j)；如果B(i,j)的值大于β，則NA(i,j)＝1，否則NA(i,j)＝0。待檢測(cè)時(shí)期的水體覆蓋范圍可以直接從當(dāng)前水體覆蓋范圍二值影像獲得。

本發(fā)明中，類似NA(i,j)這種描述均表示圖像NA中行列號(hào)為i,j的像素點(diǎn)處的取值，NT、NF等依此類推，不再特別說明。

優(yōu)選地，步驟3)中，所述正常期的水體覆蓋時(shí)長(zhǎng)(Normal Time)按如下方式確定：

A、提取水體覆蓋時(shí)長(zhǎng)(Time)

按照如下公式提取水體覆蓋時(shí)長(zhǎng)影像T(i,j)：

式中，A_k(i,j)表示的是研究區(qū)在年內(nèi)觀測(cè)時(shí)段內(nèi)的水體范圍二值影像，n₂為該時(shí)段影像個(gè)數(shù)，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)，m為每一景影像代表的時(shí)間尺度；

當(dāng)前待檢測(cè)時(shí)期的水體覆蓋時(shí)長(zhǎng)也根據(jù)此算法計(jì)算；

B、提取正常期水體覆蓋時(shí)長(zhǎng)(Normal Time)

按照如下公式提取正常期水體覆蓋時(shí)長(zhǎng)影像NT(i,j)：

式中，T_k(i,j)為研究區(qū)多年同期的水體覆蓋時(shí)長(zhǎng)影像，n₁為年數(shù)，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)。

優(yōu)選地，步驟3)中，所述正常期的水體覆蓋頻率按如下方式確定：

A、提取水體覆蓋頻率(Frequent)

按照如下公式提取水體覆蓋頻率影像F(i,j)：

式中，A_k(i,j)為年內(nèi)觀測(cè)時(shí)段內(nèi)的有時(shí)序(按時(shí)間序列排列)水體范圍二值影像，A₁(i,j)為像元點(diǎn)(i,j)在年內(nèi)觀測(cè)時(shí)段的開始值，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)，n₂為該觀測(cè)時(shí)段內(nèi)的影像個(gè)數(shù)；將每個(gè)像元點(diǎn)在觀測(cè)時(shí)段內(nèi)的像元值變化用一維序列X(i,j)表示，對(duì)于像元點(diǎn)(i,j)，X(i,j)＝{該像元點(diǎn)在影像1像元值，該像元點(diǎn)在影像2像元值，…，該像元點(diǎn)在影像n₂像元值}，各影像按時(shí)間順序排列；依次向后檢索一維序列X(i,j)，每次檢索2個(gè)連續(xù)元素，步長(zhǎng)為1；若檢索到的兩個(gè)元素為{0,1}，則f(0,1)＝1，否則f(0,1)＝0；

并根據(jù)此算法，計(jì)算待檢測(cè)時(shí)期的水體覆蓋時(shí)長(zhǎng)。

B、提取正常期水體覆蓋頻率(Normal Frequent)

按照如下公式提取正常期水體覆蓋頻率影像NF(i,j)：

式中，F(xiàn)_k(i,j)為多年同期水體覆蓋頻率數(shù)據(jù)影像，n₁為年數(shù)，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)。

優(yōu)選地，步驟4.1)中，所述的水體覆蓋范圍異常(Abnormal Area)，按照如下公式提取水體覆蓋范圍異常影像AbA(i,j)：

AbA(i,j)＝A(i,j)-(A(i,j)∪NA(i,j))

式中，A(i,j)，NA(i,j)分別為待檢測(cè)時(shí)期的水體范圍二值影像、正常期的水體最大覆蓋范圍影像，i，j為對(duì)應(yīng)圖像的行列號(hào)，若AbA(i,j)為1，則為異常點(diǎn)。

優(yōu)選地，步驟4.2)中，所述水體覆蓋時(shí)長(zhǎng)異常(Abnormal Time)，按照如下公式提取水體覆蓋時(shí)長(zhǎng)異常影像AbT(i,j)：

式中，T(i,j)、NT(i,j)分別為待檢測(cè)時(shí)期的水體覆蓋時(shí)長(zhǎng)影像、正常期水體覆蓋時(shí)長(zhǎng)影像。i，j為對(duì)應(yīng)圖像的行列號(hào)。若閾值的范圍在20％～100％之間，則該像元點(diǎn)作為異常檢出；取值越大，表示水體覆蓋時(shí)長(zhǎng)的偏離程度越大，異常程度越大。

優(yōu)選地，步驟4.2)中，所述水體覆蓋頻率異常(Abnormal Frequent)，按照如下公式提取水體覆蓋頻率異常影像AbF(i,j)：

式中，F(xiàn)(i,j)、NF(i,j)分別為待檢測(cè)時(shí)期的水體覆蓋頻率影像、正常期水體覆蓋頻率影像。i，j為對(duì)應(yīng)圖像的行列號(hào)。若AbF(i,j)>γ，閾值γ的范圍在20％～100％之間，則該像元點(diǎn)作為異常檢出；γ取值越大，表示水體覆蓋頻率的偏離程度越大，異常程度越大。

優(yōu)選地，步驟3)中，利用動(dòng)態(tài)諧波回歸模型DHR，擬合水體面積的正常季節(jié)周期規(guī)律，得到正常期的水體覆蓋面積變化趨勢(shì)規(guī)律；

步驟4.3)中，所述的水體覆蓋面積動(dòng)態(tài)趨勢(shì)異常，按如下公式判斷是否為異常信息：

式中，表示t_i時(shí)間節(jié)點(diǎn)真值與擬合值的差值，表示待檢測(cè)的t_i時(shí)間節(jié)點(diǎn)的水體覆蓋面積值，表示采用動(dòng)態(tài)諧波回歸模型DHR模擬得到的t_i時(shí)間節(jié)點(diǎn)的正常期水體覆蓋面積值，t_i時(shí)間節(jié)點(diǎn)表示在觀測(cè)時(shí)段內(nèi)的任一一個(gè)時(shí)間點(diǎn)，若則為異常。

優(yōu)選地，該方法還包括如下步驟：采用動(dòng)態(tài)諧波回歸模型(DHR)，基于步驟2)中提取的水體范圍二值影像進(jìn)行時(shí)空分析，以獲得研究區(qū)水體覆蓋面積的時(shí)空變化規(guī)律，包括水體面積的長(zhǎng)期變化趨勢(shì)，水體面積變化的季節(jié)周期規(guī)律以及擬合出的水體面積多年的正常變化趨勢(shì)規(guī)律。

本發(fā)明的有益效果是：本發(fā)明通過光譜傳感器對(duì)洪澇易發(fā)區(qū)域進(jìn)行長(zhǎng)時(shí)間動(dòng)態(tài)連續(xù)監(jiān)測(cè)，并獲取連續(xù)的遙感影像數(shù)據(jù)，建立相關(guān)模型，獲取水體覆蓋范圍、水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)以及水體覆蓋面積的變化規(guī)律這幾個(gè)水體特征的正常期變化規(guī)律，從多個(gè)角度來分析水體特征的異常，根據(jù)異常分析對(duì)觀測(cè)區(qū)發(fā)生洪澇災(zāi)害的風(fēng)險(xiǎn)進(jìn)行預(yù)警，從而有效解決洪澇災(zāi)害災(zāi)前監(jiān)測(cè)預(yù)警的困難問題。

附圖說明

圖1為實(shí)施例中所提供的洪澇災(zāi)前預(yù)警方法的總體流程圖。

圖2為實(shí)施例中多時(shí)相多角度水體異常檢測(cè)技術(shù)流程圖。

圖3為實(shí)施例中水體覆蓋范圍異常檢測(cè)的結(jié)果圖。

圖4為實(shí)施例中水體覆蓋時(shí)長(zhǎng)異常檢測(cè)的結(jié)果圖。

圖5為實(shí)施例中水體覆蓋頻率異常檢測(cè)的結(jié)果圖。

圖6為實(shí)施例中水體覆蓋面積趨勢(shì)異常檢測(cè)及時(shí)空分析結(jié)果圖。

具體實(shí)施方式

下面結(jié)合附圖和具體實(shí)施例對(duì)本發(fā)明作進(jìn)一步的詳細(xì)說明。

如圖1～2所示，本發(fā)明所提供的基于時(shí)序水范圍遙感影像的洪澇災(zāi)前預(yù)警方法，包括如下步驟：

1)輸入多時(shí)相高分辨率地表反射率遙感影像，結(jié)合歷史災(zāi)前記錄，剔除發(fā)生洪澇災(zāi)害的數(shù)據(jù)，通過配準(zhǔn)和相對(duì)輻射校正對(duì)輸入的遙感影像進(jìn)行數(shù)據(jù)預(yù)處理。

2)根據(jù)步驟1)獲得的地表反射率影像計(jì)算水體指數(shù)，采用歸一化水體指數(shù)NDWI提取水體范圍，包括如下步驟為：

A、提取水體范圍的方法有很多，現(xiàn)有的遙感影像水體提取方法主要分為以下四類：特征分割方法、模式分類方法、輪廓線檢測(cè)方法以及多源數(shù)據(jù)綜合分析方法。特征指數(shù)法是通過多光譜波段運(yùn)算獲取反映水體與背景反差的指數(shù)，進(jìn)而利用直方圖閾值分割得到提取結(jié)果，典型的水體特征指數(shù)如歸一化水體指數(shù)(NDWI)、改進(jìn)的歸一化差異水體指數(shù)(MNDWI)、新型水體指數(shù)(NWI)等，其模型簡(jiǎn)單，并且有不錯(cuò)的提取結(jié)果。在這里我們使用NDWI閾值分割法來提取水體范圍，根據(jù)輸入的多時(shí)相的地表反射率遙感影像，利用NDWI計(jì)算公式，計(jì)算每個(gè)像素對(duì)應(yīng)的NDWI值。

B、NDWI的計(jì)算表達(dá)式為：

其中，p(Green)為綠波段反射率值，p(NIR)為近紅外反射率值。

對(duì)于生成的NDWI指數(shù)影像，選取適當(dāng)?shù)拈撝郸吝M(jìn)行二值分割。選取適當(dāng)?shù)拈撝郸吝M(jìn)行二值分割，閾值α的取值區(qū)間為(-1,1)，進(jìn)一步優(yōu)選為[-0.5,0.5]，其具體合適的值還需要反復(fù)試驗(yàn)并結(jié)合實(shí)際的影像經(jīng)過目視解譯得出，以盡量準(zhǔn)確地區(qū)分出水體像元為準(zhǔn)；一般情況下，α的經(jīng)驗(yàn)值是0。指數(shù)影像中像元值大于等于閾值α的像元被提取為水體像元，其他像元為非水體像元，水體像元取值為1，非水體像元取值為0，得到的二值影像即為水體覆蓋范圍提取結(jié)果。

實(shí)施例中，根據(jù)鄱陽湖流域2009-2014年多年4-12月HJ1A衛(wèi)星數(shù)據(jù)影像，4-12月每月一景，基于NDWI計(jì)算公式，計(jì)算這些影像的NDWI值；并根據(jù)NDWI影像，提取了研究區(qū)的水體范圍。

3)根據(jù)步驟2)中的水體范圍數(shù)據(jù)，按下述方法提取待檢測(cè)時(shí)期水體特征和正常期水體特征：

3.1)根據(jù)步驟2)中的水體范圍數(shù)據(jù)，提取水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)，并獲取待檢測(cè)時(shí)期的水體特征，包括待檢測(cè)時(shí)期水體最大覆蓋范圍、待檢測(cè)時(shí)期水體覆蓋頻率、待檢測(cè)時(shí)期水體覆蓋時(shí)長(zhǎng)和待檢測(cè)時(shí)期水體覆蓋面積變化趨勢(shì)規(guī)律；

3.2)根據(jù)步驟2)中的水體范圍和步驟3.1)中的水體覆蓋頻率和水體覆蓋時(shí)長(zhǎng)，提取正常期水體特征，包括正常期水體最大覆蓋范圍、正常期水體覆蓋頻率、正常期水體覆蓋時(shí)長(zhǎng)和正常期水體覆蓋面積變化趨勢(shì)規(guī)律；

其中，正常期水體最大覆蓋范圍(Normal Area)是指根據(jù)研究區(qū)多年同期的水體覆蓋范圍影像，通過疊加分析，提取的在該時(shí)期研究區(qū)被水體覆蓋的最大正常范圍。

按如下方式確定正常期水體最大覆蓋范圍影像NA(i,j)：

A、按下式對(duì)研究區(qū)多年同期(多年的含義是6年或6年以上，同期的含義是指年內(nèi)某日的一景影像或某旬或月內(nèi)的一景多天合成影像)的水體范圍影像進(jìn)行疊加分析：

式中，A_k(i,j)為研究區(qū)多年同期的水體覆蓋范圍二值影像，n₁為年數(shù)，i，j為對(duì)應(yīng)圖像的行列號(hào)，計(jì)算時(shí)按照影像的時(shí)間順序依次相加。

B、設(shè)置適當(dāng)?shù)拈撝郸?，β的取值范圍?～20％，提取出研究區(qū)被水體覆蓋的最大正常范圍NA(i,j)，如果B(i,j)的值大于β，則NA(i,j)＝1，否則，NA(i,j)＝0。

閾值β的選取是為了剔除掉某像元點(diǎn)由于云覆蓋或偶然降水導(dǎo)致的水體像元。

實(shí)施例中，基于HJ1A衛(wèi)星數(shù)據(jù)獲得的鄱陽湖流域2009-2014年多年6月的水體覆蓋范圍數(shù)據(jù)影像，根據(jù)歷史洪澇災(zāi)害記錄，2010年6月鄱陽湖流域發(fā)生了洪澇災(zāi)害，所以剔除了2010年的數(shù)據(jù)。利用步驟3.2)中正常期水體最大覆蓋范圍的算法，將其他年份的數(shù)據(jù)進(jìn)行疊加分析，閾值a取為10，提取鄱陽湖6月的正常期最大水體覆蓋范圍。如圖3所示，左圖表示的是鄱陽湖在6月份未發(fā)生洪澇災(zāi)害情況時(shí)的最大覆蓋面積；中間表示的是鄱陽湖2015年6月的水體覆蓋范圍。

步驟3.1)中，水體覆蓋時(shí)長(zhǎng)(Time)是指研究區(qū)在年內(nèi)某觀測(cè)時(shí)段內(nèi)被水體覆蓋的總時(shí)長(zhǎng)。

按如下公式得到水體覆蓋時(shí)長(zhǎng)影像T(i,j)：

式中，A_k(i,j)表示的是研究區(qū)在年內(nèi)該觀測(cè)時(shí)段內(nèi)的時(shí)間序列水體覆蓋范圍二值影像，n₂為該時(shí)段影像個(gè)數(shù)，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)，m為每一景影像代表的時(shí)間尺度，計(jì)算時(shí)按照影像的時(shí)間順序依次相加。

步驟3.2)中，正常期的水體覆蓋時(shí)長(zhǎng)(Normal Time)是指根據(jù)研究區(qū)多年同期的水體覆蓋時(shí)長(zhǎng)影像，通過疊加分析，計(jì)算研究區(qū)各像元點(diǎn)在該觀測(cè)時(shí)段內(nèi)被水覆蓋的平均時(shí)長(zhǎng)。

按如下公式得到正常期水體覆蓋時(shí)長(zhǎng)影像NT(i,j)：

式中，T_k(i,j)為研究區(qū)多年同期水體覆蓋時(shí)長(zhǎng)影像，n₁為年數(shù)，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)。

實(shí)施例中，基于HJ1A衛(wèi)星數(shù)據(jù)所獲得的鄱陽湖流域2009-2014年4-12月的水體覆蓋范圍影像，4-12觀測(cè)時(shí)間段內(nèi)一共9景影像，每月一景，每一景影像代表的時(shí)間尺度為30天。根據(jù)步驟3.1)中的計(jì)算水體覆蓋時(shí)長(zhǎng)的算法，分別計(jì)算2009-2014年每年4-12月這一觀測(cè)時(shí)段的水體覆蓋時(shí)長(zhǎng)；再利用步驟3.2)中的正常期水體覆蓋時(shí)長(zhǎng)算法，根據(jù)2009-2014年6年的水體覆蓋時(shí)長(zhǎng)影像計(jì)算鄱陽湖4-12月的正常期水體覆蓋時(shí)長(zhǎng)。如圖4所示，圖4中的左圖表示的是在4-12月這個(gè)時(shí)間段，鄱陽湖流域在未發(fā)生洪澇災(zāi)害時(shí)，被水體覆蓋的正常時(shí)長(zhǎng)，為水體覆蓋時(shí)長(zhǎng)異常檢測(cè)提供正常期的參照基礎(chǔ)；圖4中的中間圖表示的是2015年4-12月的水體覆蓋時(shí)長(zhǎng)，為待檢測(cè)時(shí)期的水體覆蓋時(shí)長(zhǎng)。

步驟3.1)中，水體覆蓋頻率(Frequent)是指研究區(qū)在年內(nèi)觀測(cè)時(shí)段內(nèi)被水體覆蓋變動(dòng)的次數(shù)或頻次。

按如下公式得到年內(nèi)該觀測(cè)時(shí)段的水體覆蓋頻率影像F(i,j)：

式中，年內(nèi)某觀測(cè)時(shí)段內(nèi)有時(shí)序水體范圍二值影像A_k(i,j)，A₁(i,j)表示像元點(diǎn)(i,j)在年內(nèi)觀測(cè)時(shí)段的開始值，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)，n₂為該觀測(cè)時(shí)段內(nèi)的影像個(gè)數(shù)；將每個(gè)像元點(diǎn)在觀測(cè)時(shí)段內(nèi)的像元值變化用一維序列X(i,j)表示，對(duì)于像元點(diǎn)(i,j)，X(i,j)＝{該像元點(diǎn)在影像1像元值，該像元點(diǎn)在影像2像元值，…，該像元點(diǎn)在影像n₂像元值}；依次向后檢索一維序列X(i,j)，每次檢索2個(gè)連續(xù)元素，步長(zhǎng)為1；若檢索到的兩個(gè)元素為{0,1}，則f(0,1)＝1，否則f(0,1)＝0。

步驟3.2)中所述的正常期的水體覆蓋頻率NF(Normal Frequent)是指根據(jù)研究區(qū)多年同期的水體覆蓋頻率影像，通過疊加分析，計(jì)算研究區(qū)各像元點(diǎn)在該觀測(cè)時(shí)段內(nèi)被水覆蓋的平均頻率。

按如下公式得到正常期水體覆蓋頻率影像NF(i,j)：

式中，F(xiàn)_k(i,j)為研究區(qū)多年同期水體覆蓋頻率數(shù)據(jù)影像，n₁為年數(shù)，i，j為像元點(diǎn)對(duì)應(yīng)的行列號(hào)。

實(shí)施例中，基于HJ1A衛(wèi)星數(shù)據(jù)所獲得的鄱陽湖流域2009-2014年4-12月的水體覆蓋范圍影像，4-12觀測(cè)時(shí)間段內(nèi)一共9景影像，每月一景。根據(jù)步驟3.1)中計(jì)算水體覆蓋頻率的算法，分別計(jì)算2009-2014年每年年內(nèi)的水體覆蓋頻率；再利用步驟3.2)中正常期水體覆蓋頻率的算法，根據(jù)2009-2014多年的水體覆蓋頻率影像計(jì)算正常期水體覆蓋頻率。如圖5所示，圖5中的左圖表示的是鄱陽湖流域在4-12月這個(gè)時(shí)間段，未發(fā)生洪澇災(zāi)害時(shí)被水體覆蓋的正常頻率，為水體覆蓋頻率異常檢測(cè)提供正常期的參照基礎(chǔ)；圖5中的中間圖表示的是2015年4-12月這一觀測(cè)時(shí)段內(nèi)的水體覆蓋頻率，為待檢測(cè)水體覆蓋頻率。

步驟3.2)中所述的正常期水體覆蓋面積變化趨勢(shì)規(guī)律是指研究區(qū)的水體覆蓋面積在觀測(cè)時(shí)段內(nèi)的正常變化趨勢(shì)規(guī)律。

按如下方式得到正常期水體覆蓋面積變化趨勢(shì)規(guī)律：

根據(jù)研究區(qū)長(zhǎng)時(shí)間序列的水體覆蓋范圍二值影像，基于動(dòng)態(tài)諧波回歸模型(DHR)，擬合水體面積的正常變化趨勢(shì)規(guī)律。

動(dòng)態(tài)諧波回歸模型(DHR)基于UC模型發(fā)展而來，UC模型有多種表現(xiàn)形式，DHR(Dynamic Harmonic Regession)模型可以認(rèn)為是UC模型的一種特殊表現(xiàn)形式(Young等，1999)，它建立在UC模型的基礎(chǔ)上。

對(duì)于一個(gè)時(shí)間序列，DHR模型主要用于擬合時(shí)間序列的3個(gè)成分，即T_t、S_t及e_t，表達(dá)為：

y_t＝T_t+S_t+e_t e_t～N(0,σ²)

式中，y_t表示原始時(shí)間序列在t時(shí)刻的值；T_t表示為原始序列在該時(shí)刻的趨勢(shì)成分；S_t表示原始序列在該時(shí)刻的周期成分，其中S_t為季節(jié)周期成分；e_t為殘差，通常為了簡(jiǎn)便，e_t被視作為高斯白噪聲分布。本技術(shù)中的研究對(duì)象為水體面積，水體面積具有季節(jié)和周期變化規(guī)律，DHR模型的主要特點(diǎn)是對(duì)季節(jié)和周期成分進(jìn)行擬合，因此它適用于有季節(jié)或周期變化成分的時(shí)間序列分析。

本技術(shù)中的DHR模型是以MATLAB平臺(tái)的Captain Toolbox軟件包為基礎(chǔ)進(jìn)行開發(fā)建模的，該軟件包由Lancaster University的Young等設(shè)計(jì)并發(fā)展(http://www.es.lancs.ac.uk/cres/captain/)。根據(jù)開發(fā)的DHR模型功能，直接輸入長(zhǎng)時(shí)間序列水體覆蓋范圍影像，就可以擬合出研究區(qū)水體面積的正常變化趨勢(shì)規(guī)律。

本實(shí)施例中，基于HJ1A衛(wèi)星數(shù)據(jù)所獲得的鄱陽湖流域2009-2014年4-12月的水體覆蓋范圍影像，4-12觀測(cè)時(shí)間段內(nèi)一共9景影像，每月一景。根據(jù)步驟3.2)中計(jì)算正常期水體覆蓋面積變化趨勢(shì)的算法，利用動(dòng)態(tài)諧波回歸模型(Dynamic Harmonic Regession)對(duì)研究區(qū)的水體覆蓋面積進(jìn)行分析。如圖6所示有4條曲線，在本步驟中，需要的結(jié)果是擬合水體面積，表示的是研究區(qū)的真實(shí)水體面積在去除噪聲的擾動(dòng)后，擬合的水體覆蓋面積的正常變化趨勢(shì)規(guī)律，為水體面積趨勢(shì)異常檢測(cè)提供正常依據(jù)。

4)分別將待檢測(cè)時(shí)期的水體覆蓋范圍、水體覆蓋頻率、水體覆蓋時(shí)長(zhǎng)和水體覆蓋面積變化趨勢(shì)與正常期水體最大覆蓋范圍、正常期水體覆蓋頻率、正常期水體覆蓋時(shí)長(zhǎng)以及正常期水體覆蓋面積變化趨勢(shì)規(guī)律進(jìn)行疊加對(duì)比分析，檢測(cè)異常信息。

4.1)水體覆蓋范圍異常檢測(cè)是以傳統(tǒng)異常檢測(cè)方法理論為基礎(chǔ)的異常檢測(cè)方法，即在傳統(tǒng)異常檢測(cè)理論中，認(rèn)為如數(shù)據(jù)與給定統(tǒng)計(jì)分布或模型的顯著性差異超過某一特定數(shù)值或范圍即為異常。在進(jìn)行水體覆蓋范圍異常檢測(cè)時(shí)，將當(dāng)前待檢測(cè)的水體范圍影像與該時(shí)間的正常期最大水體覆蓋范圍進(jìn)行疊置變化檢測(cè)分析，提取出與正常期的水體分布不一致部分，即為異常的水體覆蓋范圍。

4.2)水體覆蓋時(shí)長(zhǎng)或頻率異常檢測(cè)是以基于傳統(tǒng)異常檢測(cè)方法理論為基礎(chǔ)的異常檢測(cè)方法，即在傳統(tǒng)異常檢測(cè)理論中，認(rèn)為如數(shù)據(jù)與給定統(tǒng)計(jì)分布或模型的顯著性差異超過某一特定數(shù)值或范圍即為異常。在進(jìn)行水體覆蓋時(shí)長(zhǎng)或頻率異常檢測(cè)時(shí)，將當(dāng)前待檢測(cè)的水體覆蓋時(shí)長(zhǎng)或頻率影像與同期的正常期特征進(jìn)行疊置閾值分析，若待檢測(cè)時(shí)期的值超過正常期特征值一定閾值之外則認(rèn)為是異常的。

4.3)水體覆蓋面積動(dòng)態(tài)趨勢(shì)異常檢測(cè)是以偏差異常檢測(cè)理論中基于輪廓的方法為基礎(chǔ)的一種異常檢測(cè)方法?；谳喞漠惓z測(cè)方法，首先給定一個(gè)時(shí)間序列的正常輪廓，將新的時(shí)間點(diǎn)或時(shí)間段與輪廓進(jìn)行比對(duì)來識(shí)別其異常與否。這里，輪廓實(shí)際上是對(duì)時(shí)間序列的歷史數(shù)據(jù)進(jìn)行平滑而得到的一個(gè)標(biāo)準(zhǔn)序列，即這類的方法的核心在于如何獲得時(shí)間序列的正常輪廓。

步驟4.1)中所述的水體覆蓋范圍異常(Abnormal Area)，按照如下公式提取異常水體范圍影像AbA(i,j)：

AbA(i,j)＝A(i,j)-(A(i,j)∪NA(i,j))

式中，A(i,j)，NA(i,j)分別為待檢測(cè)時(shí)期的水體范圍二值影像、正常期最大水體范圍二值影像，i，j為對(duì)應(yīng)圖像的行列號(hào)，若AbA(i,j)為1，則為異常點(diǎn)。

實(shí)施例中，根據(jù)步驟3.2)獲取的6月鄱陽湖流域正常期最大水體覆蓋范圍并結(jié)合2015年6月的水體范圍影像，利用步驟4.1)中檢測(cè)水體范圍異常的算法，提取的異常水體范圍如圖3中右圖所示。通過對(duì)研究區(qū)水體范圍的異常檢測(cè)，我們就能夠知道待檢測(cè)時(shí)期的水體范圍是否超出了正常期的最大范圍，若超出了，就應(yīng)當(dāng)及時(shí)預(yù)警。

步驟4.2)中所述的水體覆蓋時(shí)長(zhǎng)異常(Abnormal Time)，按照如下公式提取水體覆蓋時(shí)長(zhǎng)異常影像AbT(i,j)：

式中，T(i,j)、NT(i,j)分別為待檢測(cè)水體覆蓋時(shí)長(zhǎng)影像、正常期水體覆蓋時(shí)長(zhǎng)影像。i，j為對(duì)應(yīng)圖像的行列號(hào)。若閾值的范圍在20％-100％之間，則該像元點(diǎn)作為異常檢出；取值越大，表示水體覆蓋時(shí)長(zhǎng)的偏離程度越大，異常程度越大。

實(shí)施例中，根據(jù)步驟3.2)獲取的鄱陽湖流域4-12月的正常期水體覆蓋時(shí)長(zhǎng)并結(jié)合2015年4-12月的水體覆蓋時(shí)長(zhǎng)影像，利用步驟4.2)中檢測(cè)水體覆蓋時(shí)長(zhǎng)異常的算法，閾值選擇為20％，檢測(cè)水體覆蓋時(shí)長(zhǎng)異常，并輸出異常點(diǎn)的空間分布，結(jié)果如圖4中右圖所示。水體覆蓋時(shí)長(zhǎng)異常表示的是異常區(qū)域在研究時(shí)段內(nèi)被水體淹沒的時(shí)長(zhǎng)超出了正常期的值，表明異常區(qū)域原先是沒有水或者是較少被水覆蓋的，應(yīng)當(dāng)及時(shí)預(yù)警。

步驟4.2)所述的水體覆蓋頻率異常(Abnormal Frequent)，按照如下公式提取水體覆蓋頻率異常影像AbF(i,j)：

式中，F(xiàn)(i,j)、NF(i,j)分別為待檢測(cè)水體覆蓋頻率影像、正常期水體覆蓋頻率影像。i，j為對(duì)應(yīng)圖像的行列號(hào)。若AbF(i,j)>γ，閾值γ的范圍在20％-100％之間，則該像元點(diǎn)作為異常檢出；γ取值越大，表示水體覆蓋頻率的偏離程度越大，異常程度越大。

實(shí)施例中，根據(jù)步驟3.2)獲取的鄱陽湖流域4-12月的正常期水體覆蓋頻率并結(jié)合2015年4-12月的水體覆蓋頻率影像，利用步驟4.2中檢測(cè)水體覆蓋頻率異常的算法，閾值γ選擇為20％，檢測(cè)水體覆蓋頻率異常，并輸出異常點(diǎn)的空間分布，結(jié)果分別如圖5中右圖所示。水體覆蓋頻率異常表示的是異常區(qū)域在研究時(shí)段內(nèi)被水體頻繁淹沒，而異常表明這些異常區(qū)域的被水體淹沒的風(fēng)險(xiǎn)提高，應(yīng)當(dāng)及時(shí)預(yù)警。

步驟4.3)中所述的水體覆蓋面積動(dòng)態(tài)趨勢(shì)異常檢測(cè)，按照如下方式提取異常信息：

式中，表示t_i時(shí)間節(jié)點(diǎn)真值與擬合值的差值，表示t_i時(shí)間節(jié)點(diǎn)的真實(shí)的水體面積值，表示t_i時(shí)間節(jié)點(diǎn)的擬合的正常期水體面積值，如果則為異常，否則為非異常。

實(shí)施例中，根據(jù)步驟3.2)的算法，所獲取的鄱陽湖水體面積正常變化趨勢(shì)規(guī)律，利用步驟4.3)中檢測(cè)水體覆蓋面積動(dòng)態(tài)趨勢(shì)異常的算法，將水體面積的真實(shí)值與擬合的正常趨勢(shì)規(guī)律疊加分析，如圖6所示，圓圈部分為真值高于擬合的正常規(guī)律的時(shí)間段，若異常持續(xù)的時(shí)間過長(zhǎng)，則說明該區(qū)域被水體淹沒的面積持續(xù)增大，發(fā)生洪澇的風(fēng)險(xiǎn)增加，應(yīng)該及時(shí)預(yù)警。

步驟5，基于DHR模型，對(duì)研究區(qū)多年長(zhǎng)時(shí)間序列的水體覆蓋范圍影像進(jìn)行時(shí)空分析，探究研究區(qū)水體覆蓋面積的時(shí)空變化規(guī)律。能夠分析出研究區(qū)水體面積的多年長(zhǎng)期變化趨勢(shì)，水體面積年內(nèi)的季節(jié)周期變化規(guī)律。

實(shí)施例中，如圖6，趨勢(shì)，即研究區(qū)在這段時(shí)間內(nèi)水體面積總體變化的趨勢(shì)；周期，即研究區(qū)的水體覆蓋面積在年內(nèi)的季節(jié)周期變化規(guī)律。從水體面積變化的季節(jié)周期規(guī)律可以看出，湖水面積會(huì)隨著季節(jié)的變化呈現(xiàn)出豐水期和枯水期的交替變化，這與鄱陽湖是一個(gè)季節(jié)性吞吐湖是相符合的；另外從水體面積變化趨勢(shì)線中可以看出，研究區(qū)在這段時(shí)間水體面積的整體變化趨勢(shì)呈現(xiàn)輕微下降趨勢(shì)，和近年來鄱陽湖水體面積逐漸萎縮的實(shí)際情況相符合。