陳飛宇，胡友彬，施恩

（國防科技大學(xué)氣象海洋學(xué)院，南京 211101）

0　引言

海岸線是陸地與海洋相互交匯的部位，包括大陸海岸線和島嶼海岸線，不僅標(biāo)識了沿海地區(qū)的水陸分界線，而且蘊(yùn)含著豐富的環(huán)境信息，并對沿海的灘涂面積、濕地生態(tài)系統(tǒng)的興衰等具有重要的指示作用。近些年來，遙感技術(shù)以其觀測范圍廣、信息量大、獲取信息快、更新周期短等優(yōu)點，逐漸成為海岸線提取的重要方法和手段。

混合像元是指在傳感器的瞬時視場中包含多種地物類別的像元（如圖1所示），是遙感影像所固有的特征[1-2]。硬分類方法將混合像元完全歸為某一種地物類型，會導(dǎo)致分類結(jié)果不準(zhǔn)確。如何對遙感影像中這種廣泛存在的混合像元問題進(jìn)行處理已經(jīng)成為定量遙感分析中最為關(guān)鍵的問題[3-4]。亞像元定位（Sub-Pixel Mapping，SPM），也稱超分辨率制圖（Super-Resolution Mapping，SRM），是目前解決混合像元問題的常用手段，其目的是確定混合像元中不同地物類型的具體空間位置，獲取更高空間分辨率的地物分類圖。因此，亞像元定位可以看作是一種提高遙感影像空間分辨率的技術(shù)，即將低分辨率遙感影像的混合像元分解轉(zhuǎn)換成高分辨率的地物分類圖。海岸線的亞像元定位是指在混合像元內(nèi)部確定水、陸的空間分布，國內(nèi)外學(xué)者對此做出了深入地研究，取得了豐碩的成果。張旸等[5]使用黃河三角洲海岸Landsat衛(wèi)星遙感數(shù)據(jù)，基于研究區(qū)域低分辨率6波段的海陸類型軟分類結(jié)果及其變差函數(shù)，以高分辨率8波段的指示變差函數(shù)為精細(xì)尺度先驗信息模型，采用地統(tǒng)計學(xué)方法生成海陸類型發(fā)生概率模擬圖像，通過等值線法提取海岸線空間分布特征。凌峰等[6]考慮到亞像元定位是一個欠定反演問題，其約束條件遠(yuǎn)遠(yuǎn)少于求解參數(shù)，采用高分辨率的數(shù)字高程模型作為海岸線亞像元定位的輔助信息，有效提高了海岸線亞像元定位的精度。Foody等[7]詳細(xì)分析了對海岸線提取進(jìn)行亞像元定位分析的必要性，分別通過擬合軟分類的類別成員輪廓法和地統(tǒng)計學(xué)法進(jìn)行海岸線的亞像元定位，并與硬分類方法的提取結(jié)果進(jìn)行了對比。Muslim等[8]認(rèn)為在同類地物光譜差異較大的遙感影像中選取全局的陸地端元和水體端元，缺乏代表性，誤差較大，提出了基于局部軟分類的海岸線亞像元定位技術(shù)，該方法根據(jù)光譜差異，選取若干個樣本區(qū)域，再將提取的海岸線分成若干段，每段海岸線就近選取地類端元，并將該思想分別應(yīng)用到兩點直方圖法和像元交換法中，取得了較好的效果。Muad等[9]以多幅時間序列的遙感影像為數(shù)據(jù)源，采用Hopfield神經(jīng)網(wǎng)絡(luò)方法對湖泊邊界進(jìn)行亞像元定位，與單幅影像相比，精度顯著提高。

然而，海岸線的亞像元定位結(jié)果受地類組分誤差的影響較大，采用全局端元進(jìn)行光譜分解得到地類組分的亞像元定位方法并未考慮遙感影像中屬于同類地物的不同像元可能存在的較大光譜差異的情況。例如，在現(xiàn)實生活中，陸地可能是淤泥、巖石或者植被等多種地物，反映在遙感影像上就是光譜值的多種多樣，同樣，海水隨著溫度、鹽度等的變化，其反射率的值也各不相同。因此，選取全局的陸地和海水的端元，不能充分反映屬于同類地物的不同像元之間存在的較大的光譜差異，從而造成混合像元分解誤差，影響海岸線亞像元定位的精度。

針對上述問題，本文提出了基于GAC的局部自適應(yīng)海岸線亞像元定位（Locally Adaptive Sub-Pixel Shoreline Mapping based on GAC algorithm,LA_SPSMG）算法。該算法首先采用GAC模型提取初始海岸線，然后通過數(shù)學(xué)形態(tài)學(xué)腐蝕和膨脹操作，得到需要進(jìn)行亞像元定位的混合像元，最后采用局部自適應(yīng)選取地物端元、融合光譜信息的亞像元定位方法提取海岸線，減緩了采用全局端元引起的光譜分解誤差的影響，提高了海岸線亞像元定位的精度。

圖1　純像元與混合像元示意圖

1　基于GAC的初始海岸線提取和混合像元確定

對海岸線進(jìn)行亞像元定位，首先要確定需要進(jìn)行亞像元定位的混合像元集合。根據(jù)水、陸的空間位置關(guān)系可知，海陸交界處的像元最有可能為混合像元。因此，本文利用GAC模型得到初始的海岸線，然后通過形態(tài)學(xué)膨脹和腐蝕等操作可得到以初始海岸線為中心的向水、陸兩側(cè)擴(kuò)展一定寬度的帶狀區(qū)域，假設(shè)該帶狀區(qū)域內(nèi)的像元為混合像元，這個帶狀區(qū)域之外的像元為純水像元或者純陸像元。

1.1　初始海岸線提取

測地線活動輪廓（Geodesic Active Contour,GAC）模型[10]是一種基于變分原理的圖像分割方法，其主要思想是利用黎曼空間中的測地線概念，把尋找圖像中的邊界線問題轉(zhuǎn)化為尋找一條加權(quán)弧長最小值問題，利用能量泛函，結(jié)合圖像低層數(shù)據(jù)和高層信息（如目標(biāo)邊緣輪廓的連續(xù)性、閉合性、形狀先驗信息等）來完成分割。該模型具有對噪聲不敏感、檢測精度高等優(yōu)點，所得到的光滑連續(xù)曲線可彌補(bǔ)目標(biāo)輪廓上的噪聲、間隙以及其它不規(guī)則形狀，非常適合于海岸線這類具有復(fù)雜特征地物的輪廓提取。因此，本文以遙感影像為數(shù)據(jù)源，采用GAC模型得到初始的海岸線。

設(shè)γ(q):[0,1]→R2為平面參數(shù)曲線，I:[0,a]×[0,b]→R+為待檢測圖像，定義能量泛函

其中，g(·)為停止函數(shù)：

其中，Gσ是方差為σ的高斯核函數(shù)，因此，g為關(guān)于圖像梯度強(qiáng)度的遞減函數(shù)。

計算能量泛函式的變分，曲線演化的梯度下降流為：

其中，κ是歐氏空間曲率，N→是單位內(nèi)法向量。

利用差分法對式（3）進(jìn)行離散化求解，經(jīng)過若干次的迭代演化，可得到海岸線輪廓。

1.2　混合像元的確定

數(shù)學(xué)形態(tài)學(xué)的基本原理是將結(jié)構(gòu)元素在圖像范圍內(nèi)平移，同時施加交、并等基本集合運(yùn)算，從而保留主要形狀，刪除不相干形狀（如噪聲、毛刺等）。設(shè)初始海岸線的二值化圖像為A，模板（結(jié)構(gòu)元素）為B，則A被B膨脹表示為A⊕B，得到的膨脹后的圖像為：

A被B腐蝕表示為 AΘB，得到的腐蝕后的圖像為：

則需要進(jìn)行亞像元定位的混合像元集合表示為：

通過純陸像元、混合像元和純水像元三者之間的空間位置關(guān)系，可進(jìn)一步得到純水像元集合pw和純陸像元集合pl。

如圖2所示，圖2（a）中的研究區(qū)域是一個小島，提取初始海岸線后，將得到的二值化圖像進(jìn)行膨脹和腐蝕，中間的帶狀區(qū)域（圖2（b））即為需要進(jìn)行亞像元定位的混合像元的集合，再根據(jù)海陸之間的空間關(guān)系，易知混合像元帶狀區(qū)域往里的部分是純陸像元的集合，往外的部分是純水像元的集合，這樣就得到了混合像元集合mp、純水像元集合pw和純陸像元集合pl。

圖2　混合像元集合、純水像元集合和純陸像元集合的確定

2　基于全局端元的亞像元定位方法

假設(shè)進(jìn)行亞像元定位的低分辨率遙感影像Y共有L(L≥1)個波段，Y的每個波段圖像共包括M×N個像元，放大因子為s，Y的每個像元被分割為s2個亞像元。輸入數(shù)據(jù)為低分辨率遙感影像Y，輸出數(shù)據(jù)為亞像元定位圖像X，X包含(M×s)×(N×s)個亞像元。li(j)∈{1,…,K }代表第i個混合像元中第 j個亞像元的地物類別，其中i∈mp，1≤j≤s2，K表示地物類別數(shù)，在海岸線的亞像元定位過程中，一般取K=2，即水和陸兩種地物類型。融合光譜信息的亞像元定位模型表述如下：

其中，E是求取X的能量函數(shù)，Espectral是光譜能量函數(shù)，用于通過Y對X進(jìn)行光譜約束；Espatial是空間能量函數(shù)，通過空間相關(guān)性原理對X進(jìn)行空間分布上的約束；β是平衡參數(shù)，用于控制光譜能量函數(shù)和空間能量函數(shù)的權(quán)重。依據(jù)最大后驗概率理論，最有可能的亞像元定位結(jié)果：

由式可知，求解X，最關(guān)鍵的是確定Espectral和Espatial。

2.1　光譜能量函數(shù)

光譜能量函數(shù)的主要作用是確保亞像元定位結(jié)果滿足組分約束，假設(shè)低分辨率遙感影像的光譜值服從正態(tài)分布，采用最大似然估計構(gòu)建光譜能量函數(shù)[11]。Ri表示低分辨率遙感影像Y的第i個混合像元的光譜向量，μk代表第k類地物的光譜向量，θi(k)代表第i個混合像元中第k類地物所占的組分，其定義如下：

其中，1≤k≤K，li(j)∈{1,…,K }代表第i個混合像元中第 j個亞像元的地物類別。

由上述可知，第i個混合像元的光譜能量值可表達(dá)如下：

總的光譜能量值計算如下：

2.2　空間能量函數(shù)

空間能量函數(shù)主要是為了亞像元定位提供地物的空間分布信息。根據(jù)空間相關(guān)性原則，距離較近的亞像元和距離較遠(yuǎn)的亞像元相比，更加可能屬于同一種類型，這一理論已經(jīng)被廣泛應(yīng)用于各類亞像元定位算法中。根據(jù)這一原則，本文為待確定地物類別的亞像元選取一鄰域，根據(jù)鄰域中其他亞像元的地物類別來確定該亞像元的地物類別，具體模型如下：

其中，espatial(i)代表第i個混合像元的空間能量值，Ni(j)表示以第i個混合像元中第 j個亞像元為中心，選取大小為w×w的方形窗口所包含的亞像元組成的鄰域集合；p(·)表示亞像元的位置，λ(p(j),p(a))是各向同性空間權(quán)重函數(shù)，定義如下：

l(·)表示亞像元的地物類別，函數(shù)φ定義如下：

因此，總的空間能量值定義如下：

3　LA_SPSMG算法

3.1　模型的建立

由式可知，μk代表第k種地物的光譜向量，一般是通過在遙感影像上選取的若干個該類別的純像元的光譜向量，然后求均值得到的，并且在之后對每一個像元的光譜能量計算過程中不發(fā)生改變。然而，遙感影像上同種地類不同像元的光譜值可能差異較大，選取全局的陸地和海水的端元，不能充分反映這種屬于同類地物的不同像元之間的光譜差異性，會造成混合像元分解誤差，進(jìn)而影響海岸線亞像元定位的精度。

針對上述問題，本文提出了局部自適應(yīng)選取端元的亞像元定位模型。對任意混合像元i(i∈mp) 計算其陸地或海水端元光譜的過程如下：以該混合像元為中心，選取窗口（預(yù)設(shè)定大小W×W）范圍內(nèi)（如圖3所示）第k類地物的純像元，如果該窗口內(nèi)沒有該類地物的純像元，則增大W，擴(kuò)大窗口范圍，直至找到；記選取到的純像元為 pjk，相應(yīng)的光譜向量為R(pjk)，其中1≤j≤n，n為窗口內(nèi)所包含的該類地物純像元的個數(shù)，則對于混合像元i，第k類地物的端元光譜可表示為：

λ*()

pjk,i代表空間權(quán)重函數(shù)，其定義如下：

其中，Ω*為使的歸一化常數(shù)，d(pjk,i)代表純像元pjk與混合像元i之間的歐式距離，ω*為非線性參數(shù)，1≤k≤K。

將計算得到的混合像元i中第k類地物的光譜向量μi(k)帶入到亞像元定位模型中，因此，LA_SPSMG模型可表示為：

Espatial的表達(dá)式見式（16）。

3.2　模型的求解

根據(jù)式（16）、（19）和（20），最終的亞像元定位結(jié)果可計算如下：

采用模擬退火算法[12]對式（21）進(jìn)行求解。溫度T隨著迭代次數(shù)的變化如下：

圖3　鄰域窗口

其中參數(shù)σ∈(0,1)控制溫度下降的速率。

在初始化過程中，根據(jù)低分辨率遙感影像Y得到純陸像元集合pl、混合像元集合mp和純水像元集合pw，對每一個屬于mp的混合像元，在高分辨率亞像元定位結(jié)果圖X中對應(yīng)的s×s個亞像元，隨機(jī)分配地物類別。在每一次迭代過程中，如果該亞像元的類別改變使目標(biāo)函數(shù)值降低，則接受這次改變；否則，根據(jù)當(dāng)前的溫度，以一定的概率接受這次改變。當(dāng)達(dá)到最大的迭代次數(shù)后，算法運(yùn)行停止。整個模型求解的算法流程如下：

算法1 LA_SPSMG算法

輸入：低分辨率遙感影像Y，波段數(shù)L，地物類別數(shù)K=2，放大因子s，窗口大小w和W，非線性參數(shù)ω和ω*，平衡參數(shù) β，初始溫度T0，溫度遞減率σ，迭代次數(shù)t。

I）初始化：依據(jù)低分辨率遙感影像Y分別確定純陸像元集合 pl、混合像元集合mp和純水像元集合pw；然后依據(jù)pl、mp和 pw對高分辨率亞像元定位圖X進(jìn)行初始化，其中混合像元部分，在X中隨機(jī)分配亞像元的地物類別。

II）生成高分辨率亞像元定位結(jié)果圖

Forn=1:t

Fori=1:M×N

Ifi∈mp

a）以混合像元i為中心，構(gòu)建W×W的像元窗口，分別選取低分辨率遙感影像Y中的純陸像元和純水像元，并計算混合像元i中第k類地物的均值向量μi(k)(1≤k≤K )；如果窗口中沒有純像元，則擴(kuò)大W的值，直到找到純像元為止。

4　實驗設(shè)計與結(jié)果分析

4.1　實驗設(shè)計

為了驗證LA_SPSMG算法的有效性，選取兩個不同研究區(qū)域分別進(jìn)行海岸線的亞像元定位，并對實驗結(jié)果進(jìn)行定量評價。

實驗1的研究區(qū)域是位于青藏高原上的一個湖泊，如圖4所示。實驗選用Landsat ETM+遙感影像作為實驗數(shù)據(jù)，Landsat ETM+共有8個波段，其中第1～5波段以及第7波段空間分辨率為30m，第6波段空間分辨率為60m，第8波段空間分辨率為15m，多尺度的分辨率信息為實驗提供了豐富的數(shù)據(jù)資源。實驗所用的遙感影像成像時間為2001年9月22日，采用第6波段（空間分辨率為60m）包含目標(biāo)研究區(qū)域、141×142像元大小的遙感數(shù)據(jù)作為輸入的原始圖像Y，放大因子s=4，最終得到的亞像元定位結(jié)果為X（空間分辨率為15m），并利用空間分辨率同樣為15m的第8波段（全色波段）遙感影像作為驗證數(shù)據(jù)評價定位結(jié)果的精度。實驗中的其他參數(shù)如下：地物類別K=2，亞像元鄰域窗口大小w=7，搜尋純像元的窗口大小W=7，非線性參數(shù)ω=ω*=10，模擬退火初始溫度T0=0.2，溫度遞減率σ=0.982，迭代次數(shù)t=150。實驗結(jié)果如圖5所示，圖中黑色像素代表水，白色像素代表陸地；（a）為輸入的遙感影像Y；（b）為對Y通過閾值分割得到的提取結(jié)果；（c）為對Y利用基于全局端元的亞像元定位方法得到的提取結(jié)果；（d）為空間分辨率為15m的參考影像；（e）為利用LA_SPSMG算法得到的亞像元定位結(jié)果X；（f）為從參考影像中得到的提取結(jié)果。

圖4　實驗1研究區(qū)域

實驗2的研究區(qū)域是小金門島（如圖6所示），選用的Landsat ETM+遙感影像成像時間為1999年9月7日，將第3、4和5波段合成多光譜數(shù)據(jù)（空間分辨率為30m），并通過重采樣將其空間分辨率降為60m，截取包含目標(biāo)研究區(qū)域、105×108像元大小的多光譜數(shù)據(jù)作為輸入的原始數(shù)據(jù)Y，放大因子s=4，最終得到的亞像元定位結(jié)果為X（空間分辨率為15m），并利用第8波段（全色波段，空間分辨率為15m）遙感影像作為參考數(shù)據(jù)評價定位結(jié)果的精度。實驗中的其他參數(shù)設(shè)置與實驗1相同。實驗結(jié)果如圖7所示，圖中黑色像素代表水，白色像素代表陸地；（a）為輸入的遙感影像Y；（b）為對Y通過閾值分割得到的提取結(jié)果；（c）為對Y利用基于全局端元的亞像元定位模型得到的提取結(jié)果；（d）為空間分辨率為15m的參考影像；（e）為利用LA_SPSMG算法得到的亞像元定位結(jié)果X；（f）為從參考影像中得到的提取結(jié)果。

圖5　實驗1結(jié)果圖

圖6　實驗2研究區(qū)域

5.2　結(jié)果分析

如圖5所示，圖（a）是空間分辨率為60m的研究區(qū)影像，在水陸邊界處存在著大量的混合像元；圖（b）是采用閾值分割的方法得到的水陸邊界圖，閾值方法將整個混合像元全部劃分為水體或陸地，導(dǎo)致邊界呈鋸齒狀；圖（c）是采用基于全局端元的亞像元定位方法得到的提取結(jié)果；圖（d）是空間分辨率為15m的參考影像，相較于原始影像，其水陸邊界更加清晰；圖（e）是采用LA_SPSMG算法得到的亞像元定位結(jié)果，相較于閾值分割方法的結(jié)果，水陸邊界更加平滑，相較于基于全局端元的亞像元定位方法得到的提取結(jié)果，細(xì)節(jié)部分的保留更加完整；圖（f）是對參考影像利用GAC模型得到的水陸邊界結(jié)果。

圖7　實驗2結(jié)果圖

如圖7所示，圖（a）是研究區(qū)543波段合成的偽彩色圖，分辨率為60m，在海水和島嶼邊界處存在著大量的混合像元；圖（b）是采用閾值分割的方法得到的水陸邊界圖，邊界存在著許多毛刺，呈鋸齒狀；圖（c）是采用基于全局端元的亞像元定位方法得到的提取結(jié)果；圖（d）是分辨率為15m的參考影像，相較于原始影像，邊界混合像元較少，水陸邊界更清晰；圖（e）是LA_SPSMG算法的亞像元定位結(jié)果，其水陸邊界更加平滑，通過與參考影像的對比可知，該結(jié)果更接近于真實的水陸邊界；圖（f）是對參考影像利用GAC模型得到的水陸邊界結(jié)果。

為了定量評價海岸線提取的精度，將利用原始影像得到的海岸線結(jié)果與高空間分辨率的參考影像對比，采用總差異（Total Disagreement，TD）、數(shù)量差異（Quantity Disagreement，QD）、分配差異（Allocation Disagreement，AD）和地物百分比均方根誤差（Root Mean Square Error，RMSE）四個指標(biāo)對結(jié)果進(jìn)行評價。QD用于評價定位圖與參考圖之間由于地物百分比誤差造成的差異，AD用于評價定位圖與參考圖之間由于地物空間匹配誤差造成的差異，TD用于評價定位圖與參考圖之間的總體差異，其值為QD與AD之和。

QD計算如下：

其中，pli表示在被評價圖像（亞像元定位或硬分類圖像）上地物類別為l，而在參考圖像上地物類別為i的亞像元占總體亞像元數(shù)目的百分比。

AD的計算為：

TD的計算為：

地物百分比均方根誤差的計算為：

式中，N為像元總數(shù)，K為地物類別數(shù)，θlj為亞像元定位結(jié)果中第l(l∈(1,…,K))類地物在第 j(1≤j≤N)個像元中所占的百分比，ωlj為參考影像中第l(l∈(1,…,K))類地物在第 j(1≤j≤N )個像元中所占的百分比。

由表1和2可知，相比于閾值分割算法和基于全局端元的亞像元定位算法，LA_SPSMG算法無論是TD還是RMSE均最小，證明該算法在進(jìn)行海岸線亞像元定位時誤差更小、更有效。這主要是因為，LA_SPSMG算法局部自適應(yīng)地選取端元，每一個要進(jìn)行亞像元定位的混合像元都能在較小的空間范圍內(nèi)找到屬于自己的地類端元，在一定程度上減少了全局端元帶來的光譜分解誤差，提高了定位精度。

表1　實驗1各算法的誤差

表2　實驗2各算法的誤差

5　結(jié)語

本文提出了基于GAC的局部自適應(yīng)海岸線亞像元定位算法，該算法利用GAC模型提取初始海岸線，通過數(shù)學(xué)形態(tài)學(xué)的膨脹和腐蝕操作獲取需要進(jìn)行亞像元定位的混合像元集合，采用局部自適應(yīng)選取端元的海岸線亞像元定位算法對得到的混合像元集合進(jìn)行亞像元定位。實驗表明，該算法相比于采用全局端元的亞像元定位算法定位結(jié)果誤差更小、更準(zhǔn)確。

參考文獻(xiàn)：

[1]Li X,Ling F,Du Y,et al.A Spatial-Temporal Hopfield Neural Network Approach for Super-Resolution Land Cover Mapping with Multi-temporal Different Resolution Remotely Sensed Images[J].Isprs Journal of Photogrammetry&Remote Sensing,2014,93(7):76-87.

[2]凌峰,吳勝軍,肖飛,等.遙感影像亞像元定位研究綜述[J].中國圖象圖形學(xué)報,2011,16(8):1335-1345.

[3]Campbell J B.Introduction to Remote Sensing,3rd Edition[J].Crc Press,2002.

[4]Richards J A.Remote Sensing Digital Image Analysis:An Introduction[M].Springer Berlin Heidelberg,2006.

[5]張旸,陳沈良.結(jié)合遙感數(shù)據(jù)與地統(tǒng)計學(xué)方法的海岸線超分辨率制圖[J].遙感學(xué)報,2010,14(1):148-164.

[6]Ling F,Xiao F,Du Y,et al.Waterline Mapping at the Subpixel Scale from Remote Sensing Imagery with High-Resolution Digital Elevation Models[J].International Journal of Remote Sensing,2008,29(6):1809-1815.

[7]Giles M.Foody,Aidy M.Muslim,Peter M.Atkinson.Super Resolution Mapping of the Waterline from Remotely Sensed Data[J].International Journal of Remote Sensing,2005,26(24):5381-5392.

[8]Aidy M.Muslim,Giles M.Foody,Peter M.Atkinson.Localized Soft Classification for Super-Resolution Mapping of the Shoreline[J].International Journal of Remote Sensing,2006,27(11):2271-2285.

[9]Muad A M,Foody G M.Super-Resolution Mapping of Lakes from Imagery with a Coarse Spatial and Fine Temporal Resolution[J].International Journal of Applied Earth Observation&Geoinformation,2012,15(1):79-91.

[10]Kass M,Witkin A,Terzopolos D.Snakes:Active Contour Models.Int.J.Comput.Vis.,1998,(1):321-331.

[11]Kasetkasem T,Arora M K,Varshney P K.Super-Resolution Land Cover Mapping Using a Markov Random Field Based Approach[J].Remote Sensing of Environment,2005,96(3–4):302-314.

[12]Tolpekin V A,Stein A.Quantification of the Effects of Land-Cover-Class Spectral Separability on the Accuracy of Markov-Random-Field-Based Superresolution Mapping[J].IEEE Transactions on Geoscience&Remote Sensing,2009,47(9):3283-3297.