一种基于zy-3影像的浅水湖泊围网区提取方法
【专利摘要】本发明提出一种基于ZY?3影像的湖泊围网区快速提取方法，包括：获取研究区ZY?3影像，并进行影像预处理；根据预处理后的影像，选取典型的围网样本库；对预处理后的影像进行梯度变换；基于变换后的影像，计算围网样区的均值和方差；根据围网样区的均值和方差确定围网区阈值，并初步提取围网区；优化围网区提取结果，包括去除“椒盐”噪声、低频成分，以及利用膨胀、腐蚀去除碎斑；计算围网区的几何特征，以解译的围网区的分布及面积作为参考，对提取结果进行精度验证。利用本发明的提取方法可快速、实时、大范围和高精度的提取围网，为进一步开展浅水湖泊围网面积与水环境的关系研究提供数据基础，为水质优化提升方案制定提供依据。
【专利说明】
一种基于ZY-3影像的浅水湖泊围网区提取方法
技术领域
[0001] 本发明属于遥感应用领域，涉及一种浅水湖泊围网区的遥感提取方法。
【背景技术】
[0002] 浅水湖泊围网养殖可以取得较高的经济效益，但高密度的浅水湖泊围网养殖会加 剧浅水湖泊的富营养化进程。因此，合理科学的规划并控制湖泊围网已被湖泊管理部门列 为湖泊水质提升和优化的重点工作，而获取并掌握湖泊围网养殖面积及时空分布信息对合 理规划和控制围网具有重要的指导和现实意义。
[0003] 目前围网的提取方法仍是用人工统计、遥感影像人工解译法，费时费力且无法实 时、动态的获取围网区分布及面积信息。因此，提出一种基于遥感影像的浅水湖泊围网养殖 区快速提取方法，从而快速、准确的获取浅水湖泊围网养殖区分布及面积信息势在必行。
【发明内容】
[0004] 针对现有技术的缺陷与不足，本发明旨在提出一种基于遥感影像的浅水湖泊围网 养殖区快速提取方法，可快速、准确的获取浅水湖泊围网养殖区分布及面积信息。
[0005] 为达成上述目的，本发明的第一方面提出一种基于ZY-3影像的湖泊围网区快速提 取方法，其包括以下步骤：
[0006] 步骤一、获取研究区ZY-3影像，并进行预处理；
[0007] 步骤二、根据预处理后的影像，选取典型的围网样本库；
[0008] 步骤三、对预处理后的影像进行梯度变换；
[0009] 步骤四、基于变换后的影像，计算围网样区的均值和方差；
[0010] 步骤五、根据围网样区的均值和方差确定围网区阈值，并初步提取围网区；
[0011] 步骤六、优化围网区提取结果，包括利用中值滤波去除"椒盐"噪声，利用拉普拉斯 算子变换去除低频成分，利用膨胀、腐蚀去除碎斑；
[0012] 步骤七、计算围网区的几何特征，即面积和周长，以解译的围网区的分布及面积作 为参考，对提取结果进行精度验证。
[0013] 由以上本发明的技术方案可知，利用本发明的提取方法可快速、实时、大范围和高 精度的提取围网，为进一步开展浅水湖泊围网面积与水环境的关系研究提供数据基础，为 水质优化提升方案制定提供依据。
[0014] 应当理解，前述构思以及在下面更加详细地描述的额外构思的所有组合只要在这 样的构思不相互矛盾的情况下都可以被视为本公开的发明主题的一部分。另外，所要求保 护的主题的所有组合都被视为本公开的发明主题的一部分。
[0015] 结合附图从下面的描述中可以更加全面地理解本发明教导的前述和其他方面、实 施例和特征。本发明的其他附加方面例如示例性实施方式的特征和/或有益效果将在下面 的描述中显见，或通过根据本发明教导的【具体实施方式】的实践中得知。
【附图说明】
[0016] 附图不意在按比例绘制。在附图中，在各个图中示出的每个相同或近似相同的组 成部分可以用相同的标号表示。为了清晰起见，在每个图中，并非每个组成部分均被标记。 现在，将通过例子并参考附图来描述本发明的各个方面的实施例，其中：
[0017] 图1为本发明的基于ZY-3影像的浅水湖泊围网区提取方法的流程图。
[0018] 图2为本发明浅水湖泊围网区提取方法实施例中基于算法提取的2015年4月25日 阳澄湖的围网养殖区分布图。
[0019] 图3为本发明浅水湖泊围网区提取方法实施例中2015年4月25日的围网养殖区精 度评价结果图。
[0020] 具体的实施方式
[0021] 在本公开中参照附图来描述本发明的各方面，附图中示出了许多说明的实施例。 本公开的实施例不必定意在包括本发明的所有方面。应当理解，上面介绍的多种构思和实 施例，以及下面更加详细地描述的那些构思和实施方式可以以很多方式中任意一种来实 施，这是因为本发明所公开的构思和实施例并不限于任何实施方式。另外，本发明公开的一 些方面可以单独使用，或者与本发明公开的其他方面的任何适当组合来使用。
[0022] 结合图1所示，根据本发明的实施例，一种基于ZY_3(即资源三号卫星)影像的湖泊 围网区快速提取方法，其包括以下步骤：
[0023] 步骤一、获取研究区ZY-3影像，并进行影像预处理；
[0024]步骤二、根据预处理后的影像，选取典型的围网样本库；
[0025]步骤三、对预处理后的影像进行梯度变换；
[0026] 步骤四、基于变换后的影像，计算围网样区的均值和方差；
[0027] 步骤五、根据围网样区的均值和方差确定围网区阈值，并初步提取围网区；
[0028] 步骤六、优化围网区提取结果，包括利用中值滤波去除"椒盐"噪声，利用拉普拉斯 算子变换去除低频成分，利用膨胀、腐蚀去除碎斑；
[0029] 步骤七、计算围网的几何特征，即面积和周长，以解译结果作为参考，对提取结果 进行精度验证。
[0030] 在一些实施例中，所述步骤一的影像预处理可以采用现有的一些算法实现。本实 施例中，影像预处理包括辐射校正、几何校正和研究区裁剪，其中所述辐射校正指辐射定标 和大气校正。
[0031 ]作为可选的方式，所述步骤三中对影像的梯度变换的具体方法如下：
[0032] 梯度计算公式为：
，其模值可以表示为
[0033] 考虑到结果精确和曲线平滑这两项要求，本实施例中，选取Sobel算子进行梯度变 换，其计算公式为： i?(pai s i \cpna +1,/ -1) + 2 * (pnu +1 j) + g)n (^i 1,j + 1)]
[0034] - [(pn(i - l,j - 1) + 2 * r/>n(i - l,j) + (f)n(i - 1J + 1)]| + | [(pn{i 1,/ + 1) + 2 * + 1) + (pn{i + l,j + 1)] -[(pn(i - l,j - 1) + 2 * (pn{i,j - 1) + (pn(i + l,j - 1)]|
[0035] 其中，炉表示该图像的某一像元值，i和j分别表示该像元所在的行数和列数，n为 波段数。
[0036] 步骤四中，所述的围网样区的均值和方差，计算方法如下：
[0037] (1)平均值：以矩阵的形式来表示样区影像变换后的数据，则设获取的围网样区的 遥感影像的变换数据矩阵为X，如下：
[0039]其中，n为影像的波段数，p、q分别为样区影像的行列数。则围网样区平均值为：
[0040] (2)方差:基于围网样区平均值和数据矩阵X，通过矩阵￥"=0"-&￥^")2的平均值 计算方差，其中围网样区平均值的计算方法如（1)中所示，即
[0042] 其中n为影像的波段数，p、q分别为样区影像的行列数，〇为方差，y为矩阵Y的元素， avrl为待选取的阈值。
[0043] 其中，所述的步骤五中的围网区阈值选取为围网样本区变换后的平均值±15区间 内方差最小的值。
[0044] 其中，所述步骤六围网提取结果优化方法及操作如下：
[0045] (1)确定两个矩阵模板:6*6的A和10*10的B;
[0046] (2)使用闭操作去除拉普拉斯算子变换后留下的小孔洞，其中闭操作是先用A膨胀 图像，再用A腐蚀膨胀后的图像；
[0047] (3)使用B模板再次腐蚀图像去除碎斑。
[0048] 其中，拉普拉斯算子变换的具体操作如下：
[0049] 根据拉普拉斯算子中心差分格式，将研究区的每个波段都进行一次离散拉普拉斯 算子变换。其中，拉普拉斯算子中心差分计算公式为：
[0051] 其中，识表示该图像的某一像元值，i和j分别表示该像元所在的行数和列数，A y、 Ax分别表示纵向及横向距离，n为波段数。在变换实践过程中，一般选用3*3模板计算，则其 计算公式如下： ￥ ^V(pfj = 4-1) + (pn{i, j - 1) + (pv〇 + !,/) + (pn{i - l，i)
[0052] + cpn(i - \,j + 1) + ^ 1) + (pn(^i + l,j + l) + <pn{i - l,j - 1) - 8 * (pn(i,j)]
[0053] 其中，识表示该图像的某一像元值，i、j分别表示该像元所在的行数和列数，n为波 段数。由此，将图像矩阵中的每一个像元带入上式计算获得变换值。
[0054]腐蚀是指拖动模板在图像域移动，横向移动间隔取1个像元，纵向移动间隔取1个 扫描行。在每一个位置上，当模板的中心点平移到图像上的某一点，如果模板内的每一个像 元都与以该点为中心的相同邻域中对应像元完全相同，则保留该像元点，对于不满足条件 的像元点则全部删除。
[0055] 膨胀图像是指在每一个位置上，当模板的中心点平移到图像上的某一点，如果模 板内的有像元都与以该点为中心的相同邻域中对应像元重合那么就保留该像元点。
[0056] 其中，所述步骤七围网区域的几何特征，即面积和周长的计算和精度评价的方法 如下：
[0057] (1)几何特征计算:计算之前应先将围网区的网格状形态去除，具体操作如下：
[0058] 1)将所有值为1且距离小于围网最大长度的像元连接起来；
[0059] 2)使用最大值滤波去除空洞；
[0060] 3)使用中值滤波去除噪声；
[0061] 在去除了网格状形态之后根据围网区像元及围网区边界的像元的数量统计围网 区的几何特征。
[0062] (2)精度评价:对提取的围网区的精度进行评价前需要计算错分漏分面积。
[0063]错分漏分面积具体操作如下：1)通过Arcgis软件将算法提取的和人工解译的围网 区转变成矢量图层，假设算法提取的为A，人工解译的为B，则正确分类的围网区C = AnB，漏 分的围网区〇 = B-C，错分的围网区E=A-C，分别计算A、B、C、0、E的面积。
[0064] 2)计算误差和精度，计算公式如下：
[0065] p〇 = S(0)/S(B)；
[0066] pe = S(E)/S(A)；
[0067] pc = S(C)+S-S(A)-S(0)/S；
[0068] 其中p。是漏分误差，P(5是错分误差，p。是总体精度，S(0)是漏分的面积，S(E)是错分 的面积，S(C)是正确分类的面积，S(A)是算法提取的面积，S(B)是人工解译的面积，S是湖泊 总面积。
[0069] 下面结合附图2和附图3本发明进行进一步说明。本例中，以中国江苏典型的使用 围网的湖泊一一阳澄湖为例进行实例分析。
[0070] 阳澄湖（31° 2V~31° 3(T N、120° 3V~120° 5V E)地处江苏苏州吴中区与昆山市之 间，是古太湖的残留，南北长17千米，东西最大宽度8千米，面积119.04平方千米，蓄水量 1.67亿立方米，平均水深1.4米。经调研和资料显示，目前阳澄湖中湖和东湖有围网。下面以 阳澄湖为研究区，进行实例分析：
[0071] 步骤一、根据实测样点获取时间，从中国资源卫星应用中心的数据产品查询订购 网站(http: //218 ? 247 ? 138 ? 121/1)33?1&七;1^〇1'111/;[11(161.111：1111)订购下载2￥-3卫星影像。
[0072]步骤二、利用ENVI软件分别对影像进行大气校正(包括辐射定标和大气校正)和几 何校正，然后利用阳澄湖矢量边界裁剪出的研究区影像；
[0073]步骤三、利用ENVI软件选取典型的围网样本库；
[0074]步骤四、对影像进行变换，计算梯度变换值提取围网区；
[0075]步骤五、计算样本的平均值和方差，确定阈值；
[0076]步骤六、利用阈值初步提取围网养殖区；
[0077]步骤七、利用拉普拉斯算子变换和形态学处理方法优化提取结果；
[0078]步骤八、计算围网的几何特征(面积和周长），以人工解译结果作为参考，对提取结 果进行精度验证。
[0079]结合图2、图3所示，图2为本发明所述的方法实施例中基于2015年4月25日影像提 取的围网养殖区分布结果；图3为本发明所述的方法实施例中，以人工解译的围网养殖区为 参考结果，对本发明方法提取的围网养殖区的精度比较结果图。
[0080] 根据本发明的前述实施例的提取方法进行围网区的提取和进度评价，得到如下表 1的评价表。
[0081] 表1本发明实施例所述方法获取的阳澄湖围网养殖区提取结果的精度评价表
[0083] 表1是本发明实施例的提取方法获取的阳澄湖围网养殖区提取结果的精度评价 表。通过监测结果与实测结果的对比可知，本发明方法对阳澄湖围网养殖区的提取精度高 于90.66%，符合实际应用需求。
[0084] 尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解:在不 脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本 发明的范围由权利要求及其等同物限定。
【主权项】
1. 一种基于ZY-3影像的湖泊围网区快速提取方法，其特征在于，该方法包括W下步骤： 步骤一、获取研究区ZY-3影像，并进行影像预处理； 步骤二、根据预处理后的影像，选取典型的围网样本库； 步骤=、对预处理后的影像进行梯度变换； 步骤四、基于变换后的影像，计算围网样区的均值和方差； 步骤五、根据围网样区的均值和方差确定围网区阔值，并初步提取围网区； 步骤六、优化围网区提取结果，包括利用中值滤波去除"椒盐"噪声，利用拉普拉斯算子 变换去除低频成分，W及利用膨胀、腐蚀去除碎斑； 步骤屯、计算围网区的几何特征，W解译的围网区的分布及面积作为参考，对提取结果 进行精度验证。2. 根据权利要求1所述的方法，其特征在于，所述步骤一中的影像预处理包括福射校 正、几何校正和研究区裁剪，其中所述福射校正指福射定标和大气校正。3. 根据权利要求1所述的方法，所述的对影像的梯度变换的具体方法如下： 选取Sobel算子进行梯度变换，其计算公式为：其中，取表示图像的某一像元值，i和j分别表示像元所在的行数和列数，n为波段数。4. 根据权利要求1所述的方法，其特征在于，所述步骤四中的计算围网样区均值和方 差，计算方法如下： (1) 平均值计算：W矩阵的形式表示围网样区影像变换后的数据，则设获取的围网样区 的变换影像数据矩阵为X，如下：其中，n为影像的波段数，P和q分别为样区影像的行列数，则围网样区平均值为：(2) 方差计算:基于平均值和数据矩阵X，通过矩阵YD= (xn-avrr)2的平均值计算方差， 其中平均值的计算方法如（1)中所元.目口其中，n为影像的波段数，P和q分别为样区影像的行列数，O为方差，y为矩阵Y的元素， avr 1为待选取的阔值。5. 根据权利要求4所述的方法，其特征在于，所述步骤五中围网区阔值的选取为步骤四 中建立的围网样本区变换后的平均值±15区间内方差最小的值。6. 根据权利要求1所述的方法，其特征在于，所述步骤六的围网提取结果优化方法的操 作具体如下： (1) 确定两个矩阵模板:6*6的A和10*10的B; (2) 使用闭操作去除拉普拉斯算子变换后留下的小孔桐，其中闭操作是先用A膨胀图 像，再用A腐蚀膨胀后的图像； (3) 使用財莫板再次腐蚀图像去除碎斑； 其中，拉普拉斯算子变换的具体操作如下： 根据拉普拉斯算子中屯、差分格式，将研究区的每个波段都进行一次离散拉普拉斯算子 变换，其中拉普拉斯算子中屯、差分计算公式为：其甲，沪巧不图像W呆一像兀但，1卿J分別巧不像兀所仕W打甄卿列甄，A y、A X分别表 示纵向及横向距离，n为波段数； ^巧抢讨賴由.佛田糾3摇版if貸.则if貸公立加下.共中，沪巧不阁像的呆一像兀但，1、J分別巧不像兀所化的竹数卿列数，n刃做段数，由 此，将图像矩阵中的每一个像元带入上式计算获得变换值； 腐蚀表示拖动模板在图像域移动，横向移动间隔取1个像元，纵向移动间隔取1个扫描 行;在每一个位置上，当模板的中屯、点平移到图像上的某一点，如果模板内的每一个像元都 与W该点为中屯、的相同邻域中对应像元完全相同，则保留该像元点，对于不满足条件的像 元点则全部删除;膨胀图像是指在每一个位置上，当模板的中屯、点平移到图像上的某一点， 如果模板内的有像元都与W该点为中屯、的相同邻域中对应像元重合则就保留该像元点。7. 根据权利要求1所述的方法，其特征在于，所述步骤屯的围网区的几何特征的计算和 精度评价的方法如下： (1)几何特征计算:在计算几何特征之前先将围网区的网格状形态去除，具体去除操作 如下：1)将所有值为1且距离小于围网最大长度的像元连接起来;2)使用最大值滤波去除空 桐；3)使用中值滤波去除噪声;然后根据围网区像元及围网区边界的像元的数量统计围网 区的几何特征； (2)精度评价:对提取的围网养殖区的精度进行评价前需要计算错分漏分面积，具体操 作如下：1)通过Arcgis软件将前述步骤六优化后的围网区和解译的围网区转变成矢量图 层，假设算法提取的围网区为集合A，人工解译的为集合B，则正确分类的围网区C = AnB,漏 分的围网区O = B-C,错分的围网区E = A-C,分别计算A、B、C、0、E的面积;2)计算误差和精度， 计算公式如下： Po = S(OVSW); Pe = S(E)ZS(A); Pc = S(C)+S-S(A)-S(0)/S； 其中Po是漏分误差，Pe是错分误差，Pc是总体精度，S(O)是漏分的面积，S巧)是错分的面 积，S(C)是正确分类的面积，S(A)是算法提取的面积，S(B)是人工解译的面积，S是湖泊总面 积。
【文档编号】G06T7/00GK105913426SQ201610222606
【公开日】2016年8月31日
【申请日】2016年4月11日
【发明人】罗菊花, 马荣华, 黄帅, 李飞
【申请人】中国科学院南京地理与湖泊研究所