武易天, 陳 甫, 馬 勇, 劉建波, 李信鵬
(1.中國科學院遙感與數字地球研究所,北京 100101; 2.中國科學院大學,北京 100049)
近海水產養殖作為我國漁業經濟的重要組成部分,長期以來保持著高速增長勢頭。根據聯合國糧農組織2014年發布的數據,2012年我國水產品年產量占全球漁業市場的61.69%,其中藻類海產品總量占全球市場的53.97%[1]。然而,隨著近海水產養殖業的高速發展和養殖面積的不斷擴大,大規模養殖浮繩的散亂分布嚴重影響了海上交通; 加之人工飼養殘留在海水中的飼料會造成近海岸水體的富營養化,嚴重破壞近海岸區域的生態平衡,造成重大損失。因此,為了維護近海岸區域的生態環境以及當地水產養殖業的可持續發展,促進近海資源的合理利用和綜合管理,對近海岸水產養殖區域進行快速、準確的監測是很有必要的[2]。
衛星遙感技術以其周期短、重復觀測且能夠大面積獲取地物空間信息等特點,已成為地表動態監測的重要手段。利用衛星遙感影像能夠快速、準確地對近海養殖區域進行動態監測,在近海養殖區域的科學規劃與管理中發揮重要作用。馬艷娟等[3]基于ASTER衛星數據,通過分析影像內各波段對水產的敏感性差異,構建了歸一化水產指數(normalized difference aquaculture index,NDAI),并結合深海提取指數(marine extraction index,MEI)完成了近海水產養殖區養殖水體與其他水體的區分; 周小成等[4]通過構建卷積算子,使用鄰域分析法增強水產養殖區的空間紋理信息,對ASTER影像中的水產養殖區域進行了提??; 林桂蘭等[5]利用方差法并結合紋理分析,對廈門海灣的養殖區域進行提取,得到了養殖專題圖; 孫曉宇等[6]采用面向對象的方法,綜合運用養殖區的光譜、紋理、形狀特征及空間鄰近關系等信息,并結合多時相遙感數據,對珠江口海岸帶地區水產養殖場的變化情況進行了提取。然而,上述方法識別的都是養殖區的外部范圍,并未實現對養殖區內部的精確提取; 而養殖區外部往往有許多水域并未進行養殖生產,造成了這些方法的結果都出現了較大誤差,很難滿足監控管理的實際需求。盧業偉等[7]基于RapidEye高空間分辨率影像,提出了一種基于近海養殖區的特征指數并結合紋理信息的養殖區自動提取方法,對福建三都澳養殖水域進行高精度提取,很好地識別出養殖水域間的非養殖水域,大大提高了提取精度。
但是,由于近海岸水域受懸浮泥沙濃度和葉綠素濃度等的影響,常呈現出不同的反射特性,導致近海養殖區水域背景復雜多變; 加之云霧的影響,單一的水產特征指數在復雜水色背景下很難取得較好的提取效果。為了解決復雜背景下水產養殖區的高精度提取問題,本文首先利用基于正交子空間投影的約束能量最小化(orthogonal subspace projection-weighted constrained energy minimization,OWCEM)算法[8]對水產養殖區域進行增強,消除復雜背景對水產養殖區提取的干擾; 然后提出局部標準差自適應分割法(standard deviation adaptive segmentation,SDAS)對水產養殖區域進行分割; 最后基于分割結果實現近海養殖區的高精度自動提取。
為了驗證算法的適應性,本文選取福建寧德三都澳和山東榮成桑溝灣2處近海水域進行水產養殖區提取實驗。這2個區域的近海水產養殖區分布廣泛,養殖密度大且水產品種多樣,是我國著名的水產養殖模范基地。2個實驗區域在氣候和水文條件等因素上有一定的差異,適合用于測試水產養殖區提取算法的穩定性和可靠性,其分布位置如圖1所示。

(a) 山東桑溝灣養殖區 (b) 福建三都澳養殖區
圖1實驗區覆蓋范圍示意圖
Fig.1Sketchmapforcoverageofexperimentalareas
山東桑溝灣水產養殖區位于山東半島東端,除東面灣口與黃海相連外,其余三面被陸地環繞,主要養殖水產為貝殼類(箱式養殖)和海帶(浮繩式)。福建三都澳水產養殖區位于浙江省寧德市,是我國典型的近封閉型港灣,灣內水產養殖區域眾多,主要養殖種類為紫菜(浮網式)、海帶(浮繩式)和鮑魚貝殼類(網箱式)等。
不同種類的水產養殖方式有一定的差別。如圖2(a)所示,鮑魚和其他貝殼類常放在漁排下懸吊的養殖網箱中,而漁排多由木板和塑料材質搭建,其反射特性與水體差異明顯,在遙感影像中最容易提取,因此不是本文探討的重點。不同于鮑魚養殖,紫菜在海水中漂浮時間過久便會腐爛,在落潮后需要進行干燥晾曬,所以漁民多利用高度20 m左右的大毛竹結合浮網在沿海灘涂養殖紫菜(圖2(b))。而圖2(c)和圖2(d)反映的則是海帶的養殖方式,相比于紫菜,海帶在收割前并不需要晾曬,因此多將其系在長繩上進行養殖,而長繩要等間隔地系有漂浮球,使海帶不會漂浮到海水較深處,并保證其光合作用的效率。

(a) 鮑魚養殖區 (b) 紫菜養殖區 (c) 海帶養殖區1 (d) 海帶養殖區2
圖2水產養殖區不同的養殖方式
Fig.2Differentcultivationmodelsofaquacultureareas
遙感影像選取美國Landsat8衛星數據,通過美國地理地質調查局全球數據服務網站(http: //earthexplorer.usgs.gov/)下載,空間分辨率為30 m。結合水產的養殖周期,選取水產生長較為茂盛的時間節點進行水產養殖區提取實驗,其中山東桑溝灣海域選擇2016年4月15日獲取的數據,福建三都澳地區選擇2014年4月10日獲取的數據。而在精度驗證中則使用多光譜與全色波段融合后的Landsat8數據(15 m空間分辨率),以及從中國資源衛星中心獲取的2016年06月17日的桑溝灣海域高分二號(GF-2)數據和2014年05月08日的三都澳海域高分一號(GF-1)數據。
為實現復雜背景下的水產養殖區提取,本文首先利用OWCEM算法對水產養殖區域進行增強,然后結合紋理信息提出SDAS對水產養殖區域進行分割提取,從而實現近海養殖區的高精度自動提取。OWCEM- SDAS方法的整體流程如圖3所示。

圖3 水產養殖區提取方法(OWCEM-SDAS)流程圖Fig.3 Flow char of extraction method(OWCEM-SDAS) for aquaculture areas
對獲取的Landsat8原始影像進行以下預處理:
1)大氣校正。首先對Landsat8影像數據進行輻射定標,然后對其進行FLAASH大氣校正,以消除影像數據中輻射亮度的各種失真。
2)幾何糾正。幾何糾正的目的是通過一定數量的控制點,進行不同影像間的像元插值運算和空間變換,以達到消除遙感影像幾何畸變的目的。本文通過在遙感影像中均勻布設控制點,對遙感影像進行幾何精糾正。
3)水陸分離與影像裁剪。在以上基礎上,根據行政區邊界和全球30 m空間分辨率地表覆蓋產品[9]對遙感影像進行水陸分離,裁剪出實驗區域。
在近海岸水域,海水中葉綠素和懸浮泥沙濃度的分布差異使傳統方法所構建的水產特征指數很難得到水產養殖區增強的效果,會顯著影響后續的紋理分析和水產養殖區提取過程,而OWCEM算法能夠對多維特征空間進行目標信號增強和識別,并廣泛應用于高光譜遙感影像目標識別中[10]。
2.2.1 樣本選擇和波段擴展
對多光譜影像使用OWCEM算法前,需要進行波段擴展[11]。擴展的波段主要來源有2類: ①與目標地物有關的波段特征信息,在水產養殖區提取中主要由對水產較為敏感的波段構成; ②通過計算影像中其他像元與端元樣本光譜間的相似性所構成。
為了突出水和水產養殖區間的差別,本文在光譜特征空間內新擴展了3個水體指數,分別是改進的歸一化差異水體指數(modification of normalized difference water index,MNDWI)[12],以及Feyisa等[13]提出的水域自動提取指數(automated water extraction index,AWEI)——AWEInsh和AWEIsh,并增加了葉綠素濃度指數(chlorophyll-a concentration,CHL)和懸浮泥沙濃度指數(total suspended matter concentration,TSM)[14],用以減少水色差異對水產養殖區提取過程的干擾。在相關性特征空間中,為了計算端元樣本和其他像元間的相似性,選取光譜角距離SAD,相關系數corr,光譜距離s和光譜信息散度SID[15]這4個特征指數用于有效區分影像中的像元點與先驗樣本間的差異,提高水產養殖區提取的精度。上述9個指數的表達式如表1所示,其中,AWEInsh與AWEIsh采用文獻[8]的方法進行歸一化。

表1 波段擴展中所包含的特征指數Tab.1 Feature indexes included in band expansion

(續表)
表1光譜特征中的band1,…,band7分別為Landsat8多光譜數據中的第1波段,…,第7波段; 相關性特征中的x為先驗樣本的特征向量,y為影像中除先驗樣本外的所有像元點特征向量;p和q分別為x和y的概率密度分布; 光譜信息散度中的l為第l波段,l=1,2,…,L,L為影像的波段總數,這里取L=7;xl和yl分別為像元x和y中l波段的值,因為在數據預處理中已進行了輻射和大氣校正,因此這里的值代表了地物像元在不同波譜范圍上的地表反射率。
完成波段擴展后,得到9個新的特征波段,再加上Landsat8影像原始的7個波段,總計16個波段。這16個波段不僅包含了水產養殖區的特征信息,還包含了非養殖水域、懸浮泥沙和葉綠素濃度的信息,綜合分析這些特征信息可以更加容易地區分水產養殖區與非養殖水域,提高水產養殖區提取精度。
2.2.2 OWCEM養殖區增強
經過波段擴展后得到的特征空間并不能夠直接進行水產養殖區提取,還要經過OWCEM增強處理后結合自適應閾值進行養殖區分割和提取。OWCEM算法的數學表達式為
(1)
式中:N為像元總數;xi為影像中每一個像元的觀測信號,xi=(xi1,xi2,…,xiL)T,1≤i≤N;L為影像波段數(通常情況下L< (2) 式中P為進行正交子空間投影的系數,即 P=I-dd+, (3) 式中:I為單位矩陣;d為輸入水產樣本;d+為d的偽逆。 將R帶入式(1)可得 (4) OWCEM算法的具體求解過程如圖4所示。 圖4 OWCEM算法流程圖Fig.4 Flow chart of OWCEM algorithm 用OWCEM對近海水產養殖區域進行增強的步驟如下: 1)首先,由Landsat8影像選取原始水產樣本d,在波段擴展后利用式(2)計算正交投影加權自相關矩陣R。 2)其次,在求解出R后,將其帶入到式(3)中即可求得OWCEM能量最小解wOWCEM。 3)最后,將該系數與波段擴展后的影像相乘,便可得到增強后的近海水產養殖區域圖。 需說明的是,增強后的結果只是歸一化到(0,1)區間上的灰度影像,而水產養殖區域的精確提取還需要進行自適應閾值分割處理。 在大范圍水產養殖區遙感精確提取中,受懸浮泥沙等因素的影響,在全局范圍內很多養殖區與非養殖區在光譜特征上仍難以區分,這就需要結合空間信息進行進一步分割提取。由于水產養殖區呈條帶狀間隔分布特征,因此本文在使用OWCEM算法增強養殖區后,采用SDAS法對養殖區進行自動精確提取。 通過對OWCEM算法增強后的影像進行分析,發現水產養殖區呈固定寬度的條帶狀分布,其寬度在Landsat8影像中多為3個像元,而且相鄰養殖區域間被非養殖水域隔開?;谶@種紋理特征,本文采用5像元×5像元的窗口對增強后的影像進行分割,并計算各小窗口內部像元灰度分布的標準差。由于所選窗口寬度大于單條養殖區域,因此各窗口內像元分布只可能有養殖區域與非養殖區交錯分布或者該窗口內部不存在養殖區域2種情況。在養殖區與非養殖區交錯分布的情況下,分割窗口的標準差較大; 而在不包含養殖區域的情況下,局部窗口標準差則較小。因此,可以使用OWCEM處理中選取的樣本所在局部窗口的標準差作為閾值,剔除不包含養殖區的小窗口。對剩余的包含水產養殖區的局部窗口,使用最大類間方差分割法Otsu[16]進行水產養殖區的精確提取。 Otsu法是一種自適應的閾值選擇方法,通過使影像目標類和背景類二者的類間方差最大化,自動求得閾值,對影像中的目標和背景進行分割。在影像分割時,假設目標像元占影像總像元比例為w1,且目標像元類別的平均灰度為u1; 對應地,影像中背景像元所占影像總像元個數比例為w2,其類別平均灰度為u2,則影像的總平均灰度u為 u=w1u1+w2u2。 (5) 由此得到目標和背景二者的類別間方差g的數學表達為 g=w1(u1-u)2+w2(u2-u)2, (6) 進而可求得 g=w1w2(u1-u2)2。 (7) 通過選擇分割閾值T,使g最大,即 (8) 式中DNmax為影像的最大灰度值。 在局部標準差濾波后的窗口內自動選取閾值,對養殖區與非養殖區進行分割,可以得到近海養殖區的精確提取結果。相對于僅通過全局閾值提取近海養殖區的面向對象提取法,局部標準差自適應分割法能夠充分利用養殖區局部的空間分布特性,根據不同子區域自動求得分割閾值,從而很好地減少了懸浮泥沙和葉綠素濃度對分割閾值選擇的干擾,提高了近海水產養殖區的提取精度。 由于現有的提取算法鮮有基于Landsat8數據進行近海水產養殖區精確提取的,為了保證實驗的可靠性,本文根據Landsat8的數據特征,參考已有的其他衛星數據的水產養殖區提取方法,在2個實驗區進行了面向對象水產養殖區提取實驗,并與本文實驗進行對比分析。為了進一步提升面向對象方法的提取精度,在設置提取閾值時,本文在預設閾值的左右再進行10步(步長為0.02)搜索,以選取整體提取精度最好的閾值作為提取參數。 3.2.1 山東桑溝灣提取結果 圖5示出山東桑溝灣水產養殖區提取結果。圖5(a)是山東桑溝灣水域Landsat8影像的B4(R),B3(G),B2(B)波段合成影像; 圖5(b)是圖5(a)中所標注的3個局部矩形區域A,B,C的放大顯示圖; 圖5(c)是使用面向對象提取方法對山東桑溝灣海域進行水產養殖區提取后的結果,圖中黃色區域為數據預處理中通過掩模運算去掉的陸地區域,黑色區域為非養殖水域,白色條紋狀區域是水產養殖區。圖5(c)選擇了圖5(a)中同樣的3個矩形框,對局部區域進行放大,如圖5(d)所示。而圖5(e)則是用本文所提方法在山東桑溝灣水域提取的水產養殖分布結果,圖中同樣選了3個區域放大(圖5(f)),以便進行對比分析。 從圖5顯示的結果來看,影像中對應的3個區域A,B,C都有一定程度的懸浮泥沙干擾。其中區域A范圍較大,受懸浮泥沙影響最大,右側的部分條狀養殖區幾乎不可見,僅能看到模糊的輪廓; 區域B中的水產紋理較為密集受葉綠素濃度差異的影響,左下角灰度差異較為明顯; 而區域C雖然也有懸浮泥沙的影響,但該區域內懸浮泥沙濃度分布較為均一,對水產養殖區提取的影響相對較少。從目視效果上看,本文方法的提取結果(圖5(e)和(f))紋理較清晰,水產養殖區內部間隔的非養殖區水域也能夠被有效地剔除; 而面向對象的提取結果(圖5(c)和(d))則出現了較高的誤判率,其中在區域A中右側被懸浮泥沙影響的水產養殖區幾乎沒有提取出來,區域B中的結果中存在大量相連的區域,區域C中因全局閾值的影響,提取出的區域也相對有限。需要說明的是,受Landsat8影像空間分辨率的影響,矩形度等特征指數在面向對象提取過程中所起到的作用相對有限,這也增加了面向對象提取水產養殖區域的難度。 (a) Landsat8影像 (b) (a) 中3個區域放大顯示 (c) 面向對象提取結果 (d) (c) 中3個區域放大顯示 (e) 本文方法提取結果 (f) (e) 中3個區域放大顯示 圖5山東桑溝灣水產養殖區提取結果 Fig.5AquacultutreareaextractionresultinSanggouBay,Shandong 3.2.2 福建三都澳提取結果 福建三都澳地區的水產養殖區域提取結果如圖6所示。圖6中各子圖的位置與圖5類似,圖6(a),(c)和(e)分別是福建三都澳水域的Landsat8影像、面向對象提取結果和本文方法所提取的結果,圖6(b),(d)和(f)分別是圖6(a),(c)和(e)中3個矩形區域的放大圖。從影像中可以看出,福建三都澳海域受懸浮泥沙濃度的影響要比山東桑溝灣地區的小,但該景影像中有部分薄云的干擾,導致部分區域反射率升高,而左側的區域B未受云影響,水產和非水產養殖區的整體反射率較低。在提取結果中,面向對象的水產養殖區提取出現了很多錯分的情況,其主要原因在于大范圍區域內的某些非養殖水域反射率和部分養殖區反射率重疊,因此在設置水產養殖區提取的閾值時很難通過一個全局的閾值對大范圍水產養殖區進行提取。而本文方法在使用OWCEM進行目標信息增強后,在分割提取時采用2層不同尺度的提取,從目視效果上看獲得了較好的結果。 (a) Landsat8影像 (b) (a) 中3個區域放大顯示 (c) 面向對象提取結果 (d) (c) 中3個區域放大顯示 (e) 本文方法提取結果 (f) (e) 中3個區域放大顯示 圖6福建三都澳水產養殖區提取結果 Fig.6AquacultutreareaextractionresultinSanduaoBay,Fujian 比較圖5和圖6可以看出,與面向對象提取法相比,本文方法較好地利用了近海水產養殖區的光譜-空間信息。但本文方法在圖6的區域B中也存在部分錯分現象,并非所有養殖水域間的間隔(非養殖水域)都能很好地被提取出來。從整體而言,面向對象方法由于受全局閾值的影響,其整體分類結果較為一般; 而本文提出的分類方法提取的水產養殖區紋理清晰,精度較高,具有較高的實用價值。 驗證數據都使用ENVI軟件進行預處理。數據與處理算法精度的驗證是通過在每個實驗區內隨機生成一定數量的驗證樣本,然后結合同一地區的高空間分辨率影像數據(融合后的Landsat8和GF-1/2影像),對隨機生成的驗證樣本點進行目視解譯和分類。在進行目視解譯前,對GF數據進行了輻射和大氣校正,且與Landsat8數據進行了幾何配準。最后利用驗證樣本點與近海養殖區提取結果構建混淆矩陣,得到提取后的總體分類精度和Kappa系數以及水產養殖區的制圖精度和用戶精度。驗證點的選取和目視解譯參考影像如圖7所示。圖7(a)為山東桑溝灣的近海養殖區域圖,圖中綠點為非養殖水域的驗證樣本點,紅點為水產養殖區的驗證樣本點。圖7(b)為結合GF-2數據對圖7(a)中驗證樣本點進行目視解譯和分類的結果。同理,圖7(c)為福建三都澳水域的近海養殖區域圖,圖中綠點為非養殖水域的驗證樣本點,紅點為水產養殖區的驗證樣本點,而圖7(d)為圖7(c)中驗證樣本點在對應區域GF-2數據上的展示??梢钥闯鲈摼癎F-2數據受水色和云霧影響,雖然可以作為參考數據對Landsat8的驗證樣本點進行目視解譯,但較難直接用作水產養殖區的精確提取。在對2個實驗區驗證樣本點目視解譯后,對水產養殖區提取結果計算混淆矩陣,其結果分別如表2和表3所示。 (a) 山東桑溝灣養殖區影像(b) 山東桑溝灣養殖區部分驗證樣本點 (c) 福建三都澳養殖區影像(d) 福建三都澳養殖區部分驗證樣本點 圖7 驗證樣本點在實驗區參考影像中的分布Fig.7 Vertification samples of test regions distributed in reference images表2 山東桑溝灣區域提取結果混淆矩陣Tab.2 Confusion matrix of aquaculture area extraction result in Sanggou Bay, Shandong 表3 福建三都澳區域提取結果混淆矩陣Tab.3 Confusion matrix of aquaculture area extraction result in Sanduao Bay, Fujian 由表2和表3可以看出,面向對象方法在提取時將大量的非養殖水域像元誤提取為水產養殖區,而本文方法誤檢的像元個數相對較少; 在非養殖水域的識別上,2類方法表現都較好,誤檢像元較少。 2個實驗區2種提取方法的整體提取精度結果如表4所示。 表4 面向對象和本文方法提取精度對比Tab.4 comparison of extraction accuracy using obeject-oritented method and proposed method (%) 從表4可以看出,在山東桑溝灣的水產養殖區提取過程中,面向對象方法用戶精度達到88.78%,但其制圖精度很低,Kappa系數僅為57.03%,總體精度僅為69.69%; 本文方法提取精度較高,總體精度達到了94.86%。在福建三都澳實驗區,面向對象提取精度有所提高,但總體精度仍未達到80%; 而本文方法雖然精度相比山東桑溝灣實驗區有所下降,但Kappa系數仍在90%以上,總體精度也比面向對象方法高出15%。 綜上所述,在山東桑溝灣和福建三都澳水產養殖區提取實驗中,面向對象提取算法總體精度不足80%,而本文算法總體精度維持在93%以上。證明本文方法能夠在復雜水色背景下,對水產養殖區進行精確提取。 利用遙感技術進行近海養殖區的提取已經廣泛地應用于生產實踐,但目前傳統的提取方法在海域水色較復雜的情況下提取養殖區的精度較低,遠遠不能滿足實際應用的需求。本文在充分利用了近海養殖區的光譜和空間信息的基礎上,利用正交子空間投影的約束能量最小化(OWCEM)方法,對復雜海域下的水產養殖區進行增強; 并通過局部標準差自適應分割(SDAS)算法實現了養殖區的自動分割和提取,大大提高了養殖區提取結果的精度。與現有的近海養殖區提取方法相比,本文方法具有以下優點: 1)充分利用了影像中所有的波段信息,能夠在復雜水色背景下有效增強近海水產養殖區,降低復雜背景對提取的影響。 2)自動化程度高。除算法初期需要輸入少量樣本點(不超過10個)外,整個處理流程完全自動進行,大大減少了提取的工作量,也降低了方法使用難度。 3)基于局部標準差自適應分割方法自動調整每個子區域的分割閾值,解決了全局分割閾值不能精確提取水產養殖區的問題。 但本文方法還存在某些不足,在水產養殖區提取過程中未考慮水產生長情況對提取的影響,導致部分區域出現了一些錯分的現象,這也是下一步進行深入研究的目標,即結合多時相遙感數據更有效地對近海養殖區進行提取并監控水產的生長情況。
2.3 SDAS法養殖區提取
3 實驗結果與精度檢驗
3.1 對比實驗的選擇
3.2 實驗結果分析






3.3 精度驗證





4 結論