一种浅海地形反演方法、计算机设备及存储介质

1.本发明属于海洋遥感与浅海地形领域，尤其涉及一种基于测高卫星数据(icesat-2)与高光谱遥感影像(sentinel-2)的浅海地形反演方法、计算机设备及存储介质。


背景技术：

2.海底地形尤其是浅海地形是海洋遥感测绘的主要目标之一，根据浅海地形绘制的水深地形图可以应用于航海及航运、渔业捕捞、海洋环境治理、海洋工程和海洋资源开发等领域。
3.近年来，随着卫星测高技术的蓬勃发展，其分辨率及精确度得到了较大提升，利用卫星测高数据，能够对大地水准面及重力异常、浅海地形进行精确反演，利用测高卫星数据与光学影像成为了浅海地形反演的主要研究方法之一。
4.例如yue ma等提出了“satellite-derived bathymetry using the icesat-2lidar and sentinel-2 imagery datasets”(期刊名：isprs journal of photogrammetry and remote sensing，2020年)，其中公开了利用dbscan算法通过密度来滤除噪声提取海面海底光子的方法，该方法通过调整阈值来分别处理白天和夜晚的光子数据，但在原始光子数据中，相当一部分海底光子的密度相差很大，所以其结果会丢失部分密度与设定阈值相差较大的海底光子数据，导致海底光子数据不足，进而导致反演出来的参数不准确。
5.又例如hsiao-jou hsu等提出了“a semi-empirical scheme for bathymetric mapping in shallow water by icesat-2and sentinel-2:a case study in the south china sea”(期刊名：isprs journal of photogrammetry and remote sensing，2021年)，其中公开了通过建立光学影像蓝绿波段对数比(r值)与激光测深之间的关系来反演浅海地形的方法，该方法在处理一个像素对应多光子的情况时，采用对相应的光子取平均深度的方法，该方法在处理过程中丢失了一些高精度光子信息，使得实际参与光子数量减少，会在一定程度上造成反演的参数不够准确。
6.由于激光测高卫星(icesat-2)的光子数据普遍存在较多噪声，需要人工去噪处理才能获得置信度较高的数据，且相比光学影像，测高数据的数据量较少，以icesat-2测高卫星的光子数据为例(icesat-2是目前激光测高卫星中卫星最高的)，icesat-2卫星用3对绿色激光(每一对激光由强弱两束激光组成，两者距离90米，每对激光相距3.3千米)沿卫星轨道方向测量其到地表的距离并解算出每个光子的地表高度，其数据相对稀疏。而高光谱影像(如sentinel-2影像)的分辨率可达10m，包含的信息更全更广，如何获得高精度的测高卫星数据以及如何充分利用高精度光子数据是浅海地形反演的关键。


技术实现要素：

7.本发明的目的在于提供一种浅海地形反演方法、计算机设备及存储介质，以解决现有技术中反演出来的参数不准确的问题。本发明通过滑窗分段，先采用“双峰高斯拟合+
均值滤波”方式提取海面光子及海面以下光子，再通过“单峰高斯拟合+中值滤波”方式提取海面以下光子数据中的海底光子，即使海底光子密度不一，也能较完整地提取出海底光子，保证海底光子数据量；通过双线性插值法，为每个光子找到一个对应的r值，采用所有高精度光子数据来参与水深反演，保证高精度光子数据不丢失，保证反演参数的准确性。
8.本发明是通过如下的技术方案来解决上述技术问题的：一种浅海地形反演方法，包括以下步骤：
9.获取高光谱影像，计算所述高光谱影像中每个像素的r值；
10.获取光子数据，并对所述光子数据进行清洗，去除不在研究区域内的光子数据、非海域数据以及噪声数据；
11.对清洗后的光子数据进行分段处理，得到多段光子数据；
12.计算每段光子数据的直方图，每段光子数据对应一个直方图；
13.采用双峰高斯函数对每个所述直方图进行拟合，得到与每个所述直方图对应的最优双峰高斯函数；
14.根据最优双峰高斯函数的参数，采用高斯均值滤波法从对应的直方图中滤除海面以上光子数据，提取出海面光子数据以及海面以下光子数据；
15.采用单峰高斯函数对海面以下光子数据进行拟合，得到最优单峰高斯函数；
16.根据最优单峰高斯函数的参数，采用中值滤波法从海面以下光子数据中提取出海底光子数据；
17.对所述海底光子数据进行折射校正，得到折射校正后的海底光子数据；
18.按照经纬度将折射校正后的海底光子数据与所述高光谱影像中每个像素的r值对齐；
19.根据每个海底光子的相邻像素的r值，计算每个海底光子的r值；
20.对所有海底光子的r值和水深进行拟合，确定拟合函数的系数，所述拟合函数的系数即为反演参数。
21.优选地，采用滑窗算法对清洗后的光子数据进行分段处理，得到多段光子数据。
22.进一步地，所述计算每段光子数据的直方图的具体实现过程为：
23.根据每段光子数据内所有光子的高度范围以及光子数目计算垂直分段高度；
24.在该段光子数据的高度方向上计算该段光子数据中每个垂直分段高度内的光子数目；
25.根据所述垂直分段高度以及每个所述垂直分段高度内的光子数目形成该段光子数据的直方图。
26.进一步地，所述垂直分段高度的计算表达式为：
[0027][0028]
其中，hi为第i段光子数据的垂直分段高度，hi为第i段光子数据内所有光子的高度范围，ni为第i段光子数据内的光子数目。
[0029]
进一步地，所述双峰高斯函数的表达式为：
[0030]
[0031]
其中，f(x)1表示双峰高斯函数，a1、a2均表示振幅，u1、u2均表示期望，u1＞u2，σ1、σ2均表示标准偏差；当x表示光子高度时，最优双峰高斯函数f(x)1的振幅a1近似为对应段光子数据中落入垂直分段高度内的最大海面光子数，振幅a2近似为对应段光子数据中落入垂直分段高度内的最大海底光子数，期望u1为对应段光子数据中的近似海面高度，期望u2为对应段光子数据中的近似海底高度，σ1为海面光子波峰的标准偏差，σ2为海底光子波峰的标准偏差。
[0032]
优选地，当光子高度x＞u1+3σ1时为海面以上光子；
[0033]
当u
1-3σ1≤光子高度x≤u1+3σ1时为海面光子；
[0034]
当光子高度x＜u
1-3σ1时为海面以下光子；其中u1为最优双峰高斯函数的海面光子波峰的期望，σ1为最优双峰高斯函数的海面光子波峰的标准偏差。
[0035]
进一步地，所述单峰高斯函数的表达式为：
[0036][0037]
其中，f(x)2表示单峰高斯函数，a3表示振幅，u3表示期望，σ3表示标准偏差；当x表示光子高度时，最优单峰高斯函数f(x)2的振幅a3近似为对应段光子数据中落入垂直分段高度内的最大海底光子数，期望u3为对应段光子数据中的近似海底高度，σ3为海底光子波峰的标准偏差。
[0038]
优选地，当m
j-3σ3≤光子高度x≤mj+3σ3时为海底光子，其中σ3为最优单双峰高斯函数的标准偏差，mj为第j次中值滤波时的中值。
[0039]
进一步地，所述折射校正的具体表达式为：
[0040][0041]
其中，z为折射校正后的海底光子水深，z0为折射校正前的海底光子水深，na为激光波束在空气中的折射率，nw为激光波束在海水中的折射率。
[0042]
进一步地，采用双线性插值法计算每个海底光子的r值，具体计算表达式为：
[0043][0044][0045][0046]
其中，(x
p
,y
p
)为海底光子p的经纬度，(x1,y1)为海底光子p左下角的相邻像素的经纬度，(x1,y2)为海底光子p左上角的相邻像素的经纬度，(x2,y1)为海底光子p右下角的相邻像素的经纬度，(x2,y2)为海底光子p右上角的相邻像素的经纬度，r
11
为像素(x1,y1)的r值，r
21
为像素(x2,y1)的r值，r
12
为像素(x1,y2)的r值，r
22
为像素(x2,y2)的r值，r
p
为海底光子p的r值。
[0047]
进一步地，所述方法还包括：根据所述反演参数以及高光谱影像中每个像素的r值计算每个像素的水深，再根据每个像素的水深以及海面高程获取浅海地形。
[0048]
本发明还提供一种计算机设备，包括：存储器，用于存储计算机程序；处理器，用于
执行所述计算机程序时实现如上所述浅海地形反演方法的步骤。
[0049]
本发明还提供一种存储介质，所述存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如上所述浅海地形反演方法的步骤。
[0050]
有益效果
[0051]
与现有技术相比，本发明的优点在于：
[0052]
本发明提供的一种浅海地形反演方法、计算机设备及存储介质，采用滑窗算法对测高卫星的带状光子数据进行分段，再对分段后的光子数据进行精细处理，去除噪声，获得高精度的海面光子和海底光子，然后将海底光子水深与对应时间段的高光谱影像像素的蓝绿波段对数比(即r值)以二次函数模型建立关系，计算反演参数，最后以此反演参数对整个高光谱影像进行反演来获得高精度的浅海地形；
[0053]
本发明采用“双峰高斯拟合+均值滤波”方式稳健、准确地区分出海面以上光子(噪声数据)、海面光子、海面以下光子(包括海底光子和噪声数据)，再采用“单峰高斯拟合+中值滤波”方式提取海底光子，即使海底光子密度不一，也能较完整地提取出海底光子，保证了海底光子数据量；针对一个像素对应多个光子数据的情况，采用双线性插值法代替均值法来获取所有海底光子信息，以便充分利用所有海底光子的高精度信息反演出高精度浅海地形。
[0054]
本发明针对浅海地形反演，较其他现有反演方法具有更高的反演精度。
附图说明
[0055]
为了更清楚地说明本发明的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一个实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0056]
图1是本发明实施例中浅海地形反演方法流程图；
[0057]
图2是本发明实施例中筛选云量后经大气校正的高光谱影像；
[0058]
图3是本发明实施例中滑窗分段以及垂直分段示意图；
[0059]
图4是本发明实施例中双峰高斯函数拟合示意图；
[0060]
图5是本发明实施例中双线性插值示意图；
[0061]
图6是本发明实施例中采用不同方法拟合反演参数示意图，其中图(a)为采用均值法进行反演，参与反演光子数为769，拟合度为0.85；图(b)为采用双线性插值法进行反演，参与反演光子数为5771，拟合度为0.86；
[0062]
图7是本发明实施例中利用反演参数获得的浅海水深图。
具体实施方式
[0063]
下面结合本发明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0064]
下面以具体地实施例对本技术的技术方案进行详细说明。下面这几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例不再赘述。
[0065]
如图1所示，本实施例所提供的一种浅海地形反演方法，包括以下步骤：
[0066]
1、高光谱影像处理
[0067]
(1.1)高光谱影像获取：高光谱影像(sentinel-2的影像)可以通过公布的影像数据获得(https://scihub.copernicus.eu/dhus/#/home)。
[0068]
(1.2)高光谱影像筛选：浅海地形反演需要光学影像浅海区域像素的反射率，为了确保影像中浅海区域没有被云层遮挡且成像清晰，选择云量少于10％，尤其是浅海区域无云或少云，且光照充足的影像，尽量避免选择凌晨或傍晚时刻的影像，因为此时太阳高度角较小，阳光斜射海面，反射到卫星传感器的能量较少。
[0069]
高光谱影像筛选保证了所选影像的浅海区域没有云层遮挡且光照充足，使得浅海区域能完整成像，以便于后续利用高光谱影像浅海区域的像素的波段反射值进行反演计算。
[0070]
(1.3)大气校正：由于公布的sentinel-2数据是没有经过大气校正的，需要通过官方的sentinel-2l2a产品生成和格式化处理软件sen2cor(http://step.esa.int/main/snap-supported-plugins/sen2cor/)进行大气校正，以消除大气对遥感信号传输过程的影响，图2所示为筛选云量后经过大气校正的高光谱影像。
[0071]
(1.4)高光谱影像中每个像素的r值的计算：根据高光谱影像中每个像素蓝绿波段值计算每个像素的r值。
[0072]
2、光子数据处理
[0073]
(2.1)光子数据读取：从nsidc(national snow&ice data center，美国国家冰雪数据中心)(https://nsidc.org/data/atl03)下载的光子文件中提取出光子数据，光子数据包括光子所需关键信息，例如经纬度、高度、沿轨距离，组为一个多维列表。
[0074]
(2.2)光子数据清洗：按照经纬度和高度进行光子数据清洗，去除不在研究区域内的光子数据、高度过高和过低的非海域数据以及噪声数据。不在研究区域内的光子数据对反演过程无用，应剔除，而光子高度过高或过低的噪声数据会影响光子提取过程，也应剔除，数据清洗提高了光子提取效率和提取精度。
[0075]
(2.3)光子数据分段：采用滑窗算法对清洗后的光子数据进行分段处理，得到多段光子数据。
[0076]
如图3所示，采用滑窗在水平方向上对带状的光子数据进行分段，得到多段光子数据，每段光子数据对应一滑窗分段。对分段后的光子数据能够进行更为精细的处理，有利于噪声去除，有利于获取高精度的海面光子和海底光子。
[0077]
(2.4)计算每段光子数据的直方图，具体实现过程为：
[0078]
(2.41)根据每段光子数据内所有光子的高度范围以及光子数目计算垂直分段高度，具体计算公式为：
[0079][0080]
其中，hi为第i段光子数据的垂直分段高度，hi为第i段光子数据内所有光子的高度范围，ni为第i段光子数据内的光子数目。
[0081]
(2.42)在该段光子数据的高度方向上计算该段光子数据中每个垂直分段高度内的光子数目。如图3所示，对某段滑窗分段(即某段光子数据)，在高度方向上按照垂直分段
高度进行垂直分段，计算每个垂直分段内的光子数目。
[0082]
(2.43)根据垂直分段高度以及每个垂直分段高度内的光子数目形成该段光子数据的直方图，即根据该段光子数据中每个垂直分段内的光子数目形成该段光子数据的直方图，如图4所示。
[0083]
(2.5)对于每个直方图，采用双峰高斯函数进行拟合：采用双峰高斯函数对每个直方图进行最小二乘拟合，选择残差平方和最小的双峰高斯函数及其参数作为对应的直方图的最优双峰高斯函数，即最优拟合函数，双峰高斯函数的表达式为：
[0084][0085]
其中，f(x)1表示双峰高斯函数，a1、a2均表示振幅，u1、u2均表示期望，σ1、σ2均表示标准偏差；当x表示光子高度时，最优双峰高斯函数f(x)1的振幅a1近似为对应段光子数据中落入垂直分段高度内的最大海面光子数，振幅a2近似为对应段光子数据中落入垂直分段高度内的最大海底光子数，期望u1为对应段光子数据中的近似海面高度，期望u2为对应段光子数据中的近似海底高度，σ1为海面光子波峰的标准偏差(与海面光子的垂直范围有关)，σ2为海底光子波峰的标准偏差(与海底光子的垂直范围有关)。如图4所示，期望值大的波峰为海面光子波峰，如图4中右边的波峰，期望值小的波峰为海底光子波峰，如图4中左边的波峰。
[0086]
(2.6)采用高斯均值滤波法提取海面光子以及海面以下光子
[0087]
当光子高度x＞u1+3σ1时为海面以上光子；当u
1-3σ1≤光子高度x≤u1+3σ1时为海面光子；当光子高度x＜u
1-3σ1时为海面以下光子。
[0088]
经过海面光子提取后，剩下的海面以下光子包括海底光子和噪声，再采用单峰高斯函数提取精确的海底光子。
[0089]
(2.7)采用单峰高斯函数进行拟合：采用单峰高斯函数对海面以下光子数据进行最小二乘拟合，选择残差平方和最小的单峰高斯函数及其参数作为最优单峰高斯函数，即最优拟合函数，单峰高斯函数的表达式为：
[0090][0091]
其中，f(x)2表示单峰高斯函数，a3表示振幅，u3表示期望，σ3表示标准偏差；当x表示光子高度时，最优单峰高斯函数f(x)2的振幅a3近似为对应段光子数据中落入垂直分段高度内的最大海底光子数，期望u3为对应段光子数据中的近似海底高度，σ3为海底光子波峰的标准偏差(与海底光子的垂直范围有关)。海底通常含有较多噪声，采用单峰高斯进行拟合更准确。
[0092]
(2.8)采用三次中值滤波法提取海底光子
[0093]
当m
j-3σ3≤光子高度x≤mj+3σ3时为海底光子，其中mj为第j次中值滤波时的中值，本实施例中，j＝3，仅需较少次数的中值滤波就可以达到很好的去噪效果，得到高精度、高可靠性的海底光子。。
[0094]
(2.9)海底光子的折射校正，具体计算公式为：
[0095][0096]
其中，z为折射校正后的海底光子水深，z0为折射校正前的海底光子水深，na为激光
波束在空气中的折射率，nw为激光波束在海水中的折射率。光子数据包括光子的经纬度、高度、沿轨距离等，海底光子水深等于海面光子高度减去海底光子高度。
[0097]
3、数据对齐：按照经纬度将折射校正后的海底光子数据与高光谱影像中每个像素的r值对齐。
[0098]
4、计算海底光子的r值：根据每个海底光子的相邻像素的r值，采用双线性插值法计算每个海底光子的r值。
[0099]
高光谱影像中每个像素的r值为格网结构，海底光子为一个三维点，在数据对齐时，由于海底光子的经纬度与格网中心(即高光谱影像的像素中心)的经纬度不会完全一致，所以海底光子落入每个格网结构中后是散落在格网中心附近的，为了计算每个海底光子对应的r值，通过双线性插值法在每个海底光子的位置，利用其相邻的四个格网中心的r值计算出海底光子对应的r值，则n个海底光子经过双线性插值计算，可得n个相应的r值。
[0100]
如图5所示，设某个海底光子p在经纬度坐标系中的位置为(x
p
,y
p
)，与海底光子p相邻的四个格网中心或像素中心分别为(x1,y1)、(x1,y2)、(x2,y1)、(x2,y2)，四个网格中心对应的r值分别为r
11
、r
12
、r
21
、r
22
，对海底光子p所在位置进行双线性插值，以获取海底光子p所在位置对应的r值。
[0101]
(4.1)先在x方向对点m(x
p
,y1)、点n(x
p
,y2)进行线性插值(m、n两点分别为直线x＝x
p
与直线y＝y1、y＝y2的交点)：
[0102][0103][0104]
其中，rm为点m对应的r值，rn为点n对应的r值。
[0105]
(4.2)再在y方向对点p(x
p
,y
p
)进行线性插值，公式为：
[0106][0107]
其中，r
p
为海底光子p的r值(即图6中的横坐标r＇)。
[0108]
图6为采用不同方法拟合反演参数示意图，其中图(a)为采用均值法进行反演，参与反演光子数为769，拟合度为0.85；图(b)为采用双线性插值法进行反演，参与反演光子数为5771，拟合度为0.86。从图6中可以得出，图(b)参与反演的光子数更多，反演出的参数更准确，可信度更高。
[0109]
5、参数反演：对所有海底光子的r值和水深进行二次函数拟合，确定拟合函数的一次项系数、二次项系数以及常数项，一次项系数、二次项系数以及常数项即为反演参数。
[0110]
二次函数模型为：
[0111]zp
＝a(r
p
)2+br
p
+c
ꢀꢀꢀꢀꢀꢀ
(8)
[0112]
其中，z
p
为海底光子p的水深(即图6中的纵坐标z)，a、b和c分别为二次函数模型的二次项系数、一次项系数和常数项，即a、b和c为待求的反演参数。示例性的，如图6(a)中，a为6386.188，b为12553.308，c为6169.736；如图6(b)中，a为6742.714，b为13281.673，c为6541.305。
[0113]
6、浅海高程模型：根据反演参数以及高光谱影像中每个像素的r值计算每个像素
的水深，再根据每个像素的水深以及海面高程获取浅海地形，以及根据每个像素的水深获取浅海水深图。
[0114]
由步骤5可以求得反演参数，式(8)中的二次函数系数已知，再根据已知的每个像素的r值即可计算出每个像素的水深，像素的水深减去海面高程等于海底高程，即可获得整个研究区域的数字高程模型以及浅海水深图，如图7所示浅海水深图。
[0115]
本发明实施例还提出了一种计算机设备，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上所述浅海地形反演方法方法的步骤。
[0116]
示例性的，所述计算机程序可以被分割成一个或多个模块/单元，所述一个或者多个模块/单元被存储在所述存储器中，并由所述处理器执行，以完成本发明。所述一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述所述计算机程序在所述计算机设备中的执行过程。
[0117]
所述设备可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述设备可包括但不仅限于处理器、存储器。本领域技术人员可以理解，本发明计算机设备仅仅是设备的示例，并不构成对设备的限定，可以包括比设备更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述设备还可以包括输入输出设备、网络接入设备、总线等。
[0118]
所述处理器可以是中央处理单元(central processing unit，cpu)，还可以是其他通用处理器、数字信号处理器(digital signal processor，dsp)、专用集成电路(application specific integrated circuit，asic)、现成可编程门阵列(field-programmable gate array，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
[0119]
所述存储器可用于存储所述计算机程序和/或模块，所述处理器通过运行或执行存储在所述存储器内的计算机程序和/或模块，以及调用存储在存储器内的数据，实现所述浅海地形反演方法的每个步骤。所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序(比如声音播放功能、图像播放功能等)等；存储数据区可存储根据手机的使用所创建的数据(比如音频数据、电话本等)等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(smart media card，smc)，安全数字(secure digital，sd)卡，闪存卡(flash card)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。
[0120]
所述计算机程序被处理器执行时实现所述浅海地形反演方法的步骤。
[0121]
所述计算机设备集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、u盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(rom，read-only memory)、随机存取存储
器(ram，random access memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。
[0122]
以上所揭露的仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或变型，都应涵盖在本发明的保护范围之内。