随着社会经济的快速发展，我国的水产渔业也稳步向前, 逐渐从过度消耗天然渔业到人工放养再到围网养殖^[1].近20年来，长江流域的浅水湖泊围网养殖业快速扩张，极大地增加了经济效益，但超高密度的围网养殖也对湖泊生态造成了严重破坏.一方面，围网养殖区中的围网会阻隔水流，减缓污染物的扩散速度，且投放的饵料增加氮、磷等浓度，加剧了浅水湖泊的内源污染^[2]；另一方面，养殖的鱼蟹消耗甚至破坏了水生植物，降低湖泊的自净能力，打破了湖泊的生态平衡^[3].在此背景下，合理科学地规划并控制湖泊围网已被湖泊管理部门列为湖泊水质提升和优化的重点工作，而获取并掌握湖泊围网养殖的时空分布信息是湖泊生态管理部门在科学管理围网养殖区时制定措施的依据.

传统的围网养殖区空间分布及其面积的获取方法即人工实地调查测量法，该方法时间周期长，人力投入大.遥感技术由于其具有实时、大范围、动态和客观等优势，而逐渐成为湖泊围网养殖区监测和湖泊水环境调查的主流方法.沈芳等^[4]采用人工目视解译遥感影像后数字化的方法，提取了滆湖和洮湖围网养殖区.人工目视解译虽具有较高的精度，但仍然费时费力，很难批量自动提取.为克服该方法无法批量提取的缺陷，近年来，学者们先后基于各种遥感数据探讨了围网区的自动提取方法.如，马艳娟等^[5]基于ASTER数据直接利用光谱信息提取了烟台附近的近海围网养殖区; 朱长明等^[6]基于高分影像通过纹理信息提取了连云港附近的近海围网养殖区；卢业伟等^[7]基于RAPIDEYE数据和纹理信息提取了澳门附近的近海围网养殖区.但针对浅水湖泊的围网养殖区自动提取方法的研究相对较少，李俊杰等^[8]基于资源卫星2号初步探讨了如何运用纹理信息提取围网养殖区; 王静等^[9]基于Landsat ETM卫星数据通过对应分析的方法尝试提取围网养殖区.以上研究主要集中于使用影像的纹理信息提取围网养殖区，但图像的大多纹理特征计算复杂，运算量大，且不同时相的影像往往需要使用不同的纹理特征.梯度变换作为图像处理和模式识别的常用方法，自提出就受到了广泛运用，在此基础上提出的图像处理算法层出不穷.袁春兰等^[10]使用Sobel算子对用过Otsu算法运算后的红外图像进行梯度变换成功检测到该图像的边缘；何春华等^[11]使用自适应权值的Sobel算子进行图像的梯度变换检测图像边缘并获得较好结果；李灿灿等^[12]使用多方向Sobel算子对叶脉图像进行梯度变换以检测其边缘变化且得到很高的检测精度.以上基于不同Sobel算子的边缘检测研究都表明了梯度变换可以有效地提取图像边缘，增强图像细节^[13-14].

本文以长江下游典型的围网养殖湖泊——江苏苏州阳澄湖为研究区，以资源三号(ZY-3)高分影像为数据源，利用梯度变换尝试提出一种普适性和准确率较高的围网养殖区遥感提取算法，提取浅水湖泊围网养殖区，也可为湖泊环境的政府部门制定湖泊水质提升和围网区合理规划政策提供支持.

1 研究区概况与数据获取 1.1 研究区概况

阳澄湖(31°21′~31°30′N, 120°39′~120°51′E)地处江苏苏州吴中区与昆山市之间，是古太湖的残留，南北长17 km，东西最大宽度8 km，面积119.04 km²，蓄水量1.67×10⁸ m³，平均水深1.4 m^[15].阳澄湖平均水深适中，利于鱼蟹存活，水产资源十分丰富，湖中六宝远近闻名，其中大闸蟹更是驰名中外，也是较早开展围网养殖的浅水湖泊之一.自1992年开始，阳澄湖进行生产性围网养殖，在21世纪初形成围网养殖高峰，其中主要以大闸蟹养殖为主.

1.2 数据搜集

本文选用的遥感数据源为国产的高分辨率多光谱卫星影像数据——资源三号(ZY-3)，影像获取时间是2015年4月25日，影像多光谱数据空间分辨率为5.8 m，数据记录格式为GEOTIFF，产品总共有4个波段，分别是近红外波段(0.77~0.89 μm)、红波段(0.63~0.69 μm)、绿波段(0.52~0.59 μm)和蓝波段(0.45~0.52 μm).结合阳澄湖的矢量边界对影像进行裁切(图 1)，剪切后的研究区影像共计3147*2736个像元.考虑到背景在各波段的灰度值均为0，因此将各波段均为0的像元值赋值为NAN.

2 研究方法与提取结果 2.1 方法概述

典型围网区中围网和自然水体在ZY-3影像各波段的灰度值可以看出(图 2)，围网和自然水体对光谱吸收有差异，波长最长的近红外波段在围网区的灰度值变化率最大.因此，可以通过计算近红外波段值的空间变化率来区分围网和自然水体，而计算变化率在数学上最常用的工具是梯度.本文尝试以梯度计算为核心提出一种新算法提取围网区，具体流程如图 3.

2.2 梯度变换及优化

梯度是描述函数变化的矢量，一般用∇表示，其模值是函数在最大变化方向上的变化率大小.二维平面上，函数φ的梯度可用列向量$ \nabla \varphi = \left[\begin{array}{l} \frac{{\partial \varphi \left( {x, y} \right)}}{{\partial x}}\\ \frac{{\partial \varphi \left( {x, y} \right)}}{{\partial y}} \end{array} \right]$表示，模值为：

$ \left| {\nabla \varphi } \right| = {\left[{{{\left( {\frac{{\partial \varphi \left( {x, y} \right)}}{{\partial x}}} \right)}^2} + {{\left( {\frac{{\partial \varphi \left( {x, y} \right)}}{{\partial y}}} \right)}^2}} \right]^{\frac{1}{2}}} $

(1)

式中，x表示该点相对于参考点的横向距离，y表示该点相对于参考点的纵向距离.

实际应用中，为简化变化率，即梯度模值的计算，认为在x、y变化较小时，可以用绝对值和近似平方和的平方根表示^[16]，式(1)也可以表示为：

$ \left| {\nabla \varphi \left( {x, y} \right)} \right| \approx \left| {\frac{{\partial \varphi \left( {x, y} \right)}}{{\partial x}}} \right| + \left| {\frac{{\partial \varphi \left( {x, y} \right)}}{{\partial y}}} \right| $

(2)

由于遥感影像的DN值在空域上是离散的，计算它的梯度模值就需将公式(2)离散化.经过研究，Roberts算子、Sobel算子和Prewitt算子相继被提出并用来计算离散化的梯度模值^[17].考虑到Sobel算子计算后的结果相对于Roberts算子和Prewitt算子更精确及平滑^[18]，本文使用Sobel算子.根据算子^[10]：

$ \begin{array}{l} \frac{{\partial \varphi \left( {i, j} \right)}}{{\partial x}} = \left( {\varphi \left( {i-1, j + 1} \right) + 2\varphi \left( {i, j + 1} \right) + \varphi \left( {i + 1, j + 1} \right)} \right)-(\varphi \left( {i-1, j + 1} \right) + \\ 2\varphi \left( {i, j - 1} \right) + \varphi \left( {i + 1, j - 1} \right))\frac{{\partial \varphi \left( {i, j} \right)}}{{\partial y}} = (\varphi \left( {i + 1, j - 1} \right) + 2\varphi \left( {i + 1, j} \right) + \\ \varphi \left( {i + 1, j + 1} \right)) - \left( {\varphi \left( {i - 1, j - 1} \right) + 2\varphi \left( {i - 1, j} \right) + \varphi \left( {i - 1, j + 1} \right)} \right) \end{array} $

(3)

将式(3)代入式(2)可以得到：

$ \begin{array}{l} \left| {\nabla {\varphi ^n}\left( {i, j} \right)} \right| \cong |\left[{{\varphi ^n}\left( {i + 1, j-1} \right) + 2{\varphi ^n}\left( {i + 1, j} \right) + {\varphi ^n}\left( {i + 1, j + 1} \right)} \right] - [{\varphi ^n}\left( {i-1, j-1} \right) + \\ 2{\varphi ^n}\left( {i-1, j} \right) + {\varphi ^n}\left( {i - 1, j + 1} \right)]| + |[{\varphi ^n}\left( {i-1, j + 1} \right) + 2{\varphi ^n}\left( {i, j + 1} \right) + \\ {\varphi ^n}\left( {i + 1, j + 1} \right)] - \left[{{\varphi ^n}\left( {i-1, j-1} \right) + 2{\varphi ^n}\left( {i, j-1} \right) + {\varphi ^n}\left( {i + 1, j - 1} \right)} \right]| \end{array} $

(4)

式中, n表示波段；φ表示影像波段的DN值；i, j分别为影像的行列数.

当n=4时，变化率可较好地识别出围网和自然水体(图 2)，因此，通过对第4波段进行以上梯度变换计算，其结果如图 4.

统计梯度变换后的影像灰度值，由统计结果可知0.001~0.150区间的灰度值占比较高，为增强特征利于分类，选用幂率函数中的根变换进行灰度变换^[17].此外，自然情况下湖面浪花的镜面反射给梯度变换后的新影像引入“椒盐”噪声(图 4).为了提高围网提取的可靠性，本文采用中值滤波器去除椒盐噪声^[19]，且中值滤波器的模板大小根据噪声的空间密度而定^[20].

2.3 围网区提取及结果优化

在进行灰度变换和去噪后，进行围网区阈值的确定，具体方法为：参考目视解译的结果，在不同的区域分别选取具有代表性和典型性的围网样本作为训练样本，共选取20个样本，并保证每个样本至少包含10个围网网格，然后计算样本的数学期望将其作为围网提取阈值，为0.8.通过该阈值将围网区全部提取，结果如图 5A.但由于水生植被与水体在近红外的灰度值差异也较大，因此，此时提取的围网区还混有部分水生植被区.通过比较水生植被区和围网区的影像特征，发现围网区区域内有孔洞，且空间上孔洞分布存在周期性，但大小随空间分布没有规律(图 5A).为了剔除水生植被区，找出孔洞是关键，而标准大小的孔洞比随意大小的孔洞更易识别.根据围网的形状特征的统计结果，发现围网孔洞大小基本介于4*4像元和15*20像元之间，且频率最高的孔洞大小为12*18像元.为此，本算法先逐行扫描围网区寻找间隔在4~15像元之间的围网网格，并将这些网格通过填充(像元灰度值设置为1)或挖去(像元灰度值设置为0)的方式变成横向间隔12像元的格子，然后再逐列扫描寻找间隔在4~20像元之间的围网网格，并以同样的方式将它们变成纵向间隔18像元的格子.这样就可以将其规格范围内的孔洞全部转变为12*18像元大小(图 5B).通过对围网区的像元个数进行统计，可知阳澄湖围网区宽度基本都大于180 m即3个孔洞周期，因此本算法将筛选出的孔洞周期数大于3的区域作为围网区，结果如图 5C.最后通过将空洞处灰度值设置为1，将提取区域的孔洞填满(图 5D).

3 精度评价

采用人工解译的围网区作为参考标准，对本论文算法提取的围网区进行精度验证.具体方法为：将算法提取的围网区和人工解译的围网区转变成矢量图层，假设算法提取的围网区为集合A，人工解译的围网区为集合B，则正确分类的围网区C=A∩B，漏分的围网区O=B-C，错分的围网区E=A-C(图 6)，分别计算A、B、C、O、E的面积，并根据以下公式计算漏分误差、错分精度和总体精度^[21]：

$ {p_{\rm{o}}} = \frac{{{S_{\rm{O}}}}}{{{S_{\rm{B}}}}} \times 100\% $

(5)

$ {p_{\rm{e}}} = \frac{{{S_{\rm{E}}}}}{{{S_{\rm{A}}}}} \times 100\% $

(6)

$ {p_{\rm{e}}} = \frac{{{S_{\rm{C}}} + {S_{\rm{W}}}}}{S} = \frac{{{S_{\rm{C}}} + S-{S_{\rm{A}}}-{S_{\rm{O}}}}}{S} \times 100\% $

(7)

式中，p_o为漏分误差，p_e为错分误差，p_c为总体精度，S_O为漏分的面积，S_E为错分的围网面积，S_C为正确分类的围网面积，S_W为正确分类的水体面积，S_A为算法提取的围网面积，S_B为人工解译的围网面积，S为湖泊总面积.

基于以上验证方法得到的精度评价结果如表 1.该算法对围网区的提取精度为90.66%，错分精度为17.95%，漏分误差为15.35%.

4 结论与讨论

基于ZY-3卫星影像，根据围网区与非围网区在影像上的光谱空间变化特征，利用梯度变换提出一种围网区遥感识别方法，并利用该方法对2015年阳澄湖的围网区进行提取，并以人工目视解译结果作为标准，对提取结果进行验证和精度评价.结果表明：该方法对围网区的提取精度达到90.66%，基本满足实际应用需求.

本文以阳澄湖为研究区提出了围网提取方法，且取得较高的识别精度.理论上，该算法具有普适性，可用于其他浅水湖泊的围网提取，但考虑到不同湖泊的水环境因素差异较大，该算法能否准确提取其它的浅水湖泊(如滆湖、长荡湖等)的围网区还有待进一步验证.此外，该算法是基于高分影像ZY-3为数据源发展的，理论上该算法可以用于任何卫星影像提取浅水湖泊围网，但受空间和光谱分辨率的限制，该算法是否适用于其他卫星影像(如更低分辨率的Landsat和HJ-CCD等)提取浅水湖泊围网区，精度如何等还有待进一步研究.同时，该算法依然存在误差，理论上很多因素均会造成精度损失，如样本的选取方式、样本数量的确定和阈值的计算等，如何分析不同因素对精度的影响进而提高提取精度还有待进一步研究.