梁超,劉利,劉建強,鄒斌,鄒亞榮,崔松雪
(1. 自然資源部 國家衛星海洋應用中心,北京 100081;2. 自然資源部 空間海洋遙感與應用研究重點實驗室,北京 100081;3. 中國科學院 空天信息創新研究院,北京 100094)
紅樹林是指生長在熱帶、亞熱帶近岸潮間帶上部灘涂淺灘,以紅樹植物為主體的常綠灌木或喬木組成的潮灘濕地木本生物群落,是陸地向海洋過渡的特殊生態系統。紅樹林具有防風消浪、促淤保灘、凈化海水、保持生物多樣性等重要功能。近幾十年來,由于圍海造地、圍海養殖、砍伐等人為不合理開發活動等因素影響,使得紅樹林面積逐年銳減。紅樹林生態環境保護工作對長期監測紅樹林的分布范圍、生物量、健康水平等提出了要求,其中紅樹林分布范圍的監測是基礎。利用遙感手段可以大范圍、快速、重復監測紅樹林時空分布信息,這對于研究保護紅樹林生態系統具有重要意義。
國內外紅樹林遙感監測方面開展了大量研究,所采用的數據源通常為NOAA-AVHRR、Landsat、SPOT、Quickbird、ALOS/PALSAR等衛星以及諸如HyMap等航空高光譜遙感影像,主要方法包括人工解譯、監督分類、植被指數及紋理特征、面向對象分類等[1-14]。HY-1C是我國第一代海洋水色衛星HY-1A/B的后繼星,于2018年9月7日在太原衛星發射中心發射成功,是我國首顆業務化運行、全球開機、連續觀測的海洋水色衛星,載荷配置和性能指標、探測能力等均有明顯改進,可獲取重訪周期更短、觀測范圍更大、精度更高的水色遙感產品。HY-1C搭載的海岸帶成像儀(Coastal Zone Imager,CZI)用于獲取海陸交互區域實時數據,監測近海、海島、海岸帶等信息,其幅寬、分辨率相比上一代載荷分別提升近2倍和5倍。本文擬以廣西山口紅樹林生態自然保護區為研究區,綜合利用紅樹林與一般陸地植被光譜特征及空間分布特征,選用適當的光譜指數重構原始數據,經最小噪聲分離變換(Minimum Noise Fraction Rotation,MNF)后,建立決策樹提取紅樹林,據此探討HY-1C衛星CZI數據在紅樹林監測工作中的適用性。
本文以廣西山口紅樹林生態自然保護區為研究區(圖1),該保護區由位于廣西合浦縣東南部沙田半島東、西兩側的海域、陸域及全部灘涂組成,地處亞熱帶,是我國第二個國家級的紅樹林自然保護區,分布著發育良好、結構典型、連片較大、保存較完整的天然紅樹林,有紅海欖、木欖、秋茄、桐花樹等15種紅樹林植物。山口紅樹林國家級自然保護區具典型的大陸紅樹林海岸生態系統特征,紅樹林中還棲息著多種海洋生物和鳥類,具有重要的科學價值。

圖1 研究區位置與HY-1C衛星CZI原始影像Fig. 1 Study area location and the CZI original image of HY-1C satellite
采用2018年10月31日11時25分獲取的HY-1C衛星CZI傳感器L1B級數據,CZI主要技術指標見表1。

表1 HY-1C衛星CZI傳感器主要技術指標Table 1 The main technical indicators of HY-1C satellite CZI sensor
本文采用波段重構與光譜變換相結合的方法,以增強CZI影像上紅樹林與一般陸地植被的特征差異,提高二者可分性。圖2為研究方法流程圖,首先對CZI數據預處理,然后基于紅樹林與一般陸地植被光譜特征構建多種光譜指數,通過相關性分析篩選最佳指數組合重構原始數據,繼而采用MNF進一步增強重構數據中紅樹林與一般陸地植被的光譜及空間分布特征差異,最后建立決策樹提取紅樹林。

圖2 本文研究方法流程圖Fig. 2 The research flow chart
2.3.1 數據預處理
數據預處理包括幾何校正和輻射定標,CZI數據采用HDF5格式,幾何校正利用數據自帶的RPC參數。L1B數據像元值為輻射亮度,據公式(1)轉換為表觀反射率[15]:

式中,ρ為大氣層頂表觀反射率(無量綱),Lλ為大氣層頂進入衛星傳感器的光譜輻射亮度,單位為W/(m2·sr·μm),D為日地距離 (天文單位),ESUNλ為大氣層頂的平均太陽光譜輻照度,單位為W/(m2·μm),θ為太陽天頂角,相關參數可由HDF5數據集屬性參量中查找。
2.3.2 光譜指數
紅樹林與一般陸地植被光譜特征相似,在CZI原始影像難以區分。為突出植被、水體等信息并增強紅樹林與一般陸地植被的差異,本文擬選用適當的光譜指數重構CZI原始數據。綜合考慮大氣、紅樹林分布環境等因素的影響,本文選用了常見的8種植被指數[16-19]和1種水體指數,為進一步突出CZI影像上植被光譜特征,構造了3個波段比指數和1個光譜斜率比指數,其中波段比指數用來表達目標物在不同波段之間的歸一化差值比,可見光光譜斜率比指數用來表達不同地物在可見光波段的光譜曲線形態差異。表2列出了本文所采用的13種光譜指數。
2.3.3 最小噪聲分離變換
MNF最早由Green等[20]提出,是在主成分變換基礎上提出的一種線性變換算法。主成分變換各分量按照方差降序排列,影像主要信息集中在前幾個分量,但不能保證影像質量(信噪比)也按照降序排列,MNF用以改進主成分變換在噪聲消除及影像增強上的不足。Boardman和Kruse[21]證明了MNF等效于連續兩次主成分變換,設多光譜數據由加性噪聲和信號兩部分組成,首先對影像噪聲協方差矩陣進行變換,使變換后噪聲協方差矩陣為單位矩陣且波段間不相關,然后對變換后的數據集作標準主成分變換[22]。
第一步,對多光譜影像X的噪聲協方差矩陣進行對角化:

式中,Cxn為原始影像噪聲協方差矩陣,Λ為其特征值組成的對角陣,E為由其特征向量組成的正交矩陣。
上式前乘以 (Λ-0.5)T,后乘以Λ-0.5,則

式中,I為單位矩陣,令F=EΛ-0.5,則上式轉換為

原始影像經Y=FTX變換后的協方差矩陣為

表2 本文使用的植被光譜指數列表Table 2 List of vegetation spectral indices used in this paper

式中,Cx為原始影像協方差矩陣。
第二步,對Y應用標準主成分變換:

式中,ΛY是CY特征值降序排列的對角陣,B是其特征向量組成的正交矩陣。
綜上,MNF矩陣為

多光譜數據經過MNF后各分量按照信噪比降序排列,影像有效信息集中在前幾個分量,噪聲主要存在于后面的分量中,實現了光譜增強的同時分離數據中的噪聲,改善多光譜處理結果。
MNF最重要的一步是估計原始影像噪聲協方差,其方法一般借助影像的空間特征,據Canty[23]:

式中,x=(i,j)T表示影像X的一個像元坐標,h=(h1,h2)T表示坐標的偏移量,即噪聲協方差可以用原始影像與偏移影像差異的協方差來估計。
通過建立感興趣區,提取紅樹林與一般陸地植被樣本的表觀反射率光譜曲線(圖3),可見在CZI影像上二者均表現出典型的綠色植被光譜特征,即綠波段強反射、紅波段強吸收以及近紅外波段的高反射率,藍波段由于大氣散射的影響亦表現出較強的反射,而近紅外波段上二者光譜曲線均較為發散,這是由于中低分辨率影像上混合像元效應的影響。紅樹林與一般陸地植被光譜曲線具有較高的相似性,使得二者很難在4波段CZI數據中被區分。圖4為CZI影像上紅樹林與一般陸地植被表觀反射率分布直方圖,可見二者在各波段具有大致相近的直方圖形態和峰值分布區間,雖然在第4波段,紅樹林與一般陸地植被直方圖峰值位置略有差異,但如前述,紅樹林與陸地植被在該波段光譜曲線均較為發散,直方圖上分布范圍均較寬,從而存在較大的重疊區。綜合考慮光譜曲線及直方圖可知,單純依靠CZI表觀反射率數據很難區分紅樹林與一般陸地植被。

圖3 CZI影像上紅樹林(a)與一般陸地植被(b)光譜曲線Fig. 3 Spectral curves of mangrove (a) and general terrestrial vegetation (b) on CZI image

圖4 紅樹林(a)與一般陸地植被(b)表觀反射率分布直方圖Fig. 4 Histograms of apparent reflectance distribution of mangroves (a) and general terrestrial vegetation (b)
紅樹林與一般陸地植被在CZI原始光譜空間中較高的相似性對紅樹林檢測帶來了難度,傳統的監督/非監督分類等方法提取紅樹林效果均不甚理想。本文通過構建對植被和水體信息敏感的光譜指數組合代替原始波段數據,以增強影像中典型植被的光譜信息及其差異性。
本文以植被遙感中最常用的歸一化差異植被指數NDVI為參考基準,開展了NDVI與其他12種光譜指數的相關性分析,表3為各指數之間的相關系數。首先剔除與NDVI相關系數極低的光譜指數,如CBRI,因為極低的相關系數意味著該類指數對目標植被信息不敏感;其次,剔除與NDVI具有較高的正相關系數的光譜指數,如 SIPI、CNBI、MSAVI、EVI等,該類指數在表征植被信息方面與NDVI作用相當,從避免冗余考慮予以剔除;最后,選擇與NDVI具有中、低相關性及最大負相關性的光譜指數參與重構,如 RI、NDWI、NDGI、CBGI、ARVI和 CVSSR 等,該類指數既能從不同程度反映植被信息,又避免了信息冗余,確保重構后的新光譜空間既能最大化地增強植被及水體信息,又能最大程度減少波段之間的相關性。據此原則并經多次試驗,本文最終選取NDVI、NDWI、ARVI和CVSSR為最佳指數組合,并以之依次代替CZI 4個原始波段組成新的多光譜數據。

表3 NDVI與其他光譜指數間的相關系數Table 3 Correlation coefficients between NDVI and the other spectral indices
圖5為重構后多光譜數據中紅樹林與陸地植被的直方圖,對比圖4可見在新光譜空間中,第1波段上紅樹林與陸地植被分布大致相近,但紅樹林分布范圍更寬一些;第2波段上紅樹林的分布峰值位置約為-0.50,而一般陸地植被則小于-0.50;第3波段上二者分布特征最為接近;第4 波段上紅樹林分布峰值位置約為0.70,較之一般陸地植被偏大。這說明在重構數據光譜空間中紅樹林與一般陸地植被的差異性得到一定程度增強。
對重構后數據進行MNF,圖6為變換結果的4個分量圖像,根據MNF原理,變換后4個分量按照噪聲水平升序排列,第1、2、3分量包含了主要的空間信息,而第4分量則包含較多的噪聲信號。特別地,在MNF第3分量上,紅樹林與一般陸地植被表現出較明顯的目視差異,前者具有更高的像元值。由圖7可知,MNF后,紅樹林與一般植被直方圖分布形態差異性得到提升,峰值分布區間的重合度進一步減小。具體表現為:第1分量上,紅樹林分布范圍更寬,峰值位置約在9.30附近,陸地植被分布范圍集中在10以上,峰值位于11.52附近;第2分量上,紅樹林分布范圍亦較寬,峰值位置約在-1.23附近,陸地植被接近均值為-2.13的正態分布,主要范圍在-3.77至1.77之間;第3分量上紅樹林與一般陸地植被分布差異最為明顯,紅樹林分布中心偏右,峰值位置約為5.13,陸地植被分布峰值位置大致為3.62,且絕大部分分布在約4.36以下,而紅樹林一般高于此值(見圖7中黑色虛線);第4分量包含更多噪聲信號,二者區分不大。這說明,基于光譜指數重構數據,MNF在分離了噪聲信號的同時,由于利用了影像空間分布信息(公式(8)),進一步增強了目標地物的波譜差異,提高了紅樹林與一般陸地植被的可區分度。

圖5 光譜指數重構數據樣本直方圖Fig. 5 Histograms of the spectral index reconstruction data

圖6 重構后數據MNF的4個分量Fig. 6 The 4 MNF components of the reconstruction data

圖7 光譜指數重構數據MNF結果樣本直方圖Fig. 7 Histograms of MNF components of the spectral index reconstruction data
通過上述分析,基于光譜指數重構數據MNF結果可以較容易地構建分類決策樹,實現紅樹林提取。圖8為本文構建的決策樹,其中:
規則 1:NDVI > 0.50,區分植被與非植被。
規則2:MNF3 > 4.36,區分紅樹林與一般陸地植被。
國家衛星海洋應用中心在諸如“908”專項調查、廣西海域使用本底庫建設等工作中在該區域開展過多次遙感調查及現場核查,圖9a是基于歷史調查成果和專家經驗支持下的山口紅樹林分布范圍人工解譯結果,圖9b是基于本文決策樹的紅樹林自動提取結果,可見絕大部分紅樹林均能被正確提取,其中專家解譯的面積為9.44 km2,本文方法提取紅樹林的面積為9.65 km2,面積檢測相對準確率可達97%。進一步在研究區內隨機選取2 000個驗證樣本點,其中覆蓋紅樹林樣本點50個,本文方法檢測出樣本點44個,從而紅樹林總體檢測精度約為88%。

圖8 紅樹林提取決策樹Fig. 8 Decision tree for mangrove extracting
本文基于HY-1C衛星CZI 4波段影像,以廣西山口紅樹林自然保護區為研究區,采用光譜指數重構CZI數據,經MNF建立決策樹,實現了研究區紅樹林自動提取,研究結果表明:利用NDVI、NDWI、ARVI、CVSSR 4種指數重構CZI數據,減少了地形、大氣、傳感器等因素對植被光譜信息的影響,且重構后數據對植被、水體等更加敏感。MNF則進一步利用了影像上紅樹林與陸地植被的波譜及空間分布特征,變換結果中二者差異被顯著增強,最終通過建立決策樹提取了紅樹林分布信息。本文結合光譜指數重構和MNF方法,較好地解決了在CZI原始光譜空間中紅樹林與陸地植被難以區分的問題。因此,HY-1C CZI數據可以有效地用于紅樹林空間分布監測。
HY-1C CZI數據優勢在于重訪周期短、幅寬大,可快速重復獲取大范圍中等分辨率光學影像數據,但因波段較少,缺少對研究土壤與水分情況較重要的短波紅外、中波紅外或熱紅外等波段,且波段較寬,對典型植被波譜特征細節描述能力稍弱。本文研究方法可彌補CZI數據用于紅樹林監測的上述不足,后續還可通過形態學方法對提取結果進一步優化,提升自動檢測準確率。需要指出的是,本文決策樹規則閾值對于不同地區和時相的數據可能會有所調整,目前HY-1C衛星尚處于在軌測試階段,相信衛星正式交付后,傳感器數據的定標精度和質量穩定性會得到進一步保障,有利于后續對本文方法的驗證和完善。

圖9 研究區紅樹林提取結果分布Fig. 9 Mangrove extraction results