CN113739707A - 一种海冰信息提取方法 - Google Patents

Abstract

本发明涉及一种海冰信息提取方法，通过对获取的光子数据进行预处理，结合计算结果对海冰信息进行进一步分类，根据分类情况计算出海冰干舷高度值，并通过对现有积雪深度测量模型结果进行优化并依据其它海冰信息参数，结合阿基米德定律、线性回归函数等计算公式求得海冰厚度值。与现有技术相比，本发明在对卫星初始数据进行校正后，使用较为精确的数据进行了后续处理及筛选，进一步对数据进行了限定，提高了数据的利用率。

Description

一种海冰信息提取方法

技术领域

本发明涉及海冰信息监测领域，尤其是涉及一种海冰信息提取方法。

背景技术

极地区域的环境变化是对于全球气候变化的重要监测指标，海冰是极地最重要的自然现象之一，强烈地制约着海气热量交换，从而对北极气候乃至北半球气候起到调控作用。海冰边缘的位置及其移动、海冰密集度、海冰厚度、海冰类型是海冰动力学、热力学以及气候研究的重要参数。针对于极地海冰信息提取研究方法多样，传统的信息提取方法均存在测量范围及精度上的问题，目前应用前景较好的精度高且测量便捷的单光子激光高度计，这种激光高度计的使用能基本满足海冰信息提取的要求；

申请号CN202110277320.1的发明公开了一种基于单光子激光测高数据的海冰厚度估算方法，该发明使用激光测高卫星进行相关数据处理，提供了海冰信息的检测方法。然而地球大气会在一定程度上干扰卫星发射接收光信号，以及潮汐等地球自然现象会使得光子数据产生偏移，影响了该数据的精确性，导致无法准确反映海冰信息。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种海冰信息提取方法，解决了目前海冰信息精确性不高的问题。

本发明的目的可以通过以下技术方案来实现：

一种海冰信息提取方法，包括：

S1、使用多波束单光子激光测高卫星获取初始光子高度，对所述初始光子高度进行地球物理校正，得到校正光子高度，将所述校正光子高度减去现有平均海平面高度，得到相对光子高度。

S2、以第一设定长度中的所有光子作为一个光子单元，可得到连续的多个光子单元，对每个光子单元进行粗表面提取：

根据每个光子的相对光子高度，对每个光子单元建立第一高度直方图，在高度直方图中，每个光子的高度值形成一个光子柱，保留第一设定高度范围内的光子，得到第二高度直方图，提取所述第二高度直方图中的上下波峰的设定百分比作为第二设定高度范围，保留第二设定高度范围内的光子，得到第三高度直方图，选取出第二设定高度范围内的光子作为选中光子，计算选中光子的高度平均值：

若当前光子单元的高度平均值大于差值阈值，则舍去这个光子单元的前一半长度中的光子，以当前光子单元后一半长度的光子和下一个光子单元前一半长度的光子作为新的光子单元重新计算高度平均值。

若当前光子单元的高度平均值小于差值阈值，则将所述高度平均值作为当前光子单元的粗表面高度。

S3、对经过粗表面提取的每个光子单元数据进行精表面提取：

将当前光子单元中所有光子的相对光子高度减去粗表面提取中得到的粗表面高度，得到差值光子高度，保留第三设定高度范围内的光子，以第二设定长度中的固定数量的光子作为一个迷你光子单元，根据差值光子高度，建立第四高度直方图，提取所述第四高度直方图中每个迷你光子单元的波峰，并判断当前迷你光子单元的波峰是否在所述第四高度直方图的第一或最后一列，且光子总数大于光子阈值数：

若是，则计算差值光子高度的平均值与标准差，保留多倍标准差范围内的差值光子高度，以该数据建立第五高度直方图并进行高斯拟合，重新计算差值光子高度及展宽，将差值光子高度加上此前减去的粗表面高度，得到当前迷你光子单元的所有光子的精表面高度。

若否，则将该迷你光子单元的数据删除，对下一个迷你光子单元进行上述判断。

S4、根据检测到的光子数，结合光子偏差校正公式，得到光子实际到达多波束单光子激光测高卫星的探测器的时间。

S5、根据光子实际到达多波束单光子激光测高卫星的探测器的时间，计算每个迷你光子单元的表面光子率和背景光子率。

S6、根据季节、表面光子率、归一化的背景光子率、展宽以及太阳高度角，结合冰水判断规则对每一个光子单元中的所有迷你光子单元进行冰水分类，所述冰水分类类型中包含冰间水道。

S7、根据冰水分类结果，对每一个光子单元计算海冰干舷：

若一个光子单元中存在一个冰间水道，其中每个冰间水道构成一个冰间水道单元。

若一个光子单元存在多个冰间水道，则其中连续的多个冰间水道可构成一个冰间水道单元，以每个冰间水道单元内冰间水道的长度为权重，计算所述冰间水道单元精表面高度的加权平均值并设为其高度，再对多个冰间水道单元依据冰间水道的长度进行加权平均计算整个光子单元的高度值作为该光子单元的海冰干舷高度。

若一个光子单元中不存在冰间水道，则先判断相邻的存在冰间水道的光子单元是否有效，不存在冰间水道的光子单元若处于有效的光子单元之间，且这两个有效的光子单元中光子实际到达探测器的平均时间的时间间隔小于时间阈值以及海冰干舷高度差值绝对值小于高度差阈值，则可通过线性插值的方法计算中间光子单元的海冰干舷高度。

判断光子单元是否有效的方式为计算该光子单元到陆地的距离，当所述距离小于距离阈值或者海冰密集度小于密集度阈值时，则将当前光子单元设为无效的光子单元，反之为有效的高度单元。

S8、使用Warren Snow气候学模型计算得到积雪深度，并以其计算值的固定倍数作为初始积雪深度值；。

S9、根据所述初始积雪深度值和所述海冰干舷高度值，使用线性回归函数计算得到积雪深度截断值和海冰干舷高度截断值。

S10、根据所述积雪深度截断值和海冰干舷高度值，计算出积雪深度重建值。

S11、根据所述积雪深度重建值、所述海冰干舷高度值、水的密度值、海冰密度值和积雪密度值，结合阿基米德浮力定律计算出海冰厚度。

进一步地，所述地球物理校正公式为：

h＝h_ph-h_{ocean_tide}-h_{lpe_tide}-h_IB

式中，h表示校正光子高度，h_ph表示卫星测得的光子高度，h_{ocean_tide}、h_{lpe_tide}表示海洋潮汐校正高度值以及长周期平衡潮汐校正高度值，h_IB表示大气校正高度值。

进一步地，步骤S6中所述冰水判断规则如下表所示：

进一步地，所述每个光子单元的海冰干舷计算后与ICESat-2官方数据作比较，误差小于1％时，选用官方数据作为后续海冰干舷数据的输入。

进一步地，步骤S8中所述固定倍数为0.3～0.7。

进一步地，步骤S10中所述线性回归函数表达式为：

式中，h_sl表示初始积雪深度值，h_fl是表示海冰干舷高度值，c₁到c₅是根据冰桥计划数据计算得到的回归系数，

表示积雪深度截断值，

表示海冰干舷高度截断值。

进一步地，步骤S10中所述积雪深度重建值计算表达式为：

式中，

表示海冰干舷高度截断值，h_f表示总干舷值，大小为海冰干舷高度值与积雪深度值的和。

进一步地，步骤S3中所述多倍可取范围为2-3倍。

进一步地，步骤S11中所述阿基米德浮力定律表达式为：

式中，ρ_w表示水的密度值，h_s和ρ_s分别表示积雪深度重建值和积雪密度值，ρ_i表示海冰的密度，h_fl表示海冰干舷高度值，h_f表示总干舷高度值，大小为h_s+h_fl，h_i表示海冰厚度值。

进一步地，所述背景光子率需要进行归一化处理，计算表达式为：

r_{bkg_norm}＝r_bkg*cos(90-θ_ref)/cos(90-min[θ_s,θ_nlb])；

式中，θ_ref表示太阳高度角，θ_nlb表示背景光子率的下界，θ_s表示该光子单元对应的太阳高度角，r_bkg表示背景光子率，r_{bkg_norm}表示归一化的背景光子率。

与现有技术相比，本发明具有以下有益效果：

1、对卫星测得光子高度进行了地球物理校正，规避了地球大气潮汐等影响，并以校正后的数据进行了粗表面提取和精表面提取，进行了参考前序步骤所得数据的冰水分类，并通过相关分类方法计算不同情况下的海冰干舷高度值，同时对数据进行了筛选，提高了有效数据的利用率。

2、与现有卫星的海冰干舷数据进行对比，并同时规定了两者数据误差较小时选用卫星数据作为后续的输入，使用卫星数据作为发明数据的参考、校正，使得获取的海冰信息更准确。

3、在计算海冰干舷时，精表面提取步骤中使用了高斯拟合的方式提取高度，与ICESat-2的数据相比，简化了复杂的计算步骤，提升了计算效率。

附图说明

图1为本发明的方法流程示意图。

图2为本发明与ICESat-2提取的海冰干舷对比示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本实施例公开了一种海冰信息提取方法，流程图如图1所示，包括：

步骤S1、使用多波束单光子激光测高卫星获取初始光子高度，对初始光子高度进行地球物理校正，地球物理校正公式为：

h＝h_ph-h_{ocean_tide}-h_{lpe_tide}-h_IB

式中，h表示校正光子高度，h_ph表示卫星测得的光子高度，h_{ocean_tide}、h_{lpe_tide}表示海洋潮汐校正高度值以及长周期平衡潮汐校正高度值，h_IB表示大气校正高度值。得到校正光子高度后，将校正光子高度减去现有平均海平面高度(MSS)，得到相对光子高度。

步骤S2、以10km以内的所有光子作为一个光子单元，可得到连续的多个光子单元，对每个光子单元进行粗表面提取：

根据每个光子的相对光子高度，对每个光子单元建立第一高度直方图，在高度直方图中，每个光子的高度值形成一个光子柱，保留[-15m,15m]高度范围内的光子，得到第二高度直方图，提取第二高度直方图中的上下波峰的20％作为第二设定高度范围，保留第二设定高度范围内的光子，得到第三高度直方图，选取出第二设定高度范围内的光子作为选中光子，计算选中光子的高度平均值：

若当前光子单元的高度平均值大于3.0m，则舍去这个光子单元的前一半长度中的光子，以当前光子单元后一半长度的光子和下一个光子单元前一半长度的光子作为新的光子单元重新计算高度平均值；

若当前光子单元的高度平均值小于3.0m，则将该高度平均值作为当前光子单元的粗表面高度。

步骤S3、对经过粗表面提取的每个光子单元数据进行精表面提取：

将当前光子单元中所有光子的相对光子高度减去粗表面提取中得到的粗表面高度，得到差值光子高度，保留在[-2.0m,3.5m]高度范围内的光子，以设定长度中的固定数量的光子作为一个迷你光子单元，此处的设定长度为10-20m，固定数量为150个。根据差值光子高度，建立第四高度直方图，提取第四高度直方图中每个迷你光子单元的波峰，并判断当前迷你光子单元的波峰是否在第四高度直方图的第一或最后一列，且光子总数大于120个：

若是，则计算差值光子高度的平均值与标准差，保留2倍标准差范围内的差值光子高度，以该数据建立第五高度直方图并进行高斯拟合，重新计算差值光子高度及展宽，将差值光子高度加上此前减去的粗表面高度，得到当前迷你光子单元的所有光子的精表面高度；若否，则将该迷你光子单元的数据删除，对下一个迷你光子单元进行上述判断。

步骤S4、根据检测到的光子数，结合光子偏差校正公式，得到光子实际到达多波束单光子激光测高卫星的探测器的时间。

步骤S5、根据光子实际到达多波束单光子激光测高卫星的探测器的时间，计算每个迷你光子单元的表面光子率和背景光子率，得到背景光子率后要对其进行归一化操作，计算表达式为：

r_{bkg_norm}＝r_bkg*cos(90-θ_ref)/cos(90-min[θ_s,θ_nlb])；

θ_ref表示太阳高度角，θ_nlb表示背景光子率的下界，θ_s表示该光子单元对应的太阳高度角，r_bkg表示背景光子率，r_{bkg_norm}表示归一化的背景光子率。

步骤S6、根据季节、表面光子率、归一化的背景光子率、展宽以及太阳高度角，结合冰水判断规则对每一个光子单元中的所有迷你光子单元进行冰水分类，冰水判断规则基本可如下所示：

步骤S7、根据冰水分类结果，对每一个光子单元计算海冰干舷：

若一个光子单元中存在一个冰间水道，其中每个冰间水道构成一个冰间水道单元；

若一个光子单元存在多个冰间水道，则其中连续的多个冰间水道可构成一个冰间水道单元，以每个冰间水道单元内冰间水道的长度为权重，计算冰间水道单元精表面高度的加权平均值并设为其高度，再对多个冰间水道单元依据冰间水道的长度进行加权平均计算整个光子单元的高度值作为该光子单元的海冰干舷高度；

若一个光子单元中不存在冰间水道，则先判断相邻的存在冰间水道的光子单元是否有效，不存在冰间水道的光子单元若处于有效的光子单元之间，且这两个有效的光子单元中光子实际到达探测器的平均时间的时间间隔小于30s以及海冰干舷高度差值绝对值小于0.05m，则可通过线性插值的方法计算中间光子单元的海冰干舷高度；

判断光子单元是否有效的方式为计算该光子单元到陆地的距离，当所述距离小于25km或者海冰密集度小于50％时，则将当前光子单元设为无效的光子单元，反之为有效的高度单元；

在得到海冰干舷高度之后，可将每个光子单元的海冰干舷计算后与ICESat-2官方数据作比较，对比示意图如图2所示，当误差小于1％时，选用官方测得的数据作为后续海冰干舷数据的输入。

步骤S8、使用Warren Snow气候学模型计算得到积雪深度，可以通过以下公式计算出来：

H＝H₀+Ax+By+Cxy+Dx²+Ey² (2)

式中，x,y为经纬度坐标，H为雪深，H₀、A、B、C、D、E为最小二乘拟合参数，以其计算值的0.5倍作为初始积雪深度值。

步骤S9、根据初始积雪深度值和所述海冰干舷高度值，使用线性回归函数计算得到积雪深度截断值和海冰干舷高度截断值，线性回归函数表达式为：

式中，h_sl表示初始积雪深度值，h_fl是表示海冰干舷高度值，c₁到c₅是根据冰桥计划数据计算得到的回归系数(c₁＝0.70,c₂＝0.22，c₃＝0.16，c₄＝1.03，c₅＝0.01m)，

表示积雪深度截断值，

表示海冰干舷高度截断值。

步骤S10、根据积雪深度截断值和海冰干舷高度值，计算出积雪深度重建值，积雪深度重建值计算表达式为：

式中，

表示海冰干舷高度截断值，h_f表示总干舷值，大小为海冰干舷高度值与积雪深度值的和。

步骤S11、根据所述积雪深度重建值、所述海冰干舷高度值、水的密度值、海冰密度值和积雪密度值，结合阿基米德浮力定律计算出海冰厚度，阿基米德浮力定律表达式为：

式中，ρ_w表示水的密度值，h_s和ρ_s分别表示积雪深度重建值和积雪密度值，ρ_i表示海冰的密度，h_fl表示海冰干舷高度值，h_f表示总干舷高度值，大小为h_s+h_fl，h_i表示海冰厚度值。其中海冰密度值为915kg/m^-3，积雪密度值由Warren Snow气候学模型计算得到。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种海冰信息提取方法，其特征在于，包括：

S1、使用多波束单光子激光测高卫星获取初始光子高度，对所述初始光子高度进行地球物理校正，得到校正光子高度，将所述校正光子高度减去现有平均海平面高度，得到相对光子高度；

S2、以第一设定长度中的所有光子作为一个光子单元，可得到连续的多个光子单元，对每个光子单元进行粗表面提取：

根据每个光子的相对光子高度，对每个光子单元建立第一高度直方图，在高度直方图中，每个光子的高度值形成一个光子柱，保留第一设定高度范围内的光子，得到第二高度直方图，提取所述第二高度直方图中的上下波峰的设定百分比作为第二设定高度范围，保留第二设定高度范围内的光子，得到第三高度直方图，选取出第二设定高度范围内的光子作为选中光子，计算选中光子的高度平均值：

若当前光子单元的高度平均值大于差值阈值，则舍去这个光子单元的前一半长度中的光子，以当前光子单元后一半长度的光子和下一个光子单元前一半长度的光子作为新的光子单元重新进行计算高度平均值；

若当前光子单元的高度平均值小于差值阈值，则将所述高度平均值作为当前光子单元的粗表面高度；

S3、对经过粗表面提取的每个光子单元数据进行精表面提取：

将当前光子单元中所有光子的相对光子高度减去粗表面提取中得到的粗表面高度，得到差值光子高度，保留第三设定高度范围内的光子，以第二设定长度中的固定数量的光子作为一个迷你光子单元，根据差值光子高度，建立第四高度直方图，提取所述第四高度直方图中每个迷你光子单元的波峰，并判断当前迷你光子单元的波峰是否在所述第四高度直方图的第一或最后一列，且光子总数大于光子阈值数：

若是，则计算差值光子高度的平均值与标准差，保留多倍标准差范围内的差值光子高度，以该数据建立第五高度直方图并进行高斯拟合，重新计算差值光子高度及展宽，将差值光子高度加上此前减去的粗表面高度，得到当前迷你光子单元的所有光子的精表面高度；

若否，则将该迷你光子单元的数据删除，对下一个迷你光子单元进行上述判断；

S4、根据检测到的光子数，结合光子偏差校正公式，得到光子实际到达多波束单光子激光测高卫星的探测器的时间；

S5、根据光子实际到达多波束单光子激光测高卫星的探测器的时间，计算每个迷你光子单元的表面光子率和背景光子率；

S6、根据季节、表面光子率、归一化的背景光子率、展宽以及太阳高度角，结合冰水判断规则对每一个光子单元中的所有迷你光子单元进行冰水分类，所述冰水分类类型中包含冰间水道；

S7、根据冰水分类结果，对每一个光子单元计算海冰干舷：

若一个光子单元中存在一个冰间水道，其中每个冰间水道构成一个冰间水道单元；

若一个光子单元存在多个冰间水道，则其中连续的多个冰间水道可构成一个冰间水道单元，以每个冰间水道单元内冰间水道的长度为权重，计算所述冰间水道单元精表面高度的加权平均值并设为其高度，再对多个冰间水道单元依据冰间水道的长度进行加权平均计算整个光子单元的高度值作为该光子单元的海冰干舷高度；

若一个光子单元中不存在冰间水道，则先判断相邻的存在冰间水道的光子单元是否有效，不存在冰间水道的光子单元若处于有效的光子单元之间，且这两个有效的光子单元中光子实际到达探测器的平均时间的时间间隔小于时间阈值以及海冰干舷高度差值绝对值小于高度差阈值，则可通过线性插值的方法计算中间光子单元的海冰干舷高度；

判断光子单元是否有效的方式为计算该光子单元到陆地的距离，当所述距离小于距离阈值或者海冰密集度小于密集度阈值时，则将当前光子单元设为无效的光子单元，反之为有效的高度单元；

S8、使用沃伦·斯诺气候学模型计算得到积雪深度，并以其计算值的固定倍数作为初始积雪深度值；

S9、根据所述初始积雪深度值和所述海冰干舷高度值，使用线性回归函数计算得到积雪深度截断值和海冰干舷高度截断值；

S10、根据所述积雪深度截断值和海冰干舷高度值，计算出积雪深度重建值；

S11、根据所述积雪深度重建值、所述海冰干舷高度值、水的密度值、海冰密度值和积雪密度值，结合阿基米德浮力定律计算出海冰厚度。

2.根据权利要求1所述的一种海冰信息提取方法，其特征在于，所述地球物理校正公式为：

h＝h_ph-h_{ocean_tide}-h_{lpe_tide}-h_IB

式中，h表示校正光子高度，h_ph表示卫星测得的光子高度，h_ocean-tide、h_lpe-tide表示海洋潮汐校正高度值以及长周期平衡潮汐校正高度值，h_IB表示大气校正高度值。

3.根据权利要求1所述的一种海冰信息提取方法，其特征在于，步骤S6中所述冰水判断规则如下表所示：

4.根据权利要求1所述的一种海冰信息提取方法，其特征在于，所述每个光子单元的海冰干舷计算后与新一代冰、云和陆地高程卫星官方数据作比较，误差小于1％时，选用官方数据作为后续海冰干舷数据的输入。

5.根据权利要求1所述的一种海冰信息提取方法，其特征在于，步骤S8中所述固定倍数为0.3～0.7。

6.根据权利要求1所述的一种海冰信息提取方法，其特征在于，步骤S10中所述线性回归函数表达式为：

式中，h_sl表示初始积雪深度值，h_fl是表示海冰干舷高度值，c₁到c₅是根据冰桥计划数据计算得到的回归系数，

表示积雪深度截断值，

表示海冰干舷高度截断值。

7.根据权利要求1所述的一种海冰信息提取方法，其特征在于，步骤S10中所述积雪深度重建值计算表达式为：

式中，

表示海冰干舷高度截断值，

表示积雪深度截断值，h_f表示总干舷值，大小为海冰干舷高度值与积雪深度值的和。

8.根据权利要求1所述的一种海冰信息提取方法，其特征在于，步骤S3中所述多倍可取范围为2～3倍。

9.根据权利要求1所述的一种海冰信息提取方法，其特征在于，步骤S11中所述阿基米德浮力定律表达式为：

式中，ρ_w表示水的密度值，h_s和ρ_s分别表示积雪深度重建值和积雪密度值，ρ_i表示海冰的密度，h_fl表示海冰干舷高度值，h_f表示总干舷高度值，大小为h_s+h_fl，h_i表示海冰厚度值。

10.根据权利要求1所述的一种海冰信息提取方法，其特征在于，所述背景光子率需要进行归一化处理，计算表达式为：

r_{bkg_norm}＝r_bkg*cos(90-θ_ref)/cos(90-min[θ_s，θ_nlb])；

式中，θ_ref表示太阳高度角，θ_nlb表示背景光子率的下界，θ_s表示该光子单元对应的太阳高度角，r_bkg表示背景光子率，r_{bkg_norm}表示归一化的背景光子率。

Images