基于GAC的局部自適應海岸線亞像元定位算法

陳飛宇，胡友彬，施恩

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

0　引言

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

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

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

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

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

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

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

1.1　初始海岸線提取

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

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

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

其中，Gσ是方差為σ的高斯核函數，因此，g為關于圖像梯度強度的遞減函數。

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

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

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

1.2　混合像元的確定

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

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

則需要進行亞像元定位的混合像元集合表示為：

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

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

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

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

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

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

由式可知，求解X，最關鍵的是確定Espectral和Espatial。

2.1　光譜能量函數

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

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

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

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

2.2　空間能量函數

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

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

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

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

3　LA_SPSMG算法

3.1　模型的建立

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

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

λ*()

pjk,i代表空間權重函數，其定義如下：

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

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

Espatial的表達式見式（16）。

3.2　模型的求解

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

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

圖3　鄰域窗口

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

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

算法1 LA_SPSMG算法

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

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

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

Forn=1:t

Fori=1:M×N

Ifi∈mp

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

4　實驗設計與結果分析

4.1　實驗設計

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

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

圖4　實驗1研究區域

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

圖5　實驗1結果圖

圖6　實驗2研究區域

5.2　結果分析

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

圖7　實驗2結果圖

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

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

QD計算如下：

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

AD的計算為：

TD的計算為：

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

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

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

表1　實驗1各算法的誤差

表2　實驗2各算法的誤差

5　結語

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

參考文獻：

[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].中國圖象圖形學報,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]張旸,陳沈良.結合遙感數據與地統計學方法的海岸線超分辨率制圖[J].遙感學報,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.