一种基于遥感图像的浅海沙波信息提取方法与流程

本发明属于遥感技术应用和浅海地形监测领域，涉及一种浅海沙波信息的提取方法，具体是一种能够利用遥感图像中所含的灰度信息提取浅海水下沙波的位置、长度和波向等信息的方法。

背景技术：
海底地形测绘是探索和研究海洋的第一步，是开发和保护海洋的首要条件，也是现今海洋学研究的迫切需求。在回声测深仪发明之前，主要靠测深杆和测深锤来测量水深，测量精度较差。20世纪20年代回声测深仪出现之后，现代意义上海图的绘制才得以实现。但早期的测深仪为单波束发射，一次发射只能得到测量船正下方的水深，因而只能实现点、线测量，无法反映测线之间的地形地貌。1970年代中期出现的多波束测深技术实现了带状测量（最大宽度可达水深的7倍），显著提高了测深效率。尽管如此，限于测量周期长、人力消耗大和资金需求高等方面的劣势，在地形高变化的浅海区域进行大范围的动态变化监测任务中，仅依赖多波束测深方法仍然很难满足需求。对于浅海区域广泛分布的浅海沙波地形动态变化监测，需要掌握沙波的位置、长度和波向等基础信息。如何充分发挥遥感技术的大范围高动态观测能力，对浅海沙波区的沙波信息进行自动提取，实现动态监测，对于掌握地形演化规律和特征都有十分重要的意义。本发明针对广泛分布的浅海沙波分布区进行地形变化动态监测的需求，利用遥感技术方法的大范围覆盖和高频率重访的能力，基于遥感图像中所包含的因沙波地形对水体流场和海面粗超度的调制而在呈现的亮暗条纹特征，引入图像二维像元倾斜方向，可高效准确地提取浅海沙波的位置、长度和波向信息，可用于浅海沙波地形动态变化监测。从检索到的公开资料看，尚未有基于遥感图像，通过图像二维像元倾斜方向进行浅海沙波信息提取的相关报道。

技术实现要素：
本发明的目的是提供一种高效准确、便捷低耗的基于遥感图像的浅海沙波信息提取方法。本发明通过以下技术方案实现：一种基于遥感图像的浅海沙波信息提取方法，包括以下步骤：（1）计算含有沙波信息的遥感图像中每一像元的二维倾斜方向，获得遥感图像对应的像元倾斜方向分布图；（2）标记遥感图像中最长的亮条纹为第一条纹，最长的暗条纹为第二条纹，从第一条纹向第二条纹做第一辅助线，使其尽可能地垂直于第一条纹及第二条纹，该辅助线与正北方向的夹角记为A0，分割阈值A1=A0-T，A2=A0+T，参数T为30-45°，根据A1、A2的大小确定取值范围，遥感图像中二维倾斜方向位于取值范围内的像元的集合即为初步沙波区；（3）在遥感图像上采用矢量化方法获得与初步沙波区相应的初步沙波面状图斑；（4）对获得的初步沙波面状图斑进行面积和形状过滤，得到最终沙波面状图斑；（5）提取最终沙波面状图斑的中心线作为沙波曲线，获得沙波位置信息；（6）计算沙波曲线的长度获得沙波的长度信息；（7）用直线连接沙波曲线的两个端点作为第二辅助线，计算该辅助线与正北方向的夹角α，进而获得沙波的波向β=α+90°。上述技术方案中，步骤（1）中所述的计算遥感图像上的每一个像元的二维倾斜方向，具体方法如下：对于遥感图像最边缘的像元，不计算其二维倾斜方向，直接记为空值；对于遥感图像上的其余像元，将当前像元及其周围相邻的8个像元共计9个像元的灰度值从上到下、从左到右依次记为a、b、c、d、e、f、g、h、i，则当前像元的二维倾斜方向A的计算公式为：其中，p为临时参数，A的值域范围为0~360，单位为度。步骤（2）中对步骤（1）所计算获得的像元方向分布图进行阈值分割。考察遥感图像上沙波纹理分布特征，沙波通常呈现亮暗条纹间隔分布，根据沙波地形遥感成像机理，亮暗条纹相接处为沙波的位置。从垂直于沙波走向的剖面来看，遥感图像上呈现从背景灰度逐渐变亮（变暗），到达最亮（最暗）处，突然变成最暗（最亮），然后又逐渐变亮（变暗）到背景灰度，其中从最亮（最暗）处突变为最暗（最亮）的地方就是沙波所在位置。可见，沙波所在位置和背景区域的像元灰度变化趋势是相反的。比如，沙波位置是最亮到最暗变化时，则背景区域为从暗到亮变化；相反，沙波位置是最暗到最亮变化时，则背景区域是亮到暗变化。因此，沙波区域像元倾斜方向A的值域范围不同于背景区域。根据遥感图像上最大沙波（即目视判定最长的亮条纹和暗条纹）的位置，尽可能垂直于该最长的亮暗条纹从该亮条纹到该暗条纹绘制第一辅助线，确定该辅助线与正北方向的夹角为初始阈值A0，确定分割阈值A1和A2，A1=A0-T，A2=A0+T，其中T通常为30~45度之间的一个值。根据A1和A2的值确定取值范围，具体确定方法如下：①A1≥0°，且A2≤360°时，取值范围为[A1,A2]；②A1≥0°，且A2＞360°时，取值范围为[A1,360°]∪[0，A2-360°]；③A1＜0°，且A2≤360°时，取值范围为[360°+A1,360°]∪[0，A2]。遥感图像中二维倾斜方向位于取值范围内的像元的集合即为初步沙波区；步骤（3）中所述的矢量化方法为边缘跟踪矢量化方法，采用地理信息系统软件（ArcGIS）实现。步骤（4）中所述的对初步沙波面状图斑进行面积和形状过滤，具体方法如下：保留初步沙波面状图斑中长轴为10-20倍像元尺寸且面积为20-50个像元面积的图斑，获得最终沙波面状图斑。步骤（5）中所述的提取最终沙波面状图斑的中心线，采用欧几里德距离均分法，其具体操作步骤如下：1）在基本垂直于沙波面状图斑长轴的方向，均匀设置若干条直线辅助线，直线辅助线间的间距为遥感图像像元分辨率r的3-5倍；2）取每条直线辅助线与沙波面状图斑重叠部分的中点，记为该直线辅助线的中心点，将所有直线辅助线的中心点连接，获得沙波曲线，即为沙波位置信息。具体来说，广泛分布的浅海沙波地形能在遥感图像上成像，呈现为清晰的亮暗条幅间或分布的沙波纹理信息。其原理是：浅海沙波区起伏的地形通过水动力作用调制了海表流场，被调制后的海表流场通过波-流相互作用改变了表面波的密度谱，表面波密度谱的改变导致海面粗糙度的改变，海面粗糙度的改变引起海-气界面太阳耀光散射强度的改变，从而在遥感图像上呈现亮暗条纹间或分布的沙波波纹信息。遥感图像的亮暗条纹相间处即指示了沙波的位置所在。而沙波地形作为浅海区域经过长期水动力过程作用形成的自然地貌形态，呈现波浪形连续分布，对于特定区域范围内的沙波通常呈排列分布，沙波走向基本接近，因此遥感图像上亮暗条纹的走向也相近。上述遥感图像上的沙波波纹特征信息和沙波地形分布规律，可用于浅海沙波信息提取。本发明利用浅海沙波地形在遥感图像中呈现的特征信息（亮暗条纹间或分布）和沙波地形分布规律，引入图像二维像元倾斜方向，通过遥感图像灰度空间分析计算，获得沙波的位置、长度和波向等信息。本发明的有益效果是：浅海沙波分布区地形通常处于动态变化中，沙波位置和形态变化情况是沙波地形监测的重要内容。多波束测深虽已成为当前浅海水深测量的主要手段，但受测量周期、人力消耗和资金的限制，通常无法经常进行动态测量与监测。遥感监测是一种非常有效的沙波地形变化动态监测手段。本发明针对广泛分布的浅海沙波区的地形变化监测需求，利用遥感技术方法的大范围覆盖和高频率重访的能力，基于遥感图像中所包含的因沙波地形对水体流场和海面粗超度的调制而在呈现的亮暗条纹特征，引入图像二维像元倾斜方向，可高效准确地提取浅海沙波的位置、长度和波向信息，可用于浅海沙波管理与地形动态变化监测，是遥感信息技术用于沙波地形变化监测的一项创新，具有极大的实用价值。附图说明图1是基于遥感图像的浅海沙波信息提取方法的技术路线图；图2是计算像元的二维倾斜方向时像元分布及对应灰度值示意图；图3是遥感图像中垂直于沙波走向的典型灰度变化剖面线；图4是基于像元倾斜方向分割阈值获得的初步沙波区域；图5是经过面积和形状筛选获得的最终沙波面状图斑；图6是基于欧几里德距离均分法提取中心线示意图；图7是提取的沙波曲线信息；图8是沙波波向计算示意图。具体实施方式下面结合附图对本发明做进一步阐述。如图1所示，本发明的基于遥感图像的浅海沙波信息提取方法具体包括如下步骤：（1）计算含有沙波信息的遥感图像中每一像元的二维倾斜方向，获得遥感图像对应的像元倾斜方向分布图：收集含有沙波信息的遥感图像，针对遥感图像上的每一个像元，根据其相邻像元的灰度值，参考地形坡向，引入二维像元倾斜方向，计算获得像元倾斜方向分布图。每一个像元的二维倾斜方向，具体方法如下：对于遥感图像最边缘的像元，不计算其二维倾斜方向，直接记为空值；对于遥感图像上的其余像元，将当前像元及其周围相邻的8个像元共计9个像元的灰度值从上到下、从左到右依次记为a、b、c、d、e、f、g、h、i，如图2所示，则当前像元的二维倾斜方向A的计算公式为：其中，p为临时参数，A的值域范围为0~360，单位为度。（2）对步骤（1）获得的像元倾斜方向分布图进行阈值分割，考察遥感图像上沙波纹理分布特征，沙波通常呈现亮暗条纹间隔分布，根据沙波地形遥感成像机理，亮暗条纹相接处为沙波的位置。从垂直于沙波走向的剖面来看，如图3所示，遥感图像上呈现从背景灰度逐渐变暗（变亮），到达最暗（最亮）处，突然变成最亮（最暗），然后又逐渐变暗（变亮）到背景灰度，其中从最暗（最亮）处突变为最亮（最暗）的地方就是沙波所在位置。可见，沙波所在位置和背景区域的像元灰度变化趋势是相反的。比如，沙波位置是最亮到最暗变化时，则背景区域为从暗到亮变化；相反，沙波位置是最暗到最亮变化时，则背景区域是亮到暗变化。因此，沙波区域像元倾斜方向A的值域范围不同于背景区域。根据遥感图像上目视判定最长的亮条纹和暗条纹（即最大沙波）的位置，将该亮条纹作为第一条纹，该暗条纹作为第二条纹，尽可能垂直于第一、二条纹，从第一条纹到第二条纹绘制第一辅助线，确定该辅助线与正北方向的夹角为初始阈值A0，确定分割与阈值A1和A2，使得A1=A0-T，A2=A0+T，其中T为一个阈值范围，通常为30~45度之间的一个值。根据A1和A2的值确定取值范围，确定方法如下：①A1≥0°，且A2≤360°时，取值范围为[A1,A2]；②A1≥0°，且A2＞360°时，取值范围为[A1,360°]∪[0，A2-360°]；③A1＜0°，且A2≤360°时，取值范围为[360°+A1,360°]∪[0，A2]。遥感图像中二维倾斜方向位于取值范围内的像元的集合即为初步沙波区，如图4所示。（3）在遥感图像上采用矢量化方法获得与初步沙波区相应的初步沙波面状图斑；所述的矢量化方法可以为边缘跟踪矢量化方法，采用地理信息系统软件（ArcGIS）实现。（4）对获得的初步沙波面状图斑进行面积和形状过滤，具体方法如下：保留初步沙波面状图斑中长轴为10-20倍像元尺寸且面积为20-50个像元面积的图斑，获得最终沙波面状图斑，如图5所示。（5）采用欧几里德距离均分法提取最终沙波面状图斑的中心线作为沙波曲线，获得沙波位置信息；如图6所示，其具体操作步骤如下：1）在基本垂直于沙波面状图斑长轴的方向，均匀设置若干条直线辅助线，直线辅助线间的间距为遥感图像像元分辨率r的3-5倍；2）取每条直线辅助线与沙波面状图斑重叠部分的中点，记为该直线辅助线的中心点，将所有直线辅助线的中心点连接，获得沙波曲线，如图7所示，即获得沙波位置信息。（6）计算沙波曲线的长度获得沙波的长度信息；（7）如图8所示，用直线连接沙波曲线的两个端点作为第二辅助线，计算第二辅助线与正北方向的夹角α，进而计算沙波的波向为β=α+90°。完成沙波信息提取。以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均包含在本发明的保护范围之内。