基于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平台强大性能，该方法可用于识别其他地区的水产养殖池塘，对于科学设置区域养殖发展布局，制定环境保护措施，推动水产养殖业绿色发展具有重大意义。