一种融合海浪特征和光谱特征的浅海水深遥感提取方法与流程

本发明涉及一种海洋遥感技术应用和海底地形测绘领域，具体是一种能够利用高分辨率多光谱遥感图像中的海浪特征信息和光谱特征信息，如所希望那样提取浅海水深信息的方法。

背景技术：

海底地形测绘是探索和研究海洋的第一步，是开发和保护海洋的首要条件，也是现今海洋学研究的迫切需求。在回声测深仪发明之前，主要靠测深杆和测深锤来测量水深，测量精度较差。20世纪20年代回声测深仪出现之后，现代意义上海图的绘制才得以实现。但早期的测深仪为单波束发射，一次发射只能得到测量船正下方的水深，因而只能实现点、线测量，无法反映测线之间的地形地貌。1970年代中期出现的多波束测深技术实现了带状测量(最大宽度可达水深的7倍)，显著提高了测深效率。尽管如此，在大范围浅海水深测量任务中，仍存在测量周期长、人力消耗大和资金需求高等劣势。

随着遥感技术的出现和发展，根据光学遥感图像包含的海洋水体散射和海底反射特征信息与水深的关系，研究提出了一系列光学遥感测深方法，包括基于单波段或多波段的统计模型和半经验模型等方法。但因为该方法基于水体的光谱信息进行水深探测，因此，其适用性受水体浑浊度、海面波动和大气吸收的影响与限制，仅适用于清澈平静的海域条件，且通常需要一定数量的实测地形数据进行控制与校正。此外，当光学遥感器、太阳和观测海面呈现特定的空间几何关系时，太阳通过海面的镜面直接进入遥感器，而在遥感图像上呈现太阳耀光，此时的太阳耀光信息会影响准确提取水体光谱信息。然而，包含太阳耀光信息的光学遥感图像，能记录海面的起伏和粗糙度特征。当遥感图像的空间分辨率足够高时，海浪在光学遥感图像中呈现亮暗相间的条纹。根据近岸浅水波传输理论，海浪从水深向浅海传播时，受水深的影响，海浪波长变短，且海浪的波长、频率与水深存在定量关系。基于该原理，目前已有采用连续拍摄(间隔10秒)的多景光学遥感图像提取海浪传播速度和频率，进而进行水深测量的方法，但连续多景遥感图像的获取成本很高。

本发明针对近岸浅海地形快速、高效、经济的测绘需求，利用单景高分辨率光学遥感图像，通过太阳耀光分离，获得含有海面海浪信息的太阳耀光图像和不含耀光信息的多光谱图像。然后基于太阳耀光图像测量海浪波长及其变化信息，基于浅水波理论提取浅海水深信息。以此为地形控制数据，建立多光谱图像的水深反演模型，提取整景图像浅海水深，从而实现无现场实测数据支持的浅海水深反演。

技术实现要素：

本发明的目的是提供了一种新的融合海浪特征和光谱特征的浅海水深遥感提取方法。

本发明通过以下技术方案实现：

一种融合海浪特征和光谱特征的浅海水深遥感提取方法，其特征在于：包括以下步骤：

(1)获得包含海面海浪信息的高分辨率多光谱遥感图像，采用几何校正方法进行遥感图像定位配准，使其具有较准确的地理空间参考；

(2)分解太阳耀光信息和水体辐射信息；

(3)选定远岸深水区海浪波长测量参考线和子图；

(4)计算参考线上每一个子图的海浪波长；

(5)确定远岸深水区海浪波长；

(6)选定一定数量的测量点；

(7)计算每一个测量点的海浪波长；

(8)计算每一个测量点的水深；

(9)提取测量点的多光谱水体辐射值；

(10)建立多光谱水深反演模型；

(11)反演整景遥感图像的浅海水深信息。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(2)中进行太阳耀光信息和水体辐射信息的分解，获得太阳耀光信息图像和多光谱信息图像；具体分解方法如下：

1)选取遥感图像中浅海区域；

2)提取多光谱图像中的一个可见光波段的辐射强度l(visn)(n＝1、2、3，为可见光波段序号)和近红外波段的辐射强度l(nir)，根据l(visn)和l(nir)的对应空间位置获得二维散点图；

3)统计散点图的线性趋势线并计算倾斜角度

4)计算获得太阳耀光信息图像ls(visn)′和水体辐射信息图像lw(visn)′，分解公式分别为：

5)重复步骤2)到步骤4)，分解逐个可见光波段的太阳耀光信息和水体辐射信息，获得全部可见光波段水体辐射信息图像。其中，步骤(3)中以步骤(2)分解的太阳耀光信息图像上离海岸最远的位置为起点，沿尽可能垂直于海岸的参考线，以50-100个像素为步长(r)依次截取系列子图，直到海岸，子图的大小为n×n像素，其中n为例如64、128、256、512等2的指数。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(4)中对步骤(3)确定的每一个子图进行计算，获得其所代表的海浪波长；海浪波长计算采用空间相关法，具体计算步骤如下：

1)以子图[记为g(x，y)，大小为n像素]的中心坐标(i，j)为原点，分别向x、y方向移动距离r和c，r和c的取值为(-n/2，n/2)，获得匹配图像g′(x，y)，计算子图与匹配图像间的相关系数；其中，n是表示子图的大小，即像素个数；

两影像间的相关系数计算公式为：

其中，

2)完成全部匹配图像，即r和c的取值覆盖全部的(-n/2，n/2)，获得相关系数图，图像大小为n×n像素；

3)采用阈值分割法进行相关系数图二值化处理，确定相关系数大于等于阈值p的为1，小于阈值p的为0；p的大小根据实际相关系数值分布确定，具体方法为：先对子图相关系数图的相关系数从高到低进行统计像素数量，当累计像素数量为子图像素数量50％时所对应的相关系数作为分割阈值，取值为0.2-0.4之间；其中，相关系数最高为1，相关系数最低为0；

4)采用形态学算子进行骨架提取，得到显著代表海浪波峰的骨架线。形态a被b腐蚀的含义为：

其中，s为约束参数，a为待处理斑块的形态(如阈值分割后的一个海浪波峰所在的斑块)，b表示大小为b调节腐蚀形态，选取边长为b个像素的正方形斑块；用b调节腐蚀a的过程是先对b做关于原点的映射，再将其映像平移x后与a^c的交集至少含s个像素点，随后取结果的补集即为腐蚀结果；即用含约束参数s的b调节腐蚀a得到的集合，是b的反射与a^c至少相交s个像素点时b的原点位置集合的补集；

5)根据海浪波峰骨架线方向作垂直参考线，依次计算获得参考线与相邻两条骨架线交叉的距离，作为海浪波长的测量值，记为波长序列：v1＝{v1，v2，…，vn}；为保证数据质量，分别向左右两个方向平移参考线，进行多次采样，得到波长的序列集合v＝{v1，v2，...vn}，计算总体平均值，作为该子图的海浪波长l。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(5)中，对步骤(4)计算的每个子图的海浪波长序列，由第一个子图开始，以8-10个子图为一组进行分组(如分为t组)，计算每一组子图的海浪波长的算术平均作为该组的平均波长然后依此计算两组相邻的平均波长差的绝对值：

当显著小于时(即是的1.5倍到3倍时)，就到了深水区的边缘，即将第1到n组的子图的海浪波长的算术平均作为远岸深水区的海浪波长ld，单位：m；

如不满足此条件，则无法确定深水区，该遥感图像无法用于浅海水深反演。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(7)中根据步骤(6)确定的每一个测量点为中心点，选取大小为n×n像素的子图，其中n为2的指数，然后采用步骤(4)的方法计算获得每个测量点的海浪波长l；单位：m。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(8)中根据步骤(7)中计算获得的测量点的海浪波长l和步骤(5)确定的深水区海浪波长ld，计算测量点的水深，计算方法如下：

其中，h为测量点水深值(单位：m)，l为测量点的海浪波长，ld为步骤(5)确定深水区海浪波长。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(9)中根据步骤(6)确定的测量点的空间位置，从由步骤(2)分解处理获得的多光谱水体辐射信息图像lw(visn)′中，提取蓝光波段和绿光波段的水体辐射强度值，单位为w*m^-2*nm^-1*sr^-1，分别记为l(λ1)和l(λ2)；其中蓝光波段的中心波长为480nm，绿光波段的中心波长为545nm。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法的步骤(10)中将通过步骤(8)计算的测量点水深值和步骤(9)提取的对应位置的蓝光波段和绿光波段水体辐射值，进行对数比值模型统计分析，统计模型表达式为：

其中：k为一个修正系数，以保证kl(λ1)和kl(λ2)的值大于1，通常取1000；m0和m1是模型待定系数；h为步骤(8)计算的测量点水深；l(λ1)和l(λ2)分别为蓝光波段和绿光波段的水体辐射强度；通过对测量点的测量水深和多波段水体辐射强度值的线性拟合，确定m0和m1，从而建立该图像的水深反演模型。

作为优选，上述一种融合海浪特征和光谱特征的浅海水深遥感提取方法，其特征在于步骤(11)中将通过步骤(9)获得的水深反演模型，用于步骤(2)分解获得的整景多光谱水体辐射信息图像lw(visn)′的水深计算；依据逐像元的蓝光波段和绿光波段辐射强度值，计算获得整景遥感图像的浅海水深信息。

有益效果：多波束和单波束测深已成为当前浅海水深测量的主要手段，但受测量周期、人力消耗和资金的限制，尚有很多区域无法进行全覆盖测量，更新测量更是不足，另外对于近岸5m以浅的极浅区域，更是难以有效测量。本发明针对浅海区域的地形测绘需求，根据浅水波理论和水体多光谱辐射传输理论，利用高分辨率多光谱遥感图像上的海浪特征信息和光谱特征信息，建立了基于单景高分辨率光学遥感图像的浅海水深提取方法，可实现无现场实测水深数据支持下，仅由一景高分辨率光学遥感图像就可以提取整景遥感图像的浅海区域的水深信息，从而降低了浅海水深测量成本和对外部参数的依赖，是遥感信息技术应用方面的一项创新，是对近岸浅海水深测量的一个有益补充，具有极大的实用价值。

附图说明

图1是融合海浪特征和光谱特征的浅海水深遥感提取方法技术路线示意图；

图2是太阳耀光信息和水体辐射信息分解示意图；

图3是工作区、测量点和远岸深水区海浪波长测量位置示意图；

图4是子图空间相关性分析示意图；

图5是子图空间相关性二值化示意图；

图6是子图海浪波峰骨架线提取结果示意图；

图7是工作区水深提取结果示意图；

图2中，1—以近红外波段和可见光波段为坐标的辐射强度散点；2—全部散点的线性拟合直线；3—与2平行且过(0,0)点的参考线；4—线性拟合直线2的倾斜角。

图3中，1—陆地区域；2—工作区边界；3—工作区测量点；4—远岸深水区观测参考线；5—远岸深水区观测子图起始点；6—系列子图。

图4中，1—子图范围；2—子图的中心点(i,j)；3—匹配图像范围；4—匹配图像中心点(i+r,j+c)。

具体实施方式

下面结合附图对本发明的实施作具体说明：

实施例1

根据本发明的融合海浪特征和光谱特征的浅海水深遥感提取方法进行实验，如图1所示，具体包括如下步骤：

(1)获得包含海面海浪信息的高分辨率光学遥感图像，采用几何校正方法进行遥感图像定位配准，使其具有较准确的地理空间参考：

选取工作区的高分辨率光学遥感图像1景，检查遥感图像上海浪纹理信息的清晰度和完整性。然后进行遥感图像几何校正，通常采用地面控制点同名点方式，即在遥感图像上的陆地区域选择一定数量的相同地物点，建立控制方程。地面控制点法回避传感器成像时的实际几何状态，直接对遥感图像进行几何校正，该方法的校正精度依赖于地面控制点精度。典型的地面控制点校正法是多项式校正，将遥感图像的总体变形看作是平移、缩放、旋转、偏扭以及更高次的基本变形综合作用结果，校正前后图像相应点之间坐标关系可以用一个适当的多项式来表达。遥感图像几何校正完成后，可选取适当数量的精确测量点进行评估，确保遥感图像具有较准确的地理空间参考。

(2)分解太阳耀光信息和水体辐射信息；

分解太阳耀光信息和水体辐射信息的具体分解方法如下：1)选取遥感图像中浅海区域；2)提取多光谱图像中的一个可见光波段的辐射强度l(visn)(n＝1、2、3，为可见光波段序号)和近红外波段的辐射强度l(nir)，根据l(visn)和l(nir)的对应空间位置获得二维散点图(图2)；3)统计散点图的线性趋势线(图2中的2)并计算倾斜角度(图2中的4)；4)计算获得太阳耀光信息图像ls(visn)′和水体辐射信息图像lw(visn)′，分解公式分别为：

5)重复步骤2)到步骤4)，分解逐个可见光波段的太阳耀光信息和水体辐射信息，获得全部可见光波段的水体辐射信息图像。

(3)选定远岸深水区海浪波长测量参考线和子图；

以太阳耀光信息图像上离海岸最远的位置为起点(图3中的5)，沿尽可能垂直于海岸的参考线(图3中的4)，以一定像素(50-100个像素)为步长(r)依次截取系列子图(图3中的6)，直到海岸，子图的大小为n×n像素，其中n为例如64、128、256、512等2的指数。

(4)计算参考线上每一个子图的海浪波长；

对上一步确定的每一个子图进行空间相关分析法计算，获得其所代表的海浪波长。具体计算步骤如下：

1)以子图(图4中1)的中心坐标(i，j)为原点(图4中2)，分别向x、y方向移动r和c，r和c的取值为(-n/2，n/2)，获得匹配图像g′(x，y)(图4中3)，计算子图与匹配图像间的相关系数。两影像间的相关系数计算公式为；

其中，

2)完成全部匹配图像，即r和c的取值覆盖全部的(-n/2，n/2)，获得相关系数图，图像大小为n×n像素。

3)采用阈值分割法进行相关系数图二值化处理(图5)，确定相关系数大于等于阈值p的为1，小于阈值p的为0。p的大小根据实际相关系数值分布确定，一般为0.2-0.4之间。

4)采用形态学算子进行骨架提取，得到显著代表海浪波峰的骨架线(图6)。形态a被b腐蚀的含义为：

其中，s为约束参数，b调节腐蚀a的过程是先对b做关于原点的映射，再将其映像平移x后与a^c的交集至少含s个像素点，随后取结果的补集即为腐蚀结果。即用含约束参数s的b调节腐蚀a得到的集合，是b的反射与a^c至少相交s个像素点时b的原点位置集合的补集。

5)根据海浪波峰骨架线方向作垂直参考线，依次计算获得参考线与相邻两条骨架线交叉的距离，作为海浪波长的测量值，记为波长序列：v1＝{v1，v2，...，vn}。为保证数据质量，分别向左右两个方向平移参考线，进行多次采样，得到波长的序列集合v＝{v1，v2，...vn}，计算总体平均值，作为该子图的海浪波长l。

(5)确定远岸深水区海浪波长；

由第一个子图开始，以8-10个子图为一组进行分组(如分为t组)，计算每一组子图的海浪波长的算术平均作为该组的平均波长然后依此计算两组相邻的平均波长差的绝对值：

当显著小于时(即是的1.5倍到3倍时)，就到了深水区的边缘，即将第1到n组的子图的海浪波长的算术平均作为远岸深水区的海浪波长ld(单位：m)。如不满足此条件，则无法确定深水区，该遥感图像无法用于浅海水深提取。

(6)选定一定数量的测量点；

在太阳耀光信息图像ls(visn)′上，根据海面海浪信息分布特征，在海浪信息清晰的区域选择测量点(图3中3)，测量点的位置尽可能均匀地分布于工作区域内，且测量点的数量为20-50个。

(7)计算每一个测量点的海浪波长；

以选定的测量点为中心，选取大小为n×n像素的子图，其中n为例如64、128、256、512等2的指数。然后采用步骤(4)的海浪波长测量方法获得每个测量点的海浪波长l(单位：m)。

(8)计算每一个测量点的水深；

根据测量点的海浪波长l和深水区海浪波长ld，计算测量点的水深，计算方法如下：

其中，h为测量点水深值(单位：m)，l为测量点的海浪波长，ld为步骤(5)确定深水区海浪波长。

(9)提取测量点的多光谱水体辐射值；

根据测量点的空间位置，从多光谱水体辐射信息图像lw(visn)′中，提取蓝光波段(中心波长为480nm)和绿光波段(中心波长为545nm)的水体辐射强度值(单位为w*m^-2*nm^-1*sr^-1)，分别记为l(λ1)和l(λ2)

(10)建立多光谱水深反演模型；

将测量点水深值和对应位置的蓝光波段和绿光波段水体辐射值，进行对数比值模型统计分析，统计模型表达式为：

其中：k为一个修正系数，以保证kl(λ1)和kl(λ2)的值大于1，通常取1000；m0和m1是模型待定系数；h为测量点水深；l(λ1)和l(λ2)分别为蓝光波段和绿光波段的水体辐射强度。通过对测量点的测量水深和多波段水体辐射强度值的线性拟合，确定m0和m1，从而建立该图像的水深反演模型。

(11)反演整景遥感图像的浅海水深信息。

将获得的水深反演模型，用于分解获得的整景多光谱水体辐射信息图像lw(visn)′的水深计算。依据逐像元的蓝光波段和绿光波段辐射强度值，计算获得整景遥感图像的浅海水深信息。

实验结果如图7所示，较好地显示了工作区水下地形特征，表明本发明提取的浅海水深信息精度较高，具备实用价值。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均包含在本发明的保护范围之内。