基于GEE的廣西北部灣沿海水產養殖池塘遙感提取

文 可，姚煥玫，黃 以，陳華權，廖鵬任

（廣西大學資源環境與材料學院，南寧 530004）

0 引 言

隨著人口和糧食需求的不斷增加，水產養殖業在過去幾十年中大幅增長，特別是在中國沿海區域，水產養殖池塘已顯著擴大[1]。養殖作為廣西北部灣海岸帶的傳統產業，在地區生產中占有重要地位[2]，但在人類活動增強的影響下海岸帶生態環境不斷惡化，成為生態脆弱和災害頻發的重點區域[3]。有研究發現水產養殖池塘的擴張與河口和沿海水域的淤積、富營養化直接相關[4]。此外，為了增加養殖產量、預防疾病，過量的飼料以及抗生素也會造成嚴重的環境問題，例如沿海海域的水污染、抗生素殘留以及赤潮[5-7]。因此，了解水產養殖池塘的空間分布對于海岸帶的科學管理和漁業可持續發展具有重要意義。

衛星遙感技術具有大尺度、低成本和實時的特點，是監測和研究沿海環境的有效手段[8]，聯合國糧食及農業組織（FAO，Food and Agriculture Organization）也提議使用遙感和地理信息系統來加強水產養殖管理，近年來水產養殖區的遙感識別已成為海岸帶生態環境領域重要的研究方向之一。SPOT、WorldView-2、PlanetScope和GF-2等衛星數據因高分辨率優勢能夠在一定程度上克服混合像元對水體信息提取精度的影響，提取結果更接近養殖實際水面面積，但因影像的高成本以及衛星掃描帶寬限制，多用于單一海灣的小區域[9-11]。而Landsat系列、Sentinel-1和Sentinel-2等中分辨率遙感數據因圖幅優勢，適用于提取更大范圍的水產養殖池塘。部分研究基于養殖區域在遙感影像上的解譯標志，結合經驗和研究區域相關資料，進行目視判斷識別養殖池塘[12]，但人工成本過高且不可復制，難以推廣；部分研究通過計算特定光譜指數，增強目標地物與其他地物光譜之間的差異進行提取[13-15]；為提高識別養殖池塘的效率和準確性，一些研究提出了基于像元或者對象的養殖池塘自動提取方法，其中面向對象分類的方法被廣泛應用于識別水產養殖區域[16-21]。然而，因影像分辨率以及識別方法性能限制，目前研究在識別規模化集聚型養殖區域具有較好效果，但在識別沿岸線散亂分布的小型池塘時效果不佳。相關結果多以土地利用類型調查為導向，池塘與相鄰的堤防、其他地表水體難以分割，為養殖區域總面積，與養殖實際水面面積存在較大差異。其次，研究多以單一日期影像識別養殖池塘，未考慮廢棄池塘以及季節性水域水淹時與養殖池塘在光譜特征上高度相似，容易誤識，影響結果的準確性，僅靠單日影像難以準確識別出實際進行養殖生產活動的水體。綜上所述，大范圍的水產養殖池塘高精度提取仍然面臨挑戰。

谷歌地球引擎（GEE，Google Earth Engine）平臺的普及突破了傳統遙感數據處理方法應用于大范圍區域的工作量限制，極大降低了大數據分析系統的使用門檻[22]。本研究基于GEE平臺加載2019年全年的Sentinel-2、Sentinel-1遙感數據，提出一種適用于大范圍復雜環境的高精度水產養殖池塘識別方法，以期識別廣西北部灣海岸帶的水產養殖池塘空間分布，并結合亞米級Google高清影像評估提取方法的準確性。

1 研究區域及數據來源

1.1 研究區域

廣西北部灣地處北回歸線以南的低緯度地區，南瀕熱帶海洋地區，受海洋性季風影響，屬于熱帶季風氣候，東自洗米河口與廣東接界，西至北侖河與越南分界，涵蓋北海、欽州和防城港3個地級市[23]。因優越的自然條件，北部灣一直是我國最適宜開展水產養殖的地區之一，加之《廣西北部灣經濟區發展規劃》獲批，更進一步促進了廣西水產養殖業的發展，但隨著養殖規模擴大，在獲取經濟效益的同時，生態環境壓力也在不斷增長[24]。為探究廣西北部灣海岸帶的水產養殖池塘空間分布，基于珍珠灣至安浦港岸線生成5 km緩沖區作為研究區域。如圖1所示，研究區域海灣狹窄，地形錯綜復雜，此外，養殖用地還與其他土地利用類型存在競爭，斑塊分散，導致養殖池塘識別結果相較于中國北方平原地區具有更大的不確定性[14]，以廣西北部灣為例也可以更好檢驗本研究方法的性能。

1.2 遙感數據

本研究使用歐盟和歐洲航天局哥白尼計劃提供的Sentinel-1 C波段合成孔徑雷達地距檢測（SAR GRD，C-band Synthetic Aperture Radar Ground Range Detected）遙感數據以及Sentinel-2多光譜（MSI，Multispectral Instrument）遙感數據。Sentinel-1 SAR GRD產品由GEE平臺收集，已進行軌道復原，熱噪聲去除，地形校正和輻射定標預處理，包括兩顆衛星，重訪周期為6 d，地面采樣距離為10 m。在GEE平臺加載了2019年全年覆蓋研究區域的Sentinel-1 SAR GRD的VH交叉極化數據，共計148景。Sentinel-2 MSI包括2A、2B兩顆衛星數據，重訪周期為5 d，空間分辨率為10 m，加載2019年覆蓋研究區域的Sentinel-2 MSI（Level-1C）數據，共計1 000景，并利用GEE提供的Sentinel-2云概率產品對數據集進行去云預處理。大部分研究區域在去云處理后仍能獲得70次以上的良好觀測數據，而同期Landsat-8數據獲取量可能不足10次，Sentinel-2高時空分辨率的優勢能夠提供穩定的時序觀測數據以監測地物年內動態變化，極大提高了水產養殖池塘識別提取的準確度。

2 方 法

2.1 訓練樣本

研究區域土地覆蓋類型主要為農用地、植被、不透水面以及水產養殖池塘等。為分析各土地覆蓋類型的光譜特征，基于2019年Google高清影像目視解譯繪制了水產養殖池塘水面和堤防采樣點位各1 400個；自地球大數據共享服務平臺（http://data.casearth.cn/）獲取了劉良云等[25]創建的30 m精細地表覆蓋產品，并隨機生成1 250個不透水面點位，以及2 000個其他類型（農田、灌叢、草原、森林等）點位，由于其分類結果可能存在誤差，進一步結合Google高清影像進行了檢驗，刪除誤分點位并就近補充，保證訓練樣本點位的準確性。最后，整合所有訓練樣本點位，分布如圖2所示。

2.2 水產養殖養殖提取方法

水產養殖池塘提取流程如圖3所示，主要包括6個步驟：1）基于遙感數據計算分類特征值；2）基于訓練樣本確定最佳分割閾值；3）生成分類對象；4）基于對象面積進行再分類，堤防識別效果不佳的對象采用再分割進一步消除誤差；5）剔除其他地表水體；6）準確性檢驗。

2.2.1 水淹頻率

1）水體識別

水產養殖池塘與地表水體具有相似的光譜特征，因此將識別工作細分為水體與非水體的識別。歸一化差異水體指數（Normalized Difference Water Index，NDWI）[26]、修正的歸一化差異水體指數（modified Normalized Difference Water Index，mNDWI）[27]和自動水體提取指數（Automated Water Extraction Index，AWEI）[28]被廣泛用于水體提取。但在Sentienl-2影像中，計算mNDWI、AWEI所需的短波紅外波段（SWIR，Shortwave Infrared）空間分辨率為20 m；而計算NDWI所用的綠光和近紅外（NIR，Near Infrared）波段空間分辨率均為10 m，能夠更好地區分堤防與池塘實際養殖水面，提取效果更好。NDWI計算公式如下

式中ρGreen、ρNIR分別為Sentinel-2影像的B3、B12波段。

水產養殖池塘會在冬、春季清塘排水[29]，此時無法基于水體識別方法進行提取，為消除這一影響大多研究使用4—10月的影像數據以避開池塘干涸期[14,19-20,30]。但廣西沿岸存在水產養殖池塘廢棄的現象，廢棄池塘以及季節性水域水淹時與實際養殖池塘在光譜特征上高度相似（圖4c），易產生誤識，僅靠單日影像難以準確識別出進行養殖生產活動的水體。水體像元NDWI往往>0，如圖 4a、b所示。養殖池塘一年中因長期儲水，NDWI>0的頻率高；而廢棄池塘只在漲潮時被淹沒并不蓄水，NDWI>0的頻率低。基于此，可通過計算時序遙感數據-水淹頻率（IF，Inundation Frequency）區分出實際養殖池塘。

2）計算水淹頻率以像元NDWI是否>0作為水體、非水體判別標準，計算各像元2019年水淹頻率IF，計算公式如下

式中Nwater為像元被識別為水體的次數，N為像元的有效觀測次數。

通過分析訓練樣本的IF數值分布，設置分割閾值為0.4較合適。水產養殖池塘因長期儲水，99.7%的池塘點位IF>0.4；農用地、植被等其他類別樣本IF分布最為集中，超過98%的其他點位IF處于[0，0.1]區間。通過檢測研究區域IF柵格數據也發現，以IF作為分類特征值可以很好剔除低含水量區域，同時根據IF大小可區分出養殖池塘與季節性水域、灘涂、水稻田等動態變化的含水土地覆蓋類型。但如圖5圈出區域所示，利用IF提取養殖池塘仍面臨兩個困難：1）由于NDWI存在一定局限性，無法準確區分水體與瀝青路、陰影等低反射率地表[31]，導致計算得出的IF出現誤差（圖5a）；2）在一些高密集的養殖區域中池塘間堤防細小難以準確識別（圖5b）。因此，僅靠特征值IF無法準確提取水產養殖池塘，需加以修正。

2.2.2 水淹頻率的誤差修正

1）修正不透水面及陰影誤差

瀝青路、陰影等低反射率地表在城區出現頻繁，導致在該區域中利用IF提取水體易出現誤識。SWIR已被廣泛應用于不透水面的遙感識別[32]，針對不透水面誤差引入Sentinel-2的“B12”第二短波紅外波段（SWIR2）加以修正，計算時間序列中按波段數值升序排列的10%～90%區間的平均值（SWIR2mean(10%～90%)）。針對陰影誤差引入Sentinel-1 SAR GRD的交叉極化VH數據，計算年度平均值（VHYEARmean），該數據具有較好的水體提取能力[33]，且不易受建筑陰影噪聲影響。

通過分析訓練樣本點位處的特征值分布，在不影響識別目標水產養殖池塘的原則下，確定SWIR2mean(10%～90%)和VHYEARmean分割閾值分別為0.09和-18。但以訓練樣本為參考，“SWIR2mean(10%～90%)<0.09”并不能完全分割出不透水面，檢查遺漏點位發現多為城區建筑陰影，且VHYAERmean均高于-18。基于此，通過整合兩個特征值互補以修正不透水面及陰影帶來的誤差，如圖6a為僅使用IF進行提取的效果，圖6b為引入SWIR2mean(10%～90%)、VHYEARmean后的修正效果，圈出誤識部分得到極大改善。

2）修正堤防誤差

針對IF識別堤防效果不佳的問題，計算NDWI在時間序列中按波段數值升序排列的85%～95%區間的平均值加以修正。相較于年度平均值、中位值等，NDWImean(85%～95%)為地表含水量較大時狀態，此時池塘周邊的堤防、水壩等不含水像元會顯得更暗，水體與非水體的可分離性更強，在區分池塘間堤防上具有顯著優勢。在不影響提取養殖池塘前提下，設置NDWImean(85%～95%)分割閾值為0.12。圖7展示引入NDWImean(85%～95%)后堤防的去除效果，堤防邊界更為清晰。但是，單一分割閾值應用于大范圍復雜環境中存在一定局限性，部分區域分割效果并不理想，當閾值過小會導致養殖池塘與堤防、河流、近岸海域、排水渠道等直接相連，不能作為單獨的目標被提取，需適當上調閾值。

為篩選出需上調閾值的識別目標，基于養殖池塘的二值化圖像使用連通分割算法生成分類對象，如果像元在二值化圖像中具有相同值且連通（4向）則被劃分為同一對象，并計算其面積、周長和景觀形狀指數（LSI，Landscape Shape Index），LSI計算公式如下

式中Perimeter為對象周長，m；Area為對象面積，m2。

如圖8a所示，當多個池塘因連通被劃分為同一對象時面積會顯著偏大，因此，基于對象面積大小進行分類的方法，可用于解決堤防不易識別問題。常規養殖池塘面積一般為0.002～0.006 km2，而大型池塘則在0.02～0.03 km2左右。圖8b篩選出面積大于0.03 km2的對象，由于包含大量堤防，相較于養殖實際水面存在較大誤差，經測試在這些對象中將NDWImean(85%～95%)分割閾值適當上調至0.2，進行再次分割后效果顯著，如圖8c所示，堤防引起的誤差得到極大改善。同時干流、近岸海域等大型天然水體被分割為獨立對象，可以通過計算面積、周長和LSI快速剔除。綜上所述，基于所有分類特征值的數值分布規律將養殖池塘識別方法設置為：[(IF>0.4)和(VHYEARmean<-18)和(SWIR2mean(10%～90%)<0.09)和(NDWImean(85%～95%)>0.12或0.2)]。

2.3 剔除其他細小水體

識別結果中還包含細小支流、湖泊和排水渠道等其他非目標水體，它們在空間形狀上與養殖池塘有所差異，有研究通過設置LSI閾值加以區分[20]，該方法在規模化集聚型養殖區域效果較好。但在廣西北部灣沿岸存在大量散亂分布的小型池塘，因分辨率限制識別結果仍有部分池塘連接，導致LSI產生偏差，僅依靠LSI閾值進行篩選易產生誤分。通過檢視養殖池塘識別結果，發現不同區域池塘的空間形態也有所差異，例如茅尾海的養殖池塘平均面積為0.005 km2，平均周長為363 m，平均LSI為1.37；而廉州灣養殖池塘較小，平均面積、周長和LSI依次為0.003 km2、274 m和1.35。為剔除其他水體，首先按海灣將所有對象分類，計算區域平均面積、周長和LSI，并根據數值分布劃定各海灣區域合理值域，三個參數都異常的對象直接刪除；對于僅有一個參數接近均值的存疑對象，則進一步結合Google高清影像進行目視解譯；對于小部分仍與其他水體連接的池塘則依次進行空間分割操作，并刪除非目標部分。

3 結果與討論

3.1 水產養殖池塘空間分布及準確性檢驗

結果顯示廣西北部灣海岸帶水產養殖池塘面積共計199.3 km2，其中北海養殖面積最大為112.9 km2，防城港和欽州分別為44.2 km2、42.2 km2。如圖9所示，規模化集聚型養殖區主要分布在安鋪港、北海市下部、廉州灣、茅尾海上部以及珍珠灣左側，其中廉州灣沿岸養殖池塘分布最為密集；其他養殖池塘則散亂分布在廣西北部灣沿岸的河流入海口和灘涂區。

通過將提取結果分為養殖池塘和非養殖池塘，并隨機生成驗證點計算混淆矩陣，這是目前相關研究進行準確性檢驗的主要方法[14,19,20]。因識別方法、影像空間分辨率限制，水產養殖池塘提取誤差往往出現于水面邊界處，當點位距離池塘過遠往往難以起到有效檢驗的作用，此類驗證點如過多還會在一定程度上高估結果準確性，因此明確驗證點位空間位置對于準確性檢驗至關重要。本研究采用兩種更嚴格的檢驗方法：1）基于養殖池塘識別結果生成20 m緩沖區；在識別結果內隨機生成500個理論池塘點位，在緩沖區非目標部分隨機生成500個理論非池塘點位，以更高頻地檢驗邊界處分類結果；最后基于Google高清影像對所有驗證點位進行目視解譯，分配類別屬性，計算混淆矩陣評估準確性。2）選取4個典型區域，基于亞米級Google高清影像，通過目視解譯繪制水產養殖池塘實際水面，逐一比對本研究方法識別結果。

圖9 為各驗證點位的空間分布，表1為養殖池塘與非養殖池塘混淆矩陣，檢驗結果表明本研究算法具有較高準確性，總體精度達到0.921，Kappa系數為0.842。其中養殖池塘的生產者精度為0.908，有48個池塘點位被遺漏未被識別出；非養殖池塘的生產者精度為0.936，有31個點位被誤分為池塘，分類錯誤多出現于池塘水面的邊界。同時，表2統計4個典型區域中本研究水產養殖池塘識別結果占Google目視解譯結果的面積比例，河流入海口處散亂分布式的養殖池塘（A）、規模化集聚型養殖池塘（B）（D）的面積占比均高于90%；廉州灣沿岸養殖池塘（C）面積小且更為密集，水面邊界像元混合現象更嚴重，導致面積占比較小為80.76%，這也是該區域錯分點位較多的原因。

表1 養殖池塘與非養殖池塘混淆矩陣Table 1 Confusion matrices for two classes of aquaculture ponds and non-aquaculture ponds

表2 水產養殖池塘識別結果Table 2 Identification result of aquacwture ponds

3.2 本研究的潛力以及局限性

部分研究基于Landsat-8影像識別了廣西北部灣沿岸水產養殖池塘，2015—2017年面積依次為389.04 km2[34]、337 km2[19]和346.63 km2[14]，但由于空間分辨率限制Landsat-8影像無法剔除池塘間堤防，同時還可能包括提取區域內的引水渠、細小河流等其他水體，識別結果為養殖區域的用地總面積，多用于評估土地利用覆蓋類型，相較于養殖實際水面面積明顯偏大。在本研究中，將池塘間堤防也作為分類對象，基于10 m空間分辨率的Sentinel-1、Sentinel-2數據提取養殖實際水面，最大限度地將每個池塘識別為獨立目標，在評估水產養殖池塘的空間分布、面積上具有更大優勢。此外，本研究采用時序數據代替原始數據作為分類特征值，IF能夠反映一年之中像元尺度下的水淹狀態，消除廢棄池塘、水稻田和季節性水域的影響；NDWImean(85%～95%)增大了水體與非水體像元的可分離性，能夠更好識別養殖實際水面的邊界；不透水面在年內一般無重大變化，采用VHYEARmean、SWIR2mean(10%～90%)年度平均值可以減少斑點噪聲影響（Sentinel-2影像難以完全去除含云像元，仍有部分異常值，因此，適當剔除了部分兩端極值）。將本研究方法應用于廣西北部灣海岸帶中效果顯著，結果表明時序數據在水產養殖池塘識別上具有顯著優勢，同時基于GEE平臺的強大性能，該方法可以輕松用于其他地區的水產養殖池塘識別。

但本研究仍存在以下不足：1）廉州灣沿岸，養殖池塘面積普遍偏小且排布密集，大部分池塘間堤防寬度在5～10 m，由于影像分辨率限制，池塘邊緣部分往往因像元混合被識別為非水像元，導致該區域提取面積較實際偏小（圖10c），這也是主要誤差來源；2）小于5 m的細小堤防光譜特征不明顯往往被忽略，因分割困難會被保留與池塘水面合并；3）部分青蟹養殖池塘因長時間被水生植物所覆蓋，在遙感影像中水體特征并不明顯，難以用水體識別方法提取。

3.3 對于推進水產養殖業綠色發展的意義

廣西北部灣水產養殖池塘大多沿海岸、河流散亂分布，如養殖尾水未經處理直接排入近海會嚴重危害水環境生態健康。近年來廣西近岸海域發生赤潮的頻率呈增加趨勢，造成了較大的經濟損失和海洋環境破壞[35]。水產養殖的水體環境質量和養殖尾水排放污染防治，已成為水產養殖業可持續發展的制約因素。廣西壯族自治區農業農村廳于2019年6月提出《關于加快推進廣西水產養殖業綠色發展的實施意見》，到2022年實現水產養殖業綠色發展空間布局明顯優化，到2035年實現養殖尾水全面達標排放。本研究方法最大限度地將池塘識別為獨立目標，可用于評估區域養殖池塘聚集程度，結合衛星影像監測尾水排放，為優化養殖空間布局、規劃尾水處理設施以及養殖區域的智能監控提供有效數據支撐，對于推進水產養殖業綠色發展具有重要意義。此外，基于水產養殖池塘的空間分布特征，欽州灣、鐵山港將是攻堅難點，池塘沿整個海灣岸線分布過于散亂，要達到尾水處理設施全覆蓋較為困難，需盡早優化養殖區域的空間布局，推動養殖工廠化、集中連片化，進而建設高標準、高效率的尾水處理設施。

4 結 論

在提取大范圍復雜環境中的沿海水產養殖池塘時，針對傳統遙感識別方法中存在分類精度不高的問題，本研究提出一種基于GEE平臺和時序遙感數據，結合多閾值分割以及面向對象分類的水產養殖池塘識別方法，并將其應用于廣西北部灣海岸帶，得出如下結論：

1）廣西北部灣海岸帶水產養殖池塘面積共計199.3 km2，北海養殖面積最大為112.9 km2，防城港和欽州分別為44.2、42.2 km2。在廉州灣，沿岸養殖池塘分布最為密集；在欽州灣和鐵山港，養殖池塘沿整個海灣岸線分布較于散亂，要達到尾水處理設施全覆蓋較為困難，需盡早優化養殖區域的空間布局。

2）水淹頻率反映了像元尺度下全年的水淹狀況，能夠有效排除廢棄池塘、水稻田和季節性水域；NDWImean(85～95%)為時間序列中按像元NDWI數值升序排列的85%～90%區間的平均值，可以增強水體與非水體的可分離性，能更好地區分池塘間的堤防。相較于以單一日期影像計算分類特征值的方法，時序遙感數據在識別養殖池塘上具有顯著優勢。

3）本研究方法總體精度達到0.921，Kappa系數為0.842，提取結果更接近養殖實際水面面積，在大范圍復雜環境中仍具有較高準確性。同時，憑借GEE平臺強大性能，該方法可用于識別其他地區的水產養殖池塘，對于科學設置區域養殖發展布局，制定環境保護措施，推動水產養殖業綠色發展具有重大意義。