本发明涉及一种冰下湖体积变化的估算方法，通过获取相应时间段内的CryoSat‑2雷达高度计SARIn模式的L1B波形数据；使用波形重跟踪方法计算冰下湖区域的表面高程；对获取的表面高程数据进行预处理，剔除高程变化异常点；使用DEM对预处理后的高程数据进行坡度改正，去除坡度造成的误差；使用克里金插值最终生成100*100m的格网，然后计算相应时间段内的高程变化，对高程变化值加入后向散射能量改正，去除后向散射能量对高程变化造成的影响；通过高程变化值获取冰下湖的形状和面积，对每个格网冰下湖表面的高程变化值进行积分，最终获取相应时间段内冰下湖的体积变化。本发明能够准确、有效地反映冰下湖表面的起伏，进而完成冰下湖体积变化的估算。
第一步、获取研究时间范围内研究区域的CryoSat-2雷达高度计SARIn模式的L1B波形数据；
第二步、针对地面测量点，对L1B波形数据使用波形重跟踪方法计算获得其的冰下湖表面高程值；
第三步、对每个地面测量点的冰下湖表面高程值进行预处理，包括潮汐改正、大气改正以及异常数据的剔除；
第四步、使用DEM栅格数据对预处理后每个地面测量点的冰下湖表面高程值进行坡度校正；
第五步、对经过坡度校正后的冰下湖表面高程值进行克里金插值，然后对地面测量点的冰下湖表面高程数据进行栅格化，栅格内地面测量点的冰下湖表面高程值的平均值作为栅格的冰下湖表面高程值；
第六步、对研究时间范围内各年度的年内栅格数据做平均，得到各年度的冰下湖表面高程栅格图，研究时间范围内的第一年作为参考年份，针对每个栅格，计算各年与参考年份年之间的冰下湖表面高程变化量，并对所述冰下湖表面高程变化量进行后向散射能量的校正；
第七步、当存在任意两年与参考年份年之间的冰下湖表面高程变化量的差的绝对值大于3倍卫星间高程变化标准差,则该栅格作为冰下湖所在区域进行提取；
第八步、针对冰下湖所在区域的每个栅格，以年份为横坐标，对应年份的冰下湖表面高程为纵坐标，进行迭代最小二乘拟合，得到的直线斜率为冰下湖的高程变化趋势；
第九步、每个栅格所拟合得到的直线斜率与栅格面积相乘得到体积变化率，对每个栅格的体积变化率做积分，得到冰下湖的体积变化量。
2.根据权利要求1所述冰下湖体积变化的估算方法，其特征在于：第二步中，首先使用波形重跟踪方法获取L1B波形前缘中点的序号nretTL＝DC+a(A-DC)
式中，y(n)为L1B波形的第n个采样点的采样功率，DC为前5个采样点采样功率的平均值，a＝0.4，A为L1B波形的振幅，TL为阈值水平， 为L1B波形中第一个采样功率超过阈值水平TL采样点的序号；
c为光速，t为卫星发射和接收脉冲的时间间隔，H为卫星距离椭球体表面高度，n0为L1B波形的中心采样点序号，R为相邻采样点的采样间隔与光速的积。
3.根据权利要求2所述冰下湖体积变化的估算方法，其特征在于：第二步中，所述L1B波形的振幅A通过波形重心偏移法得到。
4.根据权利要求1所述冰下湖体积变化的估算方法，其特征在于：第四步中，计算DEM栅格数据的每个栅格的坡度α，提取每个地面测量点所对应栅格的坡度，然后使用下式进行坡度校正：hslope为坡度校正后的高程值，h为第二步计算得到的地面测量点的冰下湖表面高程值h，H为卫星到地表的距离。
5.根据权利要求4所述冰下湖体积变化的估算方法，其特征在于：第四步中，首先对研究区域的1km分辨率的DEM栅格数据插值成20m分辨率，然后再计算每个栅格的坡度。
6.根据权利要求1所述冰下湖体积变化的估算方法，其特征在于：所述第五步中，栅格化完成后，在每一个栅格内用三倍标准差的阈值进行异常数据的剔除。
7.根据权利要求1所述冰下湖体积变化的估算方法，其特征在于：第六步中，通过下述公式对每个栅格的各年与参考年份年之间的冰下湖表面高程变化量进行后向散射能量的校正：ΔHcorrected＝ΔH-Δσ·CG
ΔHcorrected为后向散射能量的校正量，ΔH为相对于参考年份的冰下湖表面高程变化量，Δσ代表相对于参考年份的后向散射能量变化量，CG为ΔH和Δσ的相关梯度系数。
8.根据权利要求7所述冰下湖体积变化的估算方法，其特征在于：对与栅格对应的所有ΔH和Δσ进行直线拟合，所得到的斜率为该栅格的相关梯度系数CG。
9.根据权利要求8所述冰下湖体积变化的估算方法，其特征在于：所述第八步中，直线
拟合后，如果R＜0.95，则剔除异常数据后再次进行直线拟合，重复该步骤，直到R≥0.95，R为时间和高程的相关系数。
一种冰下湖体积变化的估算方法
技术领域
[0001] 本发明涉及一种冰下湖体积变化的估算方法，属于环境遥感应用技术领域。
背景技术
[0002] 南极冰下湖具有储存和定期释放的能力，据估计，每年通过冰下融化产生的65GT水中，有些是由这些湖泊储存和定期释放的。确定湖泊排水和再灌水的分布、周期性和体积通量，有助于冰下水文模型和冰盖流模型建立、冰下融化速率和地热通量的估计，以及对淡水向冰下陆架输送的理解。
[0003] 早期的卫星记录由于在空间或时间上取样稀少而无法捕捉整个冰下湖周期，因此其监测冰下湖体积波动的能力有限。虽然重复测高提供了定期观测，但由于离极点越来越远的地面航迹和雷达高度计的粗分辨率，冰下湖的监测受到了限制。由于卫星轨道偏离极点，限制正在变化的区域是困难的，而且体积变化估计也是不确定的。相反，干涉合成孔径雷达可以提供冰面的详细图像，但由于捕获周期短，因此冰下湖的监测受到限制。
[0004] CryoSat-2卫星的发射提供了一种新的平台，用于监测变化活跃的冰下湖。CryoSat-2卫星相比于ICESat卫星重访周期更短，重复周期为369d，并伴随着30d的子循环(369d的周期由连续变化的30d重复模式构成)；采样密度更大，卫星地面点的足迹沿轨
300m，跨轨1.5km。冰下湖的活动会引起冰盖表面高程发生变化，假定冰下湖冰盖表面的高程变化和冰下湖体积变化是等价的，通过CryoSat-2卫星SARIn模式的高程数据建立冰下湖表面的高程模型，因为冰下湖的表面高程变化比较剧烈，所以将高程变化大于3倍卫星间高程变化标准差(0.9m)的区域作为冰下湖的区域，并根据高程变化得到其体积变化。这种方法提供了更密集、更精确的高程数据，更加准确地反映了冰下湖的体积变化。
发明内容
[0005] 本发明要解决的问题是：针对之前测高卫星在空间或时间上取样稀少，无法准确获取冰下湖的形状和体积变化的困难，提供了一种冰下湖体积变化的估算方法，能够准确、有效地反映冰下湖表面的起伏，进而完成冰下湖体积变化的估算。
[0006] 为了解决上述技术问题，本发明提出的技术方案是：一种冰下湖体积变化的估算方法，包括以下步骤：
[0007] 第一步、获取研究时间范围内研究区域的CryoSat-2雷达高度计SARIn模式的L1B波形数据；
[0008] 第二步、针对地面测量点，对L1B波形数据使用波形重跟踪方法计算获得其的冰下湖表面高程值；
[0009] 第三步、对每个地面测量点的冰下湖表面高程值进行预处理，包括潮汐改正、大气改正以及异常数据的剔除；
[0010] 第四步、使用DEM栅格数据对预处理后每个地面测量点的冰下湖表面高程值进行坡度校正；
[0011] 第五步、对经过坡度校正后的冰下湖表面高程值进行克里金插值，然后对地面测量点的冰下湖表面高程数据进行栅格化，栅格内地面测量点的冰下湖表面高程值的平均值作为栅格的冰下湖表面高程值；
[0012] 第六步、对研究时间范围内各年度的年内栅格数据做平均，得到各年度的冰下湖表面高程栅格图，研究时间范围内的第一年作为参考年份，针对每个栅格，计算各年与参考年份年之间的冰下湖表面高程变化量，并对所述冰下湖表面高程变化量进行后向散射能量的校正；
[0013] 第七步、当存在任意两年与参考年份年之间的冰下湖表面高程变化量的差的绝对值大于3倍卫星间高程变化标准差,则该栅格作为冰下湖所在区域进行提取；
[0014] 第八步、针对冰下湖所在区域的每个栅格，以年份为横坐标，对应年份的冰下湖表面高程为纵坐标，进行迭代最小二乘拟合，得到的直线斜率为冰下湖的高程变化趋势；
[0015] 第九步、每个栅格所拟合得到的直线斜率与栅格面积相乘得到体积变化率，对每个栅格的体积变化率做积分，得到冰下湖的体积变化量。
[0016] 本发明的有益效果是：
[0017] 获取南极冰下湖体积变化对确定冰下湖排水和再灌水的分布、周期性和体积通量有助于冰下水文模型和冰盖流模型的建立，冰下融化速率和地热通量的估计以及对淡水向冰下陆架输送的理解有重要的意义。本发明实现了南极冰下湖体积变化的估算方法，通过CryoSat-2雷达高度计SARIn模式的L1B波形数据，计算表面高程，对表面高程数据以及高程变化数据进行多种校正，准确获取了冰下湖的形状和面积以及表面的起伏，进而得到了冰下湖的体积变化，相对于其他方法有无可比拟的优势。具体的有益效果如下：
[0018] 第一、本发明成功获取了冰下湖的面积和形状，以及其表面高程变化和体积变化，可进一步分析确定湖泊排水和再灌水的分布、周期性和体积通量，有助于冰下水文模型和冰盖流模型的建立，冰下融化速率和地热通量的估计以及对淡水向冰下陆架输送的理解。
[0019] 第二、本发明使用的数据获取便捷，可以直接在ESA网站上下载CryoSat-2雷达高度计SARIn模式的L1B波形数据，即可进行冰下湖体积变化的研究工作。
[0020] 第三、本发明采用CryoSat-2雷达高度计SARIn模式的L1B波形数据，使用波形重跟踪的方法计算表面高程，并对表面高程进行坡度改正以及后向散射能量改正，获取冰下湖的面积和形状，使用迭代最小二乘法进行高程变化趋势的拟合，进而获取冰下湖体积变化。
[0021] 第四、本发明方法的执行步骤简单易行，效果较好，解决了冰下湖形状和面积难以确定，体积变化无法准确获取的难题。促进了冰下湖水文过程的研究。
附图说明
[0022] 下面结合附图对本发明作进一步说明。
[0023] 图1是本发明冰下湖体积变化的估算方法的流程图。
[0024] 图2是2011-2017年冰下湖表面高程。
[0025] 图3是使用三倍卫星间高程变化的标准差为阈值提取的冰下湖形状。
[0026] 图4是2011-2017年冰下湖表面的高程变化趋势。
[0027] 图5是图3中A点的高程变化。
[0028] 图6是2011-2017年CookE2冰下湖总体的体积变化趋势。
具体实施方式
[0029] 下面根据附图详细阐述本发明，使本发明的目的和效果变得更加明显。
[0030] 本发明估算南极冰下湖体积变化所采用的CryoSat-2雷达高度计SARIn模式L1B波形数据是从欧空局(ESA)网站(ftp://science–pds.cryosat.esa.int)下载，时间为2011年1月-2017年12月。
[0031] 如图1为本实施例的流程图，本实施例冰下湖体积变化的估算方法包括以下步骤：
[0032] 第一步、获取CookE2冰下湖区域的2011-2017年CryoSat-2雷达高度计SARIn模式的L1B波形数据。研究区的范围为155.1E-156.3E，72.69S-76.97S。
[0033] 第二步、针对地面测量点，对CryoSat-2 SARIn模式的L1B波形数据使用波形重跟踪方法(阈值为0.4)计算冰下湖的表面高程，具体如下：
[0034] 首先使用波形重跟踪算法获取L1B波形前缘中点的序号nret
[0035]
[0036] TL＝DC+a(A-DC)
[0037]
[0038] 式中，y(n)为L1B波形的第n个采样点的采样功率，DC为前5个采样点采样功率的平均值，a为阈值,a＝0.4，A为L1B波形的振幅，TL为阈值水平，为L1B波形中第一个采样功率超过阈值水平TL采样点的序号；
[0039] 然后计算地面测量点的冰下湖表面高程值h：
[0040] h＝c·t-H+(n0-nret)·R
[0041] c为光速，t为卫星发射和接收脉冲的时间间隔，H为卫星距离椭球体表面高度，n0为L1B波形的中心采样点序号，R为相邻采样点的采样间隔与光速的积。
[0042] 第三步、对每个地面测量点的冰下湖表面高程值进行预处理，包括潮汐改正、大气改正、干湿对流层改正、逆气压改正、电离层改正等，以及异常数据的剔除。
[0043] 第四步、使用1km分辨率的DEM栅格数据通过克里金插值生成分辨率为20m的DEM栅格数据,使用ArcGIS将插值后的DEM生成坡度数据，提取每个地面测量点所对应栅格的坡度α，然后使用下式进行坡度校正：
[0044]
[0045] hslope为坡度校正后的高程值，h为第二步计算得到的地面测量点的冰下湖表面高程值h，H为卫星到地表的距离。
[0046] 第五步、对经过坡度校正后的冰下湖表面高程值进行克里金插值，然后对地面测量点的冰下湖表面高程数据进行栅格化，栅格化的空间分辨率为100×100m，在每一个栅格内用三倍标准差的阈值进行异常数据的剔除，将剔除后的冰下湖表面高程值取平均作为栅格的冰下湖表面高程值。
[0047] 第六步、对研究时间范围内各年度的年内栅格数据做平均，得到2011-2017各年度的冰下湖表面高程栅格图，见图2。研究时间范围内的第一年作为参考年份，针对每个栅格，计算各年与参考年份年之间的冰下湖表面高程变化量，并对所述冰下湖表面高程变化量进行后向散射能量的校正。
[0048] 本步骤中，通过下述公式对每个栅格的各年与参考年份年之间的冰下湖表面高程变化量进行后向散射能量的校正：
[0049] ΔHcorrected＝ΔH-Δσ·C3
[0050] ΔHcorrected为后向散射能量的校正量，ΔH为相对于参考年份的冰下湖表面高程变化量，Δσ代表相对于参考年份的后向散射能量变化量，CG为ΔH和Δσ的相关梯度系数，对与栅格对应的所有ΔH和Δσ进行直线拟合，所得到的斜率即为相关梯度系数CG。
[0051] 第七步、当存在任意两年与参考年份年之间的冰下湖表面高程变化量的差(0.9m)的绝对值大于3倍卫星间高程变化标准差,则该栅格作为冰下湖所在区域进行提取，进而确定冰下湖的形状和面积。冰下湖提取结果见图3。
[0052] 第八步、计算2011-2017年冰下湖的高程变化趋势。针对冰下湖所在区域的每个栅格，以年份为横坐标，对应年份的冰下湖表面高程为纵坐标，进行直线拟合，直线拟合得到的直线斜率为冰下湖的高程变化趋势。
[0053] 本步骤中，直线拟合后，如果R2<0.95，则剔除异常数据后再次进行直线拟合，重复该步骤，直到R2≥0.95，R为时间和高程的相关系数。
[0054] 高程变化趋势图见图4，图中颜色越亮和越深，代表冰下湖高程年变化越大。图5为图3中A点的高程年变化图。
[0055] 第九步、每个栅格所拟合得到的直线斜率与栅格面积相乘得到体积变化率，对每个栅格的体积变化率做积分，得到整个冰下湖2011-2017年的体积变化量。如图6所示，为2011-2017年CookE2冰下湖总体的体积变化趋势。对比图5和图6，可知：2011年到2017年CookE2冰下湖存在周期性变化(先升高后降低再升高)，但总体处于微小的上升趋势。
[0056] 除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。
