CN113129320B - 一种从航空正射影像中提取养殖水塘边界的方法和系统 - Google Patents

Abstract

本发明公开了一种从航空正射影像中提取养殖水塘边界的方法，包括：1、建立候选检测点集合；2、从候选检测点集合S_c‑dtc中筛选检测点；3、根据检测点集合S_dtc提取航空DOM影像中的水塘像素；4、提取水塘像素的边界；5、提取水塘阴影；水塘阴影与水塘像素构成航空DOM影像中水塘的全部像素集合；6、对整幅影像计算水塘边界并简化。该方法能够减少阴影和其它地物对水塘边界提取的影响，提取到精确的养殖水塘边界。

Description

一种从航空正射影像中提取养殖水塘边界的方法和系统

技术领域

本发明属于航空测绘技术领域，具体涉及一种从航空数字正射影像中提取养殖水塘边界的方法和系统。

背景技术

中国是世界水产养殖第一大国，水产养殖产量占世界一半以上。在中国东部和南部的沿海和内陆省份，分布着大片的人工养殖水塘，用于饲养鱼虾蟹等水产品，是农民重要的经济来源。在海岸及滩涂附近的养殖水塘主要是海水养殖水塘，在内陆的养殖水塘是淡水养殖水塘。这些养殖水塘呈现一定的空间聚集特征，具有明显的空间关联性，在空间分布上呈团聚状态。

在航空测绘中，由于养殖水塘被水覆盖，不能获得池塘对应的地形高程，所以在制作数字高程模型(Digital Elevation Model，DEM)时，必须对养殖水塘计算其二维边界线，绘制断裂线。此外，在大比例尺数字线划图(Digital Line Graphic，DLG)的生产中，也需要对水塘绘制水涯线或岸边线。当使用光学摄影测量相机获得测区的航空数字正射影像(Digital Orthophoto Map，DOM)后，需要从中提取养殖水塘的二维边界线，用于后续DEM和DLG的生产。

从航空DOM中提取养殖水塘边界与传统的从卫星遥感影像提取水体边界存在明显不同。常见的卫星光学多光谱/高光谱遥感影像一般包含可见光波段(红、绿、蓝波段)和近红外波段，可以计算归一化水体指数，筛选归一化水体指数高于一定阈值的像素作为水体像素。此外，星载合成孔径雷达影像也可以用于水体提取，如申请号为202010668117.2的中国专利文献，公开了一种近岸水产养殖塘提取方法，其通过计算双极化水体指数SDWI来提取水体层，然后基于水体层来提取水产养殖塘。但是目前的航空DOM主要由光学相机获得，一般只有红、绿、蓝三个波段，没有近红外波段，且未经过严格的辐射定标和辐射校正，不能计算归一化水体指数。因而，卫星遥感的方法不能直接用于航空DOM的养殖水塘边界提取。

在遥感领域，地物分类和目标检测一般常用基于像素的方法和面向对象的方法。基于像素的方法一般不够稳定，容易在分类中产生椒盐噪点。面向对象的方法已经逐渐成为主流。面向对象的方法将影像分为多个不重叠的对象，认为每个对象内部是均质的，对每个对象进行分类或目标检测。由于在对象水平上可以计算出对象的光谱、纹理和几何特征，因此在分类上可用的特征较多，分类效果更优更稳定。通过设定影像分割和分类规则，可以判断哪些对象为养殖水塘，并给出边界。但是，面向对象的方法涉及到的参数和规则众多，十分复杂。例如，尺度参数、形状因子、紧致度因子、波段光谱差异指数、形状指数等参数含义抽象，其设置均需要反复调整。增加一条规则虽然可能提高用户精度，但可能降低生产者精度。

发明内容

发明目的：本发明提供一种从航空正射影像中提取养殖水塘边界的方法，以提高作业效率和自动化程度。

技术方案：本发明一方面公开了一种从航空正射影像中提取养殖水塘边界的方法，包括：

S1、对包含养殖水塘的航空DOM影像进行栅格化，每个栅格中心点所在的像素点为候选检测点；如果所述候选检测点对应的像素是有效像素，则将该像素加入候选检测点集合S_c-dtc中；

S2、从候选检测点集合S_c-dtc中筛选检测点，构成集合S_dtc，具体包括：

S21、计算S_c-dtc中每个像素点各波段灰度值之和，第k个候选检测点的各波段灰度值之和DN_k为：

其中B为航空DOM影像波段总数，k＝1,2,…,K，K为候选检测点总数；

S22、如果DN_min≤DN_k≤DN_max，则将第k个候选检测点作为检测点；所有检测点构成检测点集合S_dtc；DN_min和DN_max为水塘像素各波段灰度值之和的最小值和最大值；

S3、根据检测点集合S_dtc提取航空DOM中的水塘像素，具体包括S31-S35：

S31、将航空DOM中所有像素点的状态设置为未使用；遍历S_dtc中的每个检测点，对于S_dtc中的第w个检测点p_w，进行如下步骤：

S32、建立种子像素集合S_seed,w和水塘像素集合S_pond,w，并均初始化为空；将p_w加入S_seed,w和S_pond,w，将p_w标记为已使用，建立扩张像素集合S_expand并初始化为空；

S33、计算当前水塘像素集合S_pond,w中像素点各波段灰度值下限dn_b,low,w和上限dn_b,high,w：

如果S_pond,w内像素数目N_pond,w为1，则第b波段灰度值下限dn_b,low,w为影像数据类型可表示的最小值V_min，第b波段灰度值上限dn_b,high,w为影像数据类型可表示的最大值V_max；

如果S_pond,w内像素数目N_pond,w大于1，先对S_pond,w内的像素计算各波段灰度值的平均值

和标准差σ_b,w

如果σ_b,w大于预设的灰度标准差阈值σ_thre，则进入步骤S32处理S_dtc中的下一个检测点；否则，计算第b波段灰度值下限dn_b,low,w：

计算第b波段灰度值上限dn_b,high,w：

其中C_coef为灰度值偏离均值的最大标准差倍数，且C_coef＞0；

如果dn_b,low,w小于V_min，则令dn_b,low,w＝V_min；如果dn_b,high,w大于V_max，则令dn_b,high,w＝V_max；

S34、对S_seed,w内的像素进行八邻域水塘扩张：

遍历S_seed,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列，其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第一扩张条件，将p_m,n加入扩张像素集合S_expand。

所述第一扩张条件为：对任意波段b，b∈{1,2,…,B}，p_m,n在第b波段的值dn_m,n,b均满足dn_b,low,w≤dn_m,n,b≤dn_b,high,w，B为航空DOM影像波段数目；

S35、遍历S_seed,w中的像素点结束后，如果S_expand内的像素点数大于0，则将S_expand内的点全部加入S_pond,w，将S_seed,w更新为S_expand，并清空S_expand，重新执行步骤S33和S34，以更新后的S_seed,w进行下一次水塘扩张；

如果S_expand内的点数为0，则对S_dtc中的第w个检测点p_w的扩张结束，跳转至步骤S32进行下一个检测点的水塘扩张，直到S_dtc中所有检测点均扩张结束；

S4、提取水塘像素的边界，具体包括：

对每个水塘像素集合S_pond,w生成其对应的二值图像I_bin,w，I_bin,w的大小与航空DOM影像相同，如果S_pond,w包含第r行第c列的像素点，则将二值图像I_bin,w第r行第c列的灰度值置为1，否则置为0；对二值图像I_bin,w先进行形态学闭运算，再进行开运算，之后进行边界追踪，将位于边界处的像素作为边界像素集合S_bound,w；

S5、提取水塘阴影，具体包括：

S51、对每个水塘像素集合S_pond,w，计算阴影各波段灰度值上限，第b波段阴影灰度值上限dn_b,shad,w：

C_shad为阴影第b波段灰度值低于无阴影部分均值的最小标准差倍数，为非负数；

S52、对每个水塘像素集合S_pond,w建立阴影种子像素集合S_seed,shad,w并初始化为空，将S_bound,w内的像素全部加入S_seed,shad,w；建立阴影像素集合S_shadow并初始化为空；

S53、对S_seed,shad,w内的像素逐个进行八邻域阴影扩张：

遍历S_seed,shad,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列；其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第二扩张条件，将p_m,n加入阴影像素集合S_shadow；所述第二扩张条件为：p_m,n的所有波段中，灰度值dn_m,n,b满足dn_m,n,b≤dn_b,shad,w的波段数目N_lowband大于等于预设的波段数目阈值N_th；

S54、遍历S_seed,shad,w中的像素点结束后，如果S_shadow内的点数大于0，则将S_shadow内的点全部加入S_pond,w，将S_seed,shad,w更新为S_shadow，并清空S_shadow，重新执行步骤S53，以更新后的S_seed,shad,w进行下一次阴影扩张；

如果S_shadow内的点数为0，则水塘阴影扩张结束；

S55、计算S_pond,w对应的面积A_pond,w：A_pond,w＝N_pond,w|d_xd_y|；其中N_pond,w为S_pond,w中像素点的数量，d_x和d_y分别为每个像素在投影坐标系下横坐标和纵坐标方向代表的长度；

如果A_pond,w满足A_min≤A_pond,w≤A_max，则称S_pond,w对应有效的水塘，将S_pond,w内的全部像素加入S_allpond，其中A_min和A_max分别为水塘面积的下限和上限；

S6、对整幅影像计算水塘边界并简化，具体包括：

S61、对S_allpond内的像素生成对应的二值图像I_allbin，二值图像I_allbin的大小和DOM影像相同，如果S_allpond包含第r行第c列的元素，则将二值图像I_allbin第r行第c列的值置为1，否则置为0；

S62、对二值图像I_allbin进行边界追踪，获得水塘边界；

S63、对水塘边界进行简化，得到提取出的养殖水塘边界。

另一方面，本发明公开了实现上述提取养殖水塘边界方法的系统，包括：

候选检测点集合建立模块1，用于根据步骤S1构建候选检测点集合S_c-dtc；所述步骤S1为：

对包含养殖水塘的航空DOM影像进行栅格化，每个栅格中心点所在的像素点为候选检测点；如果所述候选检测点对应的像素是有效像素，则将该像素加入候选检测点集合S_c-dtc中；

检测点筛选模块2，用于从候选检测点集合S_c-dtc中筛选检测点，建立检测点集合S_dtc；

水塘像素提取模块3，用于根据检测点集合S_dtc提取航空DOM影像中的水塘像素；

水塘像素边界提取模块4，用于提取每个检测点对应水塘像素的边界；

水塘阴影提取模块5，用于根据每个检测点的水塘像素边界提取水塘阴影，并根据水塘像素和阴影像素构建影像中水塘全部像素集合；

影像水塘边界提取模块6，用于对整幅影像计算水塘边界。

有益效果：本发明公开的从航空DOM提取养殖水塘边界的方法中，参数仅包括水塘最小面积A_min、最大面积A_max、波段灰度值之和的最小值DN_min、波段灰度值之和的最大值DN_max、灰度标准差阈值σ_thre、水塘像素灰度值偏离均值的最大标准差倍数C_coef、阴影灰度值低于无阴影部分均值的最小标准差倍数C_shad；参数直观、数量少且容易设置，实现简单，稳定性好，较为稳健。该方法可以有效避免水塘内的增氧机、水草、浅水区、塘梗、塘梗阴影等的影响，能够较为精确地提取出养殖水塘边界。

附图说明

图1为本发明公开的养殖水塘边界提取方法的流程图；

图2为二值图像边界追踪示意图；

图3为典型水塘在某波段灰度值的剖面图；

图4为检测点设置示意图；

图5为养殖水塘边界提取系统的组成示意图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

本发明公开了一种从航空正射影像中提取养殖水塘边界的方法，如图1所示，包括：

S0、估计所述航空DOM影像中养殖水塘的面积

如果

采用人工设置检测点，加入检测点集合S_dtc，进入步骤S3，否则，执行步骤S1；

其中A_DOM为所述航空DOM覆盖的面积；ε为面积比例因子，0<ε<0.5。

本步骤对航空DOM进行粗略估计。对养殖水塘面积占比较小的影像，如果直接采用步骤S1中的栅格化方式获取候选检测点，则候选检测点大部分位于非水塘处，导致提取无效或计算效率偏低。在这种情况下，由人工设置检测点，基本可以确保检测点位于水塘处。本实施例中，令面积比例因子ε的取值为0.2。

S1、对包含养殖水塘的航空DOM进行栅格化，每个栅格中心点所在的像素点为候选检测点；

在投影坐标系下，设航空DOM的横坐标范围为[x_min,x_max]，纵坐标范围为[y_min,y_max]；设每个栅格为边长是d_cell的正方形，则横坐标方向的栅格数N_col为：

纵坐标方向的栅格数N_row为：

为向上取整运算符；第i行第j列的栅格中心点cell_ij(x_dtc,y_dtc)坐标的计算方法是：

其中0≤i<N_row，0≤j<N_col；d_cell可以取单个水塘宽度的最小值，如40m、50m等；

cell_ij(x_dtc,y_dtc)在影像中的像素行列号坐标(n_row,n_col)的计算方法为：

其中，d_x为每个像素在投影坐标系下横坐标方向代表的长度；d_y为每个像素在投影坐标系下纵坐标方向代表的长度，

为向下取整运算符。如果(n_row,n_col)对应的像素不位于影像的无数据(nodata)区域，是有效像素，则将该像素加入候选检测点集合S_c-dtc中。

S2、从候选检测点集合S_c-dtc中筛选检测点，构成集合S_dtc；

本实施例中根据水塘像素各波段灰度值之和来筛选检测点，具体为：

计算S_c-dtc中每个像素点各波段灰度值之和，第k个候选检测点的各波段灰度值之和DN_k为：

其中B为航空DOM影像波段总数，k＝1,2,…,K，K为候选检测点总数；

如果DN_min≤DN_k≤DN_max，则将第k个候选检测点作为检测点；所有检测点构成检测点集合S_dtc；

DN_min和DN_max为水塘像素各波段灰度值之和的最小值和最大值；可以通过在影像中选择多个水塘像素进行统计获得。通过限制各波段灰度值之和的范围来选择水塘检测点，是为了将水塘和塘梗、道路、建筑物等区分开，避免检测点落在塘梗、道路、建筑物等地物上；限制各波段灰度值之和比限制单个波段的灰度值更稳健。通常而言，水塘像素各波段灰度值之和要低于塘梗、道路、建筑物等地物。

S3、根据检测点集合S_dtc提取航空DOM影像中的水塘像素，具体为：

S31、将航空DOM中所有像素点的状态设置为未使用；遍历S_dtc中的每个检测点，对于S_dtc中的第w个检测点p_w，执行如下步骤：

S32、建立种子像素集合S_seed,w和水塘像素集合S_pond,w，并均初始化为空；将p_w加入S_seed,w和S_pond,w，将p_w标记为已使用，建立扩张像素集合S_expand并初始化为空；

S33、计算当前水塘像素集合S_pond,w中像素点各波段灰度值下限dn_b,low,w和上限dn_b,high,w：

如果S_pond,w内像素数目N_pond,w为1，则第b波段灰度值下限dn_b,low,w为影像数据类型可表示的最小值V_min，第b波段灰度值上限dn_b,high,w为影像数据类型可表示的最大值V_max；

例如，对于采用一个字节长度的无符号整型作为存储单像素单波段的灰度值的基本单位的影像，其V_min为0，V_max为255；

如果S_pond,w内像素数目N_pond,w大于1，先对S_pond,w内的像素计算各波段灰度值的平均值和标准差。对第b个波段，其平均值

的计算方法是：

其中dn_j,b为S_pond,w中第j个像素点在第b个波段的值，标准差σ_b,w的计算为：

如果σ_b,w大于预设的灰度标准差阈值σ_thre，则进入步骤S32处理S_dtc中的下一个检测点；否则，计算第b波段灰度值下限dn_b,low,w：

计算第b波段灰度值上限dn_b,high,w：

其中C_coef为灰度值偏离均值的最大标准差倍数，且C_coef＞0；如果dn_b,low,w小于V_min，则设dn_b,low,w＝V_min；如果dn_b,high,w大于V_max，则设dn_b,high,w＝V_max；

灰度标准差阈值σ_thre的取值可以经验性地获取，方法是在DOM上选择多个水塘，对各水塘逐个计算不同波段灰度值的标准差，寻找各波段灰度值标准差的最大值，使σ_thre略大于该最大值，将σ_thre应用到各个波段。对于采用一个字节长度的无符号整型作为存储单像素单波段的灰度值的基本单位的影像，σ_thre取4-6即可。限制σ_thre可以避免检测点位于水塘边缘或水塘外时造成无限制的扩张。

对于C_coef，如果水塘灰度值分布较为均匀，C_coef可以取[2.5,3.5]范围内的值；例如3.0，就足够(在第b波段的灰度值符合正态分布的情况下，99.73％的灰度值分布在3倍标准差的范围内)；如果水塘同时存在浅水区和深水区，灰度值不一致，此时C_coef可以取更大的值，如[6.0,8.0]范围内的值。

S34、对S_seed,w内的像素进行八邻域水塘扩张：

遍历S_seed,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列，其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第一扩张条件，将p_m,n加入扩张像素集合S_expand。

所述第一扩张条件为：对任意波段b，b∈{1,2,…,B}，p_m,n在第b波段的值dn_m,n,b均满足dn_b,low,w≤dn_m,n,b≤dn_b,high,w；

S35、遍历S_seed,w中的像素点结束后，如果S_expand内的像素点数大于0，则将S_expand内的点全部加入S_pond,w，将S_seed,w更新为S_expand，并清空S_expand，重新执行步骤S33和S34，以更新后的S_seed,w进行下一次水塘扩张；

如果S_expand内的点数为0，则对S_dtc中的第w个检测点p_w的扩张结束，跳转至步骤S32进行下一个检测点的水塘扩张，直到S_dtc中所有检测点均扩张结束；

水塘扩张结束后，S_dtc中的每个检测点对应一个水塘像素集合S_pond,w，这些水塘像素基本是影像中水塘无阴影部分的像素。

S4、提取水塘像素的边界，具体包括：

对每个水塘像素集合S_pond,w生成其对应的二值图像I_bin,w，I_bin,w的大小与航空DOM影像相同，如果S_pond,w包含第r行第c列的像素点，则将二值图像I_bin,w第r行第c列的灰度值置为1，否则置为0；

对二值图像I_bin,w先进行形态学闭运算、再进行开运算，之后进行边界追踪，获得位于边界处的像素，构成边界像素集合S_bound,w；如图2所示，图2-(a)中白色部分灰度值为0，即为非水塘像素；斜线填充部分灰度值为1，即为水塘像素；对此二值图像进行边界追踪，获取到边界的像素；图2-(b)中，深灰色部分即为位于边界处的水塘像素。

边界追踪可以采用栅格矢量化的方法，记录构成边界的每个像素；

S5、提取水塘阴影，具体包括：

S51、对每个水塘像素集合S_pond,w，计算阴影各波段灰度值上限，第b波段阴影灰度值上限dn_b,shad,w：

为步骤S33给出的S_pond,w中像素点第b个波段灰度值的平均值，σ_b,w为步骤S33给出的S_pond,w中像素点第b个波段灰度值的标准差；C_shad为阴影灰度值低于

的最小标准差倍数，为非负数；b∈{1,2,…,B}。

考虑到从无阴影部分过渡到阴影部分时，灰度值不是突然剧烈降低，而是在几个像素的宽度内逐步降低，如图3所示，所以C_shad不能取太大的值，以避免后续不能扩张到阴影部分，通常C_shad取大于0，小于1.0的值。

S52、对每个水塘像素集合S_pond,w建立阴影种子像素集合S_seed,shad,w并初始化为空，将S_bound,w内的像素全部加入S_seed,shad,w；建立阴影像素集合S_shadow并初始化为空；

S53、对S_seed,shad,w内的像素逐个进行八邻域阴影扩张：

遍历S_seed,shad,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列；其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第二扩张条件，将p_m,n加入阴影像素集合S_shadow。所述第二扩张条件为：p_m,n的所有波段中，灰度值dn_m,n,b满足dn_m,n,b≤dn_b,shad,w的波段数目N_lowband大于等于预设的波段数目阈值N_th，本实施例中，N_th＝B/2。

在第二扩张条件中不要求p_m,n所有波段的灰度值都满足dn_m,n,b≤dn_b,shad,w，是为了避免一些阴影部分像素在某波段灰度值降低的不够多，导致不能被扩张到的情况。

S54、遍历S_seed,shad,w中的像素点结束后，如果S_shadow内的点数大于0，则将S_shadow内的点全部加入S_pond,w，将S_seed,shad,w更新为S_shadow，并清空S_shadow，重新执行步骤S53，以更新后的S_seed,shad,w进行下一次阴影扩张；

如果S_shadow内的点数为0，则水塘阴影扩张结束。

S55、计算S_pond,w对应的面积A_pond,w：A_pond,w＝N_pond,w|d_xd_y|；其中N_pond,w为S_pond,w中像素点的数量；

如果A_pond,w满足A_min≤A_pond,w≤A_max，则称S_pond,w对应有效的水塘，将S_pond,w内的全部像素加入S_allpond；S_allpond中的像素为航空DOM影像中水塘的全部像素，包括无阴影部分和塘梗及其它地物在水塘造成的阴影部分。

A_min和A_max的取值可以经验性地获取，在DOM上测量多个水塘的面积，寻找最小面积a_min和最大面积a_max。如果一个水塘被阴影分为多个不相连的无阴影部分，则需要将面积最小的无阴影部分的面积作为a_min。如图4所示，某水塘被阴影S分为A1和A2两个无阴影部分，其中A2的面积较小，则取A2的面积作为最小面积a_min。设置A_min＝a_min-α，A_max＝a_max+β，α和β均为正数，且要确保A_min＞0。

S6、对整幅影像计算水塘边界并简化，具体包括：

S61、对S_allpond内的像素生成对应的二值图像I_allbin，二值图像I_allbin的大小和DOM相同，如果S_allpond包含第r行第c列的元素，则将二值图像I_allbin第r行第c列的值置为1，否则置为0；

S62、对二值图像I_allbin进行边界追踪，获得水塘边界；

边界追踪可采用栅格矢量化的方法；

S63、对水塘边界进行简化；

边界简化方法可以采用Douglas-Peucker法、Li-Openshaw算法等。

根据简化后的边界输出水塘边界的矢量文件。

经过上述步骤后，提取出了航空DOM影像中的养殖水塘。但对于某些航空正射影像，可能存在部分水塘边界没有被提取出，或提取不完整的情况。针对这样的情况，可以手工设置检测点，再次采用上述步骤二次提取边界，包括如下步骤：

S7、再次设置检测点，具体包括：

S71、将将水塘边界和DOM叠加显示，对水塘边界进行检查，寻找以下情况的水塘：

(1)未自动提取到边界；

(2)因为水面阴影等原因，水塘边界提取不完整；

S72、在S71中找到的水塘中设置检测点，并记录所设置检测点的横坐标和纵坐标；

如果单个水塘只有一个无阴影部分，只需要在水塘中心选取一个检测点；如果单个水塘被阴影分成多个不相连的无阴影部分，则水塘的每个无阴影部分都需要选取一个检测点，如图4所示，在A1区域设置检测点d1，在A2区域设置检测点d2。检测点要尽量位于无阴影部分的中心。

S8、清空检测点集合S_dtc，将步骤7设置的检测点加入S_dtc，采用步骤S3-S6再次提取水塘边界，具体包括：

S81、清空S_dtc；

S82、将步骤S7生成的检测点文件中的所有检测点加入S_dtc；

S83、按照步骤S3-S6进行自动水塘边界提取；

在步骤S83中，可以对C_coef、C_shad、A_min、A_max等参数进行调整，以获得更好的效果。

S9、将两次边界提取得到的水塘边界进行合并，得到待提取影像中的养殖水塘边界。

之后还可以由人工对提取到的边界进行编辑，包括：

(1)删除非水塘边界；

(2)删除重复的水塘边界；

(3)手工勾绘未自动提取出的水塘边界；

(4)手工调整自动提取的水塘边界；手动调整包括水塘边界点的移动、删除等。

本实施例还公开了实现上述从航空正射影像中提取养殖水塘边界方法的系统，如图5所示，包括：

候选检测点集合建立模块1，用于构建候选检测点集合S_c-dtc，具体步骤为：

S0、估计所述航空DOM中养殖水塘的面积

如果

说明水塘少，由人工设置检测点，否则，执行步骤S1；

S1、对包含养殖水塘的航空DOM影像进行栅格化，每个栅格中心点所在的像素点为候选检测点；如果所述候选检测点对应的像素是有效像素，则将该像素加入候选检测点集合S_c-dtc中；

检测点筛选模块2，用于从候选检测点集合S_c-dtc中筛选检测点，建立检测点集合S_dtc；

水塘像素提取模块3，用于按照步骤3.1-3.5根据检测点集合S_dtc提取航空DOM影像中的水塘像素；

水塘像素边界提取模块4，用于按照步骤4提取每个检测点对应水塘像素的边界；

水塘阴影提取模块5，用于按照步骤5根据每个检测点的水塘像素边界提取水塘阴影，并根据水塘像素和阴影像素构建影像中水塘全部像素集合；

影像水塘边界提取模块6，用于按照步骤6对整幅影像计算水塘边界；

检测点设置模块7，用于按照步骤7获取指定的检测点坐标；

水塘边界二次提取模块8，用于获取根据指定的检测点坐标再次提取到的水塘边界；

水塘边界合并模块9，用于将两次边界提取得到的水塘边界进行合并，得到待提取影像中的养殖水塘边界。

Claims (9)

1.一种从航空正射影像中提取养殖水塘边界的方法，其特征在于，包括：

S1、对包含养殖水塘的航空DOM影像进行栅格化，每个栅格中心点所在的像素点为候选检测点；如果所述候选检测点对应的像素是有效像素，则将该像素加入候选检测点集合S_c-dtc中；

S2、从候选检测点集合S_c-dtc中筛选检测点，具体包括：

S21、计算S_c-dtc中每个像素点各波段灰度值之和，第k个候选检测点的各波段灰度值之和DN_k为：

其中B为航空DOM影像波段总数，k＝1,2,…,K，K为候选检测点总数；

S22、如果DN_min≤DN_k≤DN_max，则将第k个候选检测点作为检测点；所有检测点构成检测点集合S_dtc；DN_min和DN_max为水塘像素各波段灰度值之和的最小值和最大值；

S3、根据检测点集合S_dtc提取航空DOM影像中的水塘像素，具体包括S31-S35：

S31、将航空DOM中所有像素点的状态设置为未使用；遍历S_dtc中的每个检测点，对于S_dtc中的第w个检测点p_w，进行如下步骤：

S32、建立种子像素集合S_seed,w和水塘像素集合S_pond,w，并均初始化为空；将p_w加入S_seed,w和S_pond,w，将p_w标记为已使用，建立扩张像素集合S_expand并初始化为空；

S33、计算当前水塘像素集合S_pond,w中像素点各波段灰度值下限dn_b,low,w和上限dn_b,high,w：

如果S_pond,w内像素数目N_pond,w为1，则第b波段灰度值下限dn_b,low,w为影像数据类型可表示的最小值V_min，第b波段灰度值上限dn_b,high,w为影像数据类型可表示的最大值V_max；

如果S_pond,w内像素数目N_pond,w大于1，先对S_pond,w内的像素计算各波段灰度值的平均值

和标准差σ_b,w；

如果σ_b,w大于预设的灰度标准差阈值σ_thre，则进入步骤S32处理S_dtc中的下一个检测点；否则，计算第b波段灰度值下限dn_b,low,w：

计算第b波段灰度值上限dn_b,high,w：

其中C_coef为灰度值偏离均值的最大标准差倍数，且C_coef>0；

如果dn_b,low,w小于V_min，则令dn_b,low,w＝V_min；如果dn_b,high,w大于V_max，则令dn_b,high,w＝V_max；

S34、对S_seed,w内的像素进行八邻域水塘扩张：

遍历S_seed,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列，其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第一扩张条件，将p_m,n加入扩张像素集合S_expand；

所述第一扩张条件为：对任意波段b，b∈{1,2,…,B}，p_m,n在第b波段的值dn_m,n,b均满足dn_b,low,w≤dn_m,n,b≤dn_b,high,w，B为航空DOM影像波段数目；

S35、遍历S_seed,w中的像素点结束后，如果S_expand内的像素点数大于0，则将S_expand内的点全部加入S_pond,w，将S_seed,w更新为S_expand，并清空S_expand，重新执行步骤S33和S34，以更新后的S_seed,w进行下一次水塘扩张；

如果S_expand内的点数为0，则对S_dtc中的第w个检测点p_w的扩张结束，跳转至步骤S32进行下一个检测点的水塘扩张，直到S_dtc中所有检测点均扩张结束；

S4、提取水塘像素的边界，具体包括：

对每个水塘像素集合S_pond,w生成其对应的二值图像I_bin,w，I_bin,w的大小与航空DOM影像相同，如果S_pond,w包含第r行第c列的像素点，则将二值图像I_bin,w第r行第c列的灰度值置为1，否则置为0；对二值图像I_bin,w先进行形态学闭运算，再进行开运算；之后对二值图像I_bin,w进行边界追踪，获得位于边界处的像素，构成边界像素集合S_bound,w；

S5、提取水塘阴影，具体包括：

S51、对每个水塘像素集合S_pond,w，计算阴影各波段灰度值上限，第b波段阴影灰度值上限dn_b,shad,w：

C_shad为阴影第b波段灰度值低于无阴影部分均值

的最小标准差倍数，为非负数；

S52、对每个水塘像素集合S_pond,w建立阴影种子像素集合S_seed,shad,w并初始化为空，将S_bound,w内的像素全部加入S_seed,shad,w；建立阴影像素集合S_shadow并初始化为空；

S53、对S_seed,shad,w内的像素逐个进行八邻域阴影扩张：

遍历S_seed,shad,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列；其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第二扩张条件，将p_m,n加入阴影像素集合S_shadow；所述第二扩张条件为：p_m,n的所有波段中，灰度值dn_m,n,b满足dn_m,n,b≤dn_b,shad,w的波段数目N_lowband大于等于预设的波段数目阈值N_th；

S54、遍历S_seed,shad,w中的像素点结束后，如果S_shadow内的点数大于0，则将S_shadow内的点全部加入S_pond,w，将S_seed,shad,w更新为S_shadow，并清空S_shadow，重新执行步骤S53，以更新后的S_seed,shad,w进行下一次阴影扩张；

如果S_shadow内的点数为0，则水塘阴影扩张结束；

S55、计算S_pond,w对应的面积A_pond,w：A_pond,w＝N_pond,w|d_xd_y|；其中N_pond,w为S_pond,w中像素点的数量；d_x和d_y分别为每个像素在投影坐标系下横坐标和纵坐标方向的长度；

如果A_pond,w满足A_min≤A_pond,w≤A_max，则称S_pond,w对应有效的水塘，将S_pond,w内的全部像素加入S_allpond，其中A_min和A_max分别为水塘面积的下限和上限；

S6、对整幅影像计算水塘边界并简化，具体包括：

S61、对S_allpond内的像素生成对应的二值图像I_allbin，二值图像I_allbin的大小和DOM影像相同，如果S_allpond包含第r行第c列的元素，则将二值图像I_allbin第r行第c列的值置为1，否则置为0；

S62、对二值图像I_allbin进行边界追踪，获得水塘边界；

S63、对水塘边界进行简化，得到提取出的养殖水塘边界。

2.根据权利要求1所述的提取养殖水塘边界的方法，其特征在于，在步骤S1前包括：

S0、估计所述航空DOM影像中养殖水塘的面积

如果

由人工设置检测点，否则,执行步骤S1；

其中A_DOM为所述航空DOM的面积；ε为面积比例因子，0<ε<0.5。

3.根据权利要求1所述的提取养殖水塘边界的方法，其特征在于，所述灰度标准差阈值σ_thre的取值方法为：

在航空DOM上选择多个水塘，对各水塘逐个计算不同波段灰度值的标准差，寻找各波段灰度值标准差的最大值，令σ_thre大于所述最大值。

4.根据权利要求1所述的提取养殖水塘边界的方法，其特征在于，所述灰度值偏离均值的最大标准差倍数C_coef的取值范围为：

如果水塘灰度值较为均匀，C_coef取值范围为[2.5,3.5]；

如果水塘同时存在浅水区和深水区，C_coef取值范围为[6.0,8.0]。

5.根据权利要求1所述的提取养殖水塘边界的方法，其特征在于，所述步骤S6中采用Douglas-Peucker法或Li-Openshaw算法进行边界简化。

6.根据权利要求1所述的提取养殖水塘边界的方法，其特征在于，步骤S6后还包括：

S7、再次设置检测点，具体包括：

S71、将水塘边界和DOM叠加显示，对水塘边界进行检查，寻找以下情况的水塘：

(1)未自动提取到边界；

(2)水塘边界提取不完整；

S72、在S71中找到的水塘中设置检测点，并记录所设置检测点的横坐标和纵坐标；

S8、清空检测点集合S_dtc，将步骤7设置的检测点加入S_dtc，采用步骤S3-S6再次提取水塘边界；

S9、将两次边界提取得到的水塘边界进行合并，得到待提取影像中的养殖水塘边界。

7.根据权利要求6所述的提取养殖水塘边界的方法，其特征在于，还包括：由人工对提取到的边界进行编辑，包括：

(1)删除非水塘边界；

(2)删除重复的水塘边界；

(3)手工勾绘未自动提取出的水塘边界；

(4)手工调整自动提取的水塘边界。

8.一种从航空正射影像中提取养殖水塘边界的系统，其特征在于，包括：

候选检测点集合建立模块(1)，用于根据步骤S1构建候选检测点集合S_c-dtc；所述步骤S1为：

对包含养殖水塘的航空DOM影像进行栅格化，每个栅格中心点所在的像素点为候选检测点；如果所述候选检测点对应的像素是有效像素，则将该像素加入候选检测点集合S_c-dtc中；

检测点筛选模块(2)，用于从候选检测点集合S_c-dtc中筛选检测点，建立检测点集合S_dtc；

水塘像素提取模块(3)，用于根据检测点集合S_dtc提取航空DOM影像中的水塘像素；

水塘像素边界提取模块(4)，用于提取每个检测点对应水塘像素的边界；

水塘阴影提取模块(5)，用于根据每个检测点的水塘像素边界提取水塘阴影，并根据水塘像素和阴影像素构建影像中水塘全部像素集合；

影像水塘边界提取模块(6)，用于对整幅影像计算水塘边界；

所述水塘像素提取模块(3)根据检测点集合S_dtc提取航空DOM中的水塘像素的步骤为：

S31、将航空DOM中所有像素点的状态设置为未使用；遍历S_dtc中的每个检测点，对于S_dtc中的第w个检测点p_w，进行如下步骤：

S32、建立种子像素集合建立种子像素集合S_seed,w和水塘像素集合S_pond,w，并均初始化为空；将p_w加入S_seed,w和S_pond,w，将p_w标记为已使用，建立扩张像素集合S_expand并初始化为空；

S33、计算当前水塘像素集合S_pond,w中像素点各波段灰度值下限dn_b,low,w和上限dn_b,high,w：

如果S_pond,w内像素数目N_pond,w为1，则第b波段灰度值下限dn_b,low,w为影像数据类型可表示的最小值V_min，第b波段灰度值上限dn_b,high,w为影像数据类型可表示的最大值V_max；

如果S_pond,w内像素数目N_pond,w大于1，先对S_pond,w内的像素计算各波段灰度值的平均值

和标准差σ_b,w；

如果σ_b,w大于预设的灰度标准差阈值σ_thre，则进入步骤S32处理S_dtc中的下一个检测点；否则，计算第b波段灰度值下限dn_b,low,w：

计算第b波段灰度值上限dn_b,high,w：

其中C_coef为灰度值偏离均值的最大标准差倍数，且C_coef>0；

如果dn_b,low,w小于V_min，则令dn_b,low,w＝V_min；如果dn_b,high,w大于V_max，则令dn_b,high,w＝V_max；

S34、对S_seed,w内的像素进行八邻域水塘扩张：

遍历S_seed,w中的像素点，设当前待扩张像素点为p_r,c，位于航空DOM影像的第r行第c列，其八邻域像素为p_m,n，|m-r|≤1，|n-c|≤1；如果p_m,n状态为已使用，则跳过p_m,n，处理p_r,c的下一个八邻域像素；否则，将p_m,n状态设置为已使用，如果p_m,n满足第一扩张条件，将p_m,n加入扩张像素集合S_expand；

所述第一扩张条件为：对任意波段b，b∈{1,2,…,B}，p_m,n在第b波段的值dn_m,n,b均满足dn_b,low,w≤dn_m,n,b≤dn_b,high,w，B为航空DOM影像波段总数；

S35、遍历S_seed,w中的像素点结束后，如果S_expand内的像素点数大于0，则将S_expand内的点全部加入S_pond,w，将S_seed,w更新为S_expand，并清空S_expand，重新执行步骤S33和S34，以更新后的S_seed,w进行下一次水塘扩张；

如果S_expand内的点数为0，则对S_dtc中的第w个检测点p_w的扩张结束，跳转至步骤S32进行下一个检测点的水塘扩张，直到S_dtc中所有检测点均扩张结束。

9.根据权利要求8所述的从航空正射影像中提取养殖水塘边界的系统，其特征在于，还包括：

检测点设置模块(7)，用于获取指定的检测点坐标；

水塘边界二次提取模块(8)，用于获取根据指定的检测点坐标再次提取到的水塘边界；

水塘边界合并模块(9)，用于将两次边界提取得到的水塘边界进行合并，得到待提取影像中的养殖水塘边界。

Images