浅海遥感水深反演研究及其软件设计与实现

2023-12-01廖旋芝熊显名

软件导刊 2023年11期

廖旋芝，邱骞，熊显名

（1.北部湾大学电子与信息工程学院，广西钦州 535011；2.广西幼儿师范高等专科学校初等教育学院，广西南宁 530022；3.桂林电子科技大学广西光电信息处理重点实验室，广西桂林 541000）

0 引言

海洋在地球表面的覆盖面积达到71%，浅海又是人类与海洋接触最频繁的海域，因此无论从人类发展出发，还是从军事、国防角度来看，人们都期望利用海洋资源来解决能源短缺问题。浅海水下信息探测则是海洋资源探测中的一项重要内容［1-3］。不同于传统船载和机载探测水深信息方法的高成本、地域限制、不适用大规模探测等特点［4-6］，卫星遥感技术具有不受地理空间限制、大范围实时同步，以及全天时、全天候多波段成像等优势［7-8］。高分辨率的遥感影像具有空间与时间分辨率高、信息量大的特点，可为浅海水深信息探测提供有效的探测手段［9］。

针对浅海遥感水深反演理论研究方面，在国外，Bierwirth 等［10］提出一种计算水体底部反射率的方法，分析卫星图像获取到的辐射亮度与水深及对应区域海底反射率的关系，得出参数关系，从而求出水深；Lyzenga 等［11］提出利用主成分分析提取水深信息的方法，建立卫星遥感影像光谱值与对应海域水深值之间的关系模型，计算出水深数据；Stumpf 等［12］利用波段比值与水深的线型回归得出水深，在一定程度上消除反演区域不同海底地质的影响。在国内，党福星等［13］依据海底地质的光谱特性建立多光谱卫星影像在海岛浅海的水深反演模型，获得了较好结果；程洁等［14］采用国产高分卫星2 号多光谱影像对香港平洲岛进行水深反演，反演模型适合中等浅水区域，取得了较好的反演效果；李宁宁等［15］在对岛礁周边水深进行反演模型分析时，得出对数波段比值模型具有效率优势；楚森森等［16］基于光谱分层研究浅海水深，解决了水深反演模型难以适应不同光谱层的局限，提高了水深反演精度。

针对软件应用方面，近年来，国内外水深遥感测量方面的研究越来越多，但能够进行浅海水深遥感反演的软件很少。已有的软件主要分为以下几种类型：基于高级语言VC 的图像处理系统、基于MATLAB 软件的遥感模型和仿真系统，以及基于IDL 或ENVI 等可视化数据语言的遥感影像处理系统。但大多数遥感软件价格昂贵，并且仅仅作为遥感影像处理的通用平台。IDL（Interactive Data Language）作为第四代编程语言，具有强大的数据与图像处理、可视化以及交互式处理和操作的功能。其提供了大量函数和对象，可以更直观、高效地处理遥感数据，因此得到了越来越多遥感影像处理用户的青睐［17］。软件系统依托交互式数据语言IDL 平台开发，可实现对浅海海域水深反演的要求。本项目的研究将为浅海遥感影像水深探测应用提供快速、高效的水深信息可视化处理软件平台。

综上所述，目前一般利用卫星遥感影像进行浅海水深反演，其卫星影像受大气和传感器自身精度影响，在反演时需进行繁琐的参数纠正，增加了反演难度。同时，实测水深数据无法做到与反演水深数据大量匹配，无法进行多区域的模型验证，反演模型存在一定的区域适用性。本研究将以卫星图像获取的DN 值为基础构建水深反演模型，在以IDL 为编程语言的软件平台上解析影像头文件参数，自动配置纠正参数，而无需过多的繁琐操作。在具备少量的实测水深数据时，能依据反演模型进行绝对水深反演。同时，在缺乏实测水深数据时能进行相对水深反演，获取浅海水下的相对地势变化。软件平台的开发将缩短浅海水深反演周期，降低使用人员的学习成本，以快速、高效地进行浅海水深反演。

1 遥感水深反演理论分析

1.1 水深遥感理论

水深遥感探测的原理主要基于太阳光在水体的反射过程，探测器接收到含有水深信息的反射光能量，其辐射传输过程如图1 所示（彩图扫OSID 码可见，下同）。太阳光线照射到水面，其中一小部分光线被水面直接反射到大气中，形成水面散射光；未被反射的光线进入水体，以折射和透射的方式在水中传播。其中部分光线被水分子吸收和散射，部分被水体中含有的其他悬浮物质、水体微生物、水下动植物等吸收、反射、散射和衍射，形成后向散射光。光线到达水体底部被反射，并且由水体折射进大气中，探测器接受到的辐射亮度与水深信息有关。而且随着水体所含的杂质增多，浑浊度增加，水深与水底反射的光强呈负相关关系，即水体越浑浊，探测器接收到的含有水深信息的辐射亮度越低。因此，遥感测深包括太阳光在大气和水体中的两个辐射传输过程。卫星传感器所接收的信号主要包括水—大气界面的反射信号、水体的后向散射信号以及水底的反射信号，其中水底的反射信号具有水深信息。

Fig.1 Radiative transfer process of sunlight in the atmosphere and water bodies图1 太阳光在大气与水体中的辐射传输过程

1.2 水深光谱特性

不同的水体由于所含物质（叶绿素浓度、含泥沙量、浮游生物等）不同，在可见光波段范围内对水体的衰减程度也不同，如图2 所示。通常水质越混浊，衰减系数越大，即水质越清，穿透性越好。从图2 中还可得知对于同一水体，不同波段的衰减程度不同。但对于不同水体，衰减系数的极小值都出现在0.45～0.60 μm 左右的波段，说明蓝绿波段对水体的穿透性最好。但随着泥沙浓度的增大，水体的后向散射信号会增强，即出现“红移现象”。若水质为均一状态时与水深正线性相关，则混浊度较大的沿岸水一般呈黄绿色；延伸至深海，水则呈蓝绿色［18］。

Fig.2 Spectral attenuation characteristics of water图2 水的光谱衰减特性

1.3 水深遥感模型

可见光波段在水体中的传播传输过程可用经典的辐射传输方程来描述［19］，则该方程可表达为：

其中，Z为水深，θ为海底反射光到海面的入射角，φ为太阳光入射到海水时的折射角，L为离水辐射亮度。

假设辐射场具有各向同性性质，所有水体的光学参数只是水深Z的函数，大颗粒散射中的前向散射系数远远大于后向散射系数［19］，则公式（1）可转换为［20］：

式中，Ed(Z)为水中向下的辐亮度，ɑ为吸收系数，bb为背向散射系数。此公式说明向下辐亮度Ed(Z)随着深度Z按指数衰减，衰减系数是ɑ+bb。一般背向散射系数bb较小，所以Ed(Z)近似按e-ɑZ的规律衰减。当水体足够清澈，底质较均一，大气条件较好时，可根据遥感影像灰度值与可见光水体衰减系数之间的线性关系，依据式（2）推算出水深值。

（1）单波段模型。假设光进入水体后衰减系数和底质反射率为常量，其理论基础是布格定理［21］，如式（3）所示［22］。

式中，Lsi为深水区辐亮度，τi为与太阳辐亮度、大气和水面透过率、水面折射率等有关的常数，γBi为水体底质反射率，ki为水体衰减系数，f为水体路径长度［18］，Z为水深［23］。式（3）通过变形可得：

式中，Xi=ln(Li-Lsi)，令ɑ=-1/fki，b=(lnτiγBi)/fki。

（2）双波段比值模型。水体的衰减系数和底质反射率随着水体与底质类型的不同表现出较大差异，但不同波段的底质反射率之比却保持不变，且两波段衰减系数的差值也基本保持不变［24］。由式（3）可得第1 个波段的辐亮度L1和第2个波段的辐亮度L2：

由于不同波段的底质反射率之比保持不变，即γA1/γA2=γB1/γB2…，且对于不同水体，两波段衰减系数差值(k1-k2)基本保持不变。因此，有：

其中，ɑ=-1/(k1-k2)f，b=(lnτ1γB1/lnτ2γB2)/f(k1-k2)。

（3）线性多波段模型。线性多波段模型是由双波段模型扩展到n 个不同波段和n 种不同底质类型而得到的［21］。

式中，Xi=ln(Li-Lsi)，ω1、ω2为两个波段的权重因子，且满足ω1+ω2=1。将双波段模型中的波段扩展至n 个不同波段，底质类型扩展至n种不同类型，则有：

其中，有：

依据式（10）最后得到：

其中，ɑ0，ɑ1，······ɑn为待定系数。在实际计算中，选择n+1 个点，由实测获取n+1 组Z 值，回归计算出ɑ0，ɑ1，······ɑn的值。

上述模型在水深反演的应用中，采用卫星遥感影像的DN 值，结合卫星自身的纠正参数，可快速获取反演水深。同时，经过前人学者的研究和数据验证，上述模型的精确度也普遍较高，能够达到对浅海遥感水深反演软件快速、准确性的要求。

2 浅海遥感水深反演软件设计

2.1 软件总体设计

本文研究的浅海遥感水深反演软件系统作为水深反演的专用平台，主要包含四大基本功能模块：输入显示、水深反演、输出保存、其他辅助功能模块。每个基本功能模块可分别实现相应功能，同时模块之间相互联系，后者建立在前者之上。模块的基本功能如下：

（1）输入显示模块。系统可以灵活地载入多种卫星遥感影像数据，准确地提取影像头文件中所包含的地理信息、坐标系、分辨率、投影方式等信息，并将影像默认为1：1的显示，也可对影像进行缩放和漫游，以便用户查看。

（2）水深反演模块。由于遥感影像获取途径特殊，为保证反演结果的精度，在参与反演之前，首先要经过预处理。不同影像需要做的预处理工作可能不同，但大都要经过辐射校正、几何校正、图像增强和图像融合。此处理一般由卫星数据服务商完成，也可获取相关参数，在ENVI 软件中完成。在该软件中，经过预处理后的影像数据还需要进行水陆分割、模型参数获取、水深计算和等深线绘制等处理，最终获得影像的水深反演结果。

（3）输出保存模块。通过对水深反演生成的等深线进行编辑修正，获取最终的水深反演专题图，为海上作业提供水深参考。该模块包含等深线编辑、水深反演图输出、数据保存等。

（4）其他辅助功能模块。为了使软件功能更加丰富、实用，软件还包含一些辅助工具，如深度剖面分析、水深点提取、影像裁剪等。

根据软件的总体布局及其功能分析，浅海遥感水深反演软件系统组成结构如图3所示。

Fig.3 Software system composition structure图3 软件系统组成结构

2.2 软件系统功能设计

2.2.1 人机交互界面设计

浅海水深反演软件是一个能连接各大功能模块的完整系统，软件的主界面是用户与计算机交流的主要平台，直接地体现了软件设计水平。因此，简洁、清晰、合理的布局是主界面设计的基础，可使人机交互简单、方便。本系统主界面设计采用基于过程的GUI 设计方法，创建基本容器（base 组件），然后在该容器内创建其他基本容器组件，从而形成一个完整、自上而下的顶层结构。

2.2.2 输入显示模块设计

输入显示模块的功能主要是将数据输入软件系统，之后进行影像显示方面的参数设置。主要包括是否显示经纬度网、使用何种定标方式、确定图像显示的拉伸系数等，最后在图像显示区将图像显示出来。在本软件中，可输入的软件数据分为以下5 种：高分辨率的卫星遥感影像数据（*.tif 格式）、水陆分割数据（*.dat 格式）、参数标定数据（*，txt）、水深数据（*.tif 格式）和矢量等深线数据（*.shp 格式）。在载入影像时，自动解析影像头文件，如xx_metadata.txt 文件，解析获得各波段的定标系数。调用IDL 具备的大气校正函数库，消除大气散射、吸收、反射对卫星实际接受到的辐射亮度的影响。软件自动进行上述操作。

2.2.3 水深反演模块设计

水深反演模块是本系统的核心部分，包括水陆分割、模型参数获取、水深计算和等深线绘制四大功能。

（1）水陆分割。该功能可将影像分成两部分，实现水体白色和陆地黑色的二值化图，并且在进行水深计算时会将陆地部分剔除，降低计算量，提高程序的运行速度。

（2）模型参数获取。在有实测数据的情况下，选取适合本研究区域的水深反演理论模型，结合实测点数据，读取水深反演因子的灰度值，并进行最小二乘法回归计算，以确定反演模型参数。之后建立水深反演模型，把该模型代入到其他水深点可以得到相应的实际水深值。

（3）水深计算。除采用双波段比值模型外，本软件拟设计模型包含单波段模型、双波段减法比值模型和多波段自定义模型，这些模型的设计可为软件的数据处理提供多样化的选择。

（4）等深线绘制。采用矩形网格法绘制等深线，先将离散数据转换成网格数据，然后计算各条等值线与网格边的交点，按一定规律追踪并有序排列所有等值点，最后平滑绘制等值线。

2.2.4 输出保存模块设计

输出保存模块包括等深线编辑、反演图输出和数据保存功能。其中，等深线编辑主要对生成的等深线进行修正，以达到最佳的水深反演效果；反演图输出包含等深线图输出、水深灰度图输出；数据保存包含水陆分割图保存、等深线保存和水深反演数据保存。

2.2.5 其他辅助模块设计

软件还包含一些辅助工具，如深度剖面分析、水深点提取、影像裁剪等。深度剖面图分析的开发是为了让用户更好地分析感兴趣的水域，从而了解水域的深度走势；水深点提取即用户输入兴趣点的经纬度坐标，获取该点的水深值，方便评测人员将反演水深值与海图或实测数据进行比较；影像剪裁不仅可减少计算量，提升软件运行速度，而且可让感兴趣区域的显示更加醒目，方便用户观察。

3 浅海遥感水深反演软件实现

3.1 IDL开发语言选择

交互式数据语言IDL 是美国Exelis VIS 公司开发并投向市场的用于数据可视化研究与应用开发的第四代计算机语言。IDL 语言的最大特点是面向矩阵，对于数组的操作可以不通过循环而直接对数组进行运算。此外，基于IDL 语言开发的遥感图像处理平台ENVI 与IDL 语言能够方便地进行数据交互和函数调用，极大地提高了IDL 语言处理遥感图像的能力［17，25］。因此，浅海遥感水深反演软件采用IDL 语言开发，可方便调用图像处理函数库，且可视化的开发环境能大大缩短开发周期。

3.2 软件系统功能实现

3.2.1 人机交互界面实现

软件主界面分为软件菜单栏、影像属性栏、地理信息栏、系统状态显示栏以及影像显示区、鹰眼显示区等部分。软件菜单栏负责功能界面的激活，包括影像载入、参数设置、水陆分割、反演参数获取、水深反演及各种辅助功能等；影像属性栏用于显示被加载影像的信息，包括影像路径、影像的波段数和卫星影像基本信息；地理信息栏用于显示影像的地理信息，包括投影方式、影像4 个顶点和中心的经纬度、分辨率和坐标系；系统状态显示栏用于记录当前系统的运行状态，显示当前鼠标位置的经纬度、DN 值和水深值；影像显示区用于显示当前加载的卫星影像；鹰眼显示区用于存放影像的缩略图，通过鼠标点击确定鹰眼区图像，并在影像显示区显示影像的具体位置。软件系统界面如图4所示。

Fig.4 Software system interface图4 软件系统界面

3.2.2 输入显示模块实现

IDL 提供了直接输入和间接输入两种数据输入方式。直接输入是指在GUI 界面采用键盘输入的方式，该方式能满足常规的数据交互需求。间接输入是指读取格式文件的操作方式，该方式无需通过键盘输入，可直接读取内存中已有的数据。该方式能够以任意格式，用文件的形式进行大容量数据的输入，操作简洁、方便。在系统内，多数的数据输入都是应用该方法。影像显示包括影像属性显示、地理信息显示、影像鹰眼区域显示、影像原比例显示、彩色合成显示和系统运行状态栏显示。

3.2.3 水深反演模块实现

水深反演模块主要分为4 部分：水陆分割、模型参数获取、水深计算和等深线绘制。水陆分割是指把影像按水域和陆地分为两部分，生成二值图；模型参数获取是指通过影像的实测点数据，采用最小二乘法获得模型参数用于水深计算；水深计算是指在水陆分割的前提下，把影像陆地部分屏蔽，按照一定的算法模型和反演参数计算出水体深度；等深线绘制是指在已得到水深值的情况下，将水深值相等的点用等深线连接起来，便于理解和观察。重点说明以下关键技术：

（1）水陆分割。软件进行水陆分割的目的是在水深反演时将陆地信息屏蔽，从而减少计算量，降低内存占用以提高软件的运行速度。在软件的菜单栏中设置了“水陆分割”栏，单击后可生成两个下级菜单：RGBN 四波段水陆分割和近红外波段水陆分割。

（2）水深计算。软件为水深计算设计了单独的功能界面，用户可选择任意模型计算影像水深。选择的模型不同，操作界面也有所不同，包括双波段比值水深模型、单波段模型、双波段比值减法模型以及多波段自定义模型，使得软件具备更强的应用性。

现阶段主要采用双波段比值模型反演水深信息。图5为双波段比值模型界面，用户通过界面的下拉列表选取分子波段和分母波段。由前面介绍的遥感水深理论可知，蓝光和绿光在水中吸收和衰减较小，反演结果也证实了蓝绿光的反演效果较好。所以一般把蓝光作为分子波段，绿光作为分母波段计算水深值。在没有水域实际深度值的情况下，系统会给出默认值（a=-30，b=0）。按此计算出来的水深为相对水域深度，在显示上不影响等深线绘制和对水体深度变化走势的观察。

Fig.5 Interface of dual band ratio model parameter setting图5 双波段比值模型参数设置界面

（3）等深线绘制。采用基于格网序列法的矩形格网法生成等深线图，等深线的条数及每条线的深度值根据不同区域和用户的不同要求由用户输入。等深线绘制的功能界面如图6 所示。界面上方为水深统计直方图，中间标注了水深的范围，通过“编辑深度值”“添加等深线”“删除等深线”按钮可以对深度值信息进行编辑。在编辑等深线子界面中可以对等深线的值、标号、颜色进行设置。

Fig.6 Interface of contour drawing图6 等深线绘制界面

3.2.4 输出保存模块实现

经过等深线编辑的影像会生成SHP 文件格式的矢量数据，配上原始影像即可做成影像等深线图。在软件重新启动时，无需进行参数配置，此等深线图能被直接加载显示，可展示水深反演效果，方便阅览。

3.2.5 其他辅助模块实现

软件包含了一些方便用户操作和观察影像的辅助工具。

（1）深度剖面图分析。深度剖面图分析工具是为了观察直线水深段的具体水深变化而设计的，深度剖面图界面如图7 所示。图中的斜直紫线表示观察区域，在实际操作中，先点击紫线的右下方，然后点击左上方，最后会弹出深度剖面图功能界面。可以看出，水体深度曲线有明显的逐渐下滑趋势，与实际情况类似。

Fig.7 Interface of depth profile图7 深度剖面图界面

（2）水深点提取。用户输入测试点的经纬度，点击确认后，在不超过影像范围的情况下，深度值一栏返回测试点的水深值。

（3）影像裁剪。利用影像裁剪工具可以在保留影像地理信息的基础上提取目标海域。

4 软件数据处理

在有实测水深数据时，软件可以水深反演计算出绝对水深。在缺失实测水深数据时，软件可以计算出相对水深，以此反映水深变化。

4.1 反演区域数据说明

绝对水深反演区域为新竹市南寮附近海港，截取范围为24° 49′40.74″N-24° 51′49.42″N，120° 53′45.11″E-120°56′50.90″E 的浅海区域。采用IKONOS 卫星于2005 年6 月4 日拍摄的高分辨率多光谱遥感影像，多光谱共含有4 个波段，分别为0.45-0.53μm（蓝）、0.52-0.61μm（绿）、0.64-0.72μm（红）、0.76-0.86 μm（近红外），分辨率为4 m。全色谱段为450-900 nm，分辨率为1 m。本研究采用同时期获取的电子海图数据，以电子海图上的水深信息作为实测数据。在电子海图上提取水深样点45处，水深在31 m 以内。

相对水深反演区域为南沙群岛郑和群礁南部边缘的鸿庥岛，截取范围为10°10′19.39″N-10°11′1.89″-N，114°21′52.71″E-114°22′53.68″。采用GeoEye-1 卫星于2009 年4 月8 日拍摄的高分辨率多光谱遥感影像，多光谱有4 个波段，分别是0.45-0.51μm、0.52-0.58μm、0.655-0.69μm 和0.78-0.92μm，分辨率为1.64 m。全色谱段为450-900 nm，分辨率为0.41 m。因相对水深反演没有实测水深数据，本研究采用逐步逼近法，假定反演模型参数值，运用水深反演模型获取反演水深，在此基础上通过线性回归获取反演模型参数值，再通过同一水深反演模型二次获得反演水深。将前后两次反演水深进行对比分析，计算软件反演效果。

为进行潮汐校正，本研究通过卫星过境同期潮汐表获取该时刻的潮位值。卫星遥感影像获取的瞬时水深为实测水深与该时刻的潮位值之和。因用电子海图水深作为实测水深，则此瞬时水深为电子海图水深与该时刻潮位值相加后的水深。

4.2 数据处理

经过浅海区域水质的光学特性分析后，发现蓝绿波段比值模型适合当前水深反演。因此，在对反演区域进行绝对和相对水深反演时，均采用蓝绿波段比值模型，即Z=ɑ*ln(蓝/绿)+b。经软件计算，获取反演模型参数值后，在绝对水深反演下，所得公式为Z=-168.094*ln(蓝/绿)-1.074 14；在相对水深反演下，所得公式为Z=-53.178 8*ln(蓝/绿)-13.569 6。在绝对水深反演下，将电子海图水深数据与绝对水深反演数据进行拟合，并进行评价参数计算。评价参数包括决定系数R2、均方根误差RMSE（m）、平均相对误差MRE（%）和平均绝对误差MAE（m）。同理，在相对水深反演下，将两次回归前后的水深反演数据进行拟合，如图8、图9 所示，并进行评价参数计算，如表1所示。

Table 1 Inversion fitting and evaluation parameters表1 反演拟合及评价参数

Fig.8 Fitting of absolute water depth inversion data图8 绝对水深反演数据拟合

4.3 数据分析

由表1 可知，在两种情况下的水深反演，决定系数R2均在0.98 以上，拟合效果好。RMSE、MRE 和MAE 在绝对水深反演与相对水深反演下，分别为0.572 6m、6.659 3%、0.479 2m 和1.076 1m、5.496 2%、0.659 2m，反演效果良好。在有实测水深数据时，软件能较好地反演真实水深，误差偏差较小。在无实测水深数据时，软件能较好地反映水深相对变化，回归前后的数据拟合度好，误差变化小，说明软件的稳定性较好。在软件中绘制出等深线和深度剖面图，效果如图10、图11所示。