CN106197383A - 一种海冰体积的遥感估算方法 - Google Patents

Abstract

本发明公开了一种海冰体积的遥感估算方法，步骤包括从CryoSat‑2数据中提取经纬度及对应浮冰出水高度、积雪深度和海冰密集度信息；下载海冰类型数据，并提取经纬度和海冰类型信息；将上述数据进行匹配并转换为矢量数据；通过投影变换、栅格转换和空间重采样、掩膜提取处理，获得浮冰出水高度、积雪深度、海冰密集度和海冰类型栅格数据；对一年冰和多年冰赋予不同的海冰密度参数，通过流体静力学模型估算每个像元的海冰厚度；像元的实际面积大小乘以该像元对应的海冰密集度以及海冰厚度，最终得到整个海域的海冰体积。本方法考虑了积雪深度、积雪密度、一年冰和多年冰密度的差别对海冰厚度反演结果的影响，进而提高了海冰体积的估算精度。

Description

一种海冰体积的遥感估算方法

技术领域

本发明涉及一种海冰体积的遥感估算方法，属于遥感地学应用技术领域。

技术背景

两极海冰作为地球系统的冷源，对全球大气和海洋环流有重要和长期的影响，是全球气候变化的重要指示器。近年来，全球气候变化的影响越来越明显，人们越来越需要准确了解两极海冰正在发生的变化。海冰厚度和体积的变化日益受到科学界的关注。

随着空间技术的发展，利用卫星和航空航天对地观测技术获取海冰遥感资料，再通过模型算法得出海冰厚度与形态特征已成为重要的途径。2010年4月，欧空局发射的专门为冰冻圈设计的CryoSat-2卫星，克服了传统测高卫星ERS-1、ERS-2、ENVISAT轨道倾角的限制，观测数据范围最高达到北极88°N的区域。另外，其独特的多视处理和高分辨率弥补了传统遥感观测手段和测高卫星的不足。该卫星上携带的合成孔径干涉雷达高度计(SIRAL)的垂直测量精度达到1～3cm，相比于先前欧空局雷达高度计(RA)约10km的卫星足迹，它采用延迟多普勒雷达高度计(DDA)技术将卫星地面足迹减小到沿轨约为0.3km，跨轨约为1.5km。另外，SIRAL对地表点进行多视处理以减少雷达斑点引起的噪声。其测量海平面高度数据精度约为传统雷达高度计的2倍。

欧洲气象卫星应用组织(EUMETSAT)海洋海冰应用中心(OSI SAF)的核心任务是监测南北极冰雪覆盖的变化。自2005年以来，提供全球海冰类型、密集度、外缘线和反射率等数据产品。其海冰类型产品是使用贝叶斯算法和大气校正模型，基于被动微波散射计(SSMIS)的亮度温度和主动微波散射计(ASCAT)的后向散射数据，反演得到的空间分辨率为10km的极地立体投影网格数据。两者的有效结合，对反演海冰厚度和估算海冰体积有重要意义。

发明内容

本发明要解决的技术问题是：克服现有技术的上述缺陷，提供一种基于CryoSat-2卫星测高数据和OSI SAF海冰类型数据反演海冰厚度，进而估算海冰体积的方法。本发明海冰体积估算方法考虑了积雪深度、多年冰和一年冰海冰密度的差别，大大提高了海冰厚度和体积的估算精度。

为了解决上述技术问题，本发明提出的技术方案是：一种海冰体积的遥感估算方法，包括以下步骤：

第一步、准备数据，并对数据进行提取和筛选，包括如下两个方面的内容：

Ⅰ.下载CryoSat-2卫星SAR模式L2I级数据产品，并提取经纬度坐标、浮冰出水高度、雪深和海冰密集度信息，并根据指示性参数剔除不可用的数据；

II.下载对应时间和空间范围的欧洲气象卫星应用组织海洋海冰应用中心的海冰类型数据，获取海冰类型信息以及对应的经纬度信息，所述海冰类型分为一年冰和多年冰；

第二步、分别将上述提取的浮冰出水高度、雪深、海冰密集度以及海冰类型信息与对应的经纬度坐标进行匹配，并将匹配好的数据转换成矢量点数据；

第三步、获取浮冰出水高度栅格数据、雪深栅格数据、海冰密集度栅格数据和海冰类型栅格数据，并将上述栅格数据转换为25km×25km的空间分辨率；

第四步、以浮冰出水高度栅格数据为掩膜对海冰类型栅格数据进行处理，获得CryoSat-2数据覆盖范围内浮冰的海冰类型；

第五步、通过流体静力平衡公式计算每个像元的海冰厚度值H_i，并生成海冰厚度的栅格数据，公式表达如下，

$H_i=\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{F}_i+\frac{{\mathrm{\ρ}}_s}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{H}_s$

其中，ρ_w、ρ_i和ρ_s分别是海水密度、海冰密度和海冰上覆盖积雪的密度，F_i为浮冰出水高度，H_s为海冰上的积雪深度，一年冰和多年冰具有各自的海冰密度值；

第六步、每个像元所代表的实际面积乘以该像元对应的海冰密集度及海冰厚度值，得到该像元对应的海冰体积，所有像元的海冰体积累加，得到整个海域的海冰体积，公式表达如下：

$V=\underset{j=1}{\overset{n}{\Σ}}(S\×C_j\×H_j)$

其中V表示整个海域的海冰体积，C_j、H_j分别代表第j个像元的海冰密集度和海冰厚度，S为像元的实际面积，即像元的空间分辨率，n为像元总个数。

本方法数据准备阶段利用MATLAB软件对数据进行了提取、筛选、匹配和矢量转换，以便利用ArcGIS进行后续的处理，使得处理过程更加高效准确。

本发明海冰体积的遥感估算方法，还具有如下特征：

1、所述第一步中，CryoSat-2卫星SAR模式L2I级数据产品为DBL格式，读取DBL格式的原始文件获得数据集，从数据集中提取经纬度、浮冰出水高度、雪深和海冰密集度信息；欧洲气象卫星应用组织海洋海冰应用中心海冰类型数据格式为NetCDF4格式。提取数据后，剔除异常值，减小了数据量，提高了数据处理效率。

2、所述第二步中，使用MATLAB软件中shapewrite函数将匹配好的数据转换成带经纬度坐标的矢量点数据，存储为shapefile文件。将提取的数据转换成带经纬度坐标的矢量点，有利于ArcGIS的后续处理。

3、所述第三步具体包括如下步骤：

a.对第二步中获得的矢量点数据在WGS_1984地理坐标系下进行极地立体方位投影；

b.将极地立体方位投影后的矢量点数据进行栅格转换，获得浮冰出水高度栅格数据、积雪深度栅格数据、海冰密集度栅格数据和海冰类型栅格数据；

c.通过空间重采样将上述栅格数据的空间分辨率采样至相同分辨率。

获取同一坐标系下统一空间分辨率的浮冰出水高度、雪深、海冰类型和海冰密集度栅格数据，有利于海冰类型和浮冰信息相匹配，从而确定浮冰类型。

4、所述第四步中，海水密度ρ_w取1023.8kg/m³，一年海冰和多年海冰的密度分别取916.7kg/m³和882.0kg/m³，积雪密度取400.0kg/m³。

海冰密度误差在海冰厚度反演中占据主要地位，一年冰密度大于多年冰。对一年冰和多年冰赋予不同的密度参数，提高了厚度反演的精度。海水和海冰的密度参考了其他文献，例如Wadhams等人2010年发表于The Cryosphere的“The relation between sea icethickness and freeboard in the Arctic”通过689个观测站点的数据得到一年冰的平均密度是916.7±35.7kg/m³，多年冰的平均密度是882±23kg/m³，当测量值较多时，将一年冰和多年冰的密度设定为916.7kg/m³和882.0kg/m³是可行的。Rachel等人2015年发表于Nature的“Increased Arctic sea ice volume after anomalously low melting in2013”将海水密度设定为1023.8kg/m³，一年冰和多年冰的密度分别为916.7kg/m³和882.0kg/m³。

5、所述第五步中，通过ArcGIS软件中的Raster Calculator工具，根据流体静力学平衡模型反演得到海冰厚度。

本发明利用了CryoSat-2卫星数据和OSI SAF海冰类型数据，数据获取简单。CryoSat-2SIRAL测高数据能提供精确的浮冰出水高度信息，另外其辅助测量数据中还包含了积雪深度、海冰密集度的信息。OSI SAF海冰类型数据是通过贝叶斯算法和大气校正模型反演得到的空间分辨率为10km的极地立体投影网格数据。两者的有效结合提高了海冰厚度的反演精度，为海冰体积的估算提供了数据基础。

本发明数据提取和矢量转换过程均通过MATLAB编程自动实现，大大减少了工作量，同时减少了人工误差。后续栅格计算通过ArcGIS中的Raster Calculator工具完成。

综上，本发明方法的执行步骤简单易行，效果较好。目前海冰厚度大面积观测资料较少，传统的经验方法估算海冰厚度精度不高。本发明使用高精度大范围的SIRAL浮冰出水高度数据和OSI SAF海冰类型产品相结合，根据流体静力学平衡原理反演海冰厚度，并估算海冰体积。对精确量化南北极海冰物质平衡变化有深远的科学意义。

附图说明

下面结合附图对本发明作进一步的说明。

图1是本发明海冰体积遥感估算方法流程图。

图2(a)为提取的浮冰出水高度、雪深、海冰密集度信息。

图2(b)为海冰类型信息。

图3(a)、(b)、(c)、(d)分别为25km空间分辨率的浮冰出水高度栅格数据、雪深栅格数据、海冰密集度栅格数据、海冰类型栅格数据的图像。

图4为海冰厚度结果显示图。

图5为一年冰和多年冰海冰体积的估算结果。

具体实施方式

下面根据附图详细阐述本发明，使本发明的技术路线和操作步骤更加清晰。本发明实例采用的数据是CryoSat-2卫星L2I级SAR模式基线C数据和OSI SAF10km分辨率网格数据。CryoSat-2数据获取时间为2016年1月1日——2016年1月31日，OSI SAF数据获取时间为2016年1月15日。

图1为海冰体积遥感估算方法流程图，具体步骤如下:

第一步、准备数据阶段，具体操作步骤如下：

Ⅰ.依次读取CryoSat-2数据3597个DBL格式原始文件，获得数据集。从数据集中提取每个经纬度坐标点的浮冰出水高度、雪深和海冰密集度信息，并提出不可用数据；

II.下载对应时间的欧洲气象卫星应用组织(EUMETSAT)海洋海冰应用中心(OSISAF)海冰类型数据。读取OSI SAF数据NetCDF4格式文件，获取海冰类型以及对应的经纬度信息。

第二步、分别将上述提取的浮冰、雪深、密集度以及海冰类型信息与对应的经纬度匹配。使用shapewrite函数将上述匹配好的数据转换成带经纬度坐标的矢量点数据。分别为矢量文件a和b。其中矢量文件a包含了浮冰出水高度、积雪深度、海冰密集度信息。矢量文件b包含了海冰类型信息。矢量数据的属性值如图2(a)、图2(b)所示。其中图2(b)中属性值Ice_type为2代表一年冰，为3代表多年冰。

第三步、使用ArcGIS软件获取浮冰出水高度、雪深、海冰类型以及海冰密集度栅格数据，并将上述数据统一为25km空间分辨率的栅格数据集，结果如图3所示。用海冰类型鉴别浮冰。

a.定义投影。将矢量文件a和b加载到ArcGIS中，定义WGS_1984地理坐标系，进行极地立体方位投影。

b.获得栅格数据集。利用ArcGIS工具Conversion Tools->To Raster->Point toRaster，计算方式选择mean，栅格大小为默认值。获得浮冰出水高度、积雪深度、海冰密集度和海冰类型栅格数据。

c.统一空间分辨率。利用ArcGIS工具Data Management Tools->Raster->RasterProcessing->Resample，选取最邻近分配法，将上述栅格数据统一为25km×25km空间分辨率。其中浮冰出水高度、雪深和海冰密集度栅格数据分别如图3(a)、3(b)、3(c)。

第四步、掩膜提取。以浮冰出水高度数据掩膜提取海冰类型，从而获得CryoSat-2数据覆盖范围浮冰的海冰类型，如图3(d)所示。本例中，利用ArcGIS工具Spatial AnalystTools->Extraction->Extract by Mask。

第五步、通过流体静力平衡公式计算每个像元的海冰厚度值H_i，并生成栅格图层，如图4所示。公式表达如下，

$H_i=\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{F}_i+\frac{{\mathrm{\ρ}}_s}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{H}_s$

式中，ρ_w、ρ_i和ρ_s分别是海水密度、海冰密度和海冰上覆盖积雪的密度，F_i为浮冰出水高度，H_s为海冰上的积雪深度。本例中，海水密度ρ_w取1023.8kg/m³，一年海冰和多年海冰的密度分别取916.7kg/m³和882.0kg/m³，积雪密度取400.0kg/m³。

第六步、估算海冰体积。每个像元所代表的实际面积乘以该像元对应的海冰密集度及海冰厚度值，得到该像元对应的海冰体积，所有像元对应的海冰体积相加，得到整个海域的海冰体积，估算结果如图5所示，其中一年冰体积约为16.33×10³km³，多年冰体积约为7.22×10³km³。具体公式表达如下：

$V=\underset{j=1}{\overset{n}{\Σ}}(S\×C_j\×H_j)$

其中V表示整个海域的海冰体积，C_j、H_j分别代表第j个像元的海冰密集度和海冰厚度，S为像元的实际面积，即像元的空间分辨率，n为像元总个数。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。

Claims (6)

1.一种海冰体积的遥感估算方法，包括以下步骤：

第一步、准备数据，并对数据进行提取和筛选，包括如下两个方面的内容：

Ⅰ.下载CryoSat-2卫星SAR模式L2I级数据产品，并提取经纬度坐标、浮冰出水高度、雪深和海冰密集度信息，并根据指示性参数剔除不可用的数据；

II.下载对应时间和空间范围的欧洲气象卫星应用组织海洋海冰应用中心的海冰类型数据，获取海冰类型信息以及对应的经纬度信息，所述海冰类型分为一年冰和多年冰；

第二步、分别将上述提取的浮冰出水高度、雪深、海冰密集度以及海冰类型信息与对应的经纬度坐标进行匹配，并将匹配好的数据转换成矢量点数据；

第三步、获取浮冰出水高度栅格数据、雪深栅格数据、海冰密集度栅格数据和海冰类型栅格数据，并将上述栅格数据转换为25km×25km的空间分辨率；

第四步、以浮冰出水高度栅格数据为掩膜对海冰类型栅格数据进行处理，获得CryoSat-2数据覆盖范围内浮冰的海冰类型；

第五步、通过流体静力平衡公式计算每个像元的海冰厚度值H_i，并生成海冰厚度的栅格数据，公式表达如下，

$H_i=\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{F}_i+\frac{{\mathrm{\ρ}}_s}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{H}_s$

其中，ρ_w、ρ_i和ρ_s分别是海水密度、海冰密度和海冰上覆盖积雪的密度，F_i为浮冰出水高度，H_s为海冰上的积雪深度，一年冰和多年冰具有各自的海冰密度值；

第六步、每个像元所代表的实际面积乘以该像元对应的海冰密集度及海冰厚度值，得到该像元对应的海冰体积，所有像元的海冰体积累加，得到整个海域的海冰体积，公式表达如下：

$V=\underset{j=1}{\overset{n}{\Σ}}(S\×C_j\×H_j)$

其中V表示整个海域的海冰体积，C_j、H_j分别代表第j个像元的海冰密集度和海冰厚度，S为像元的实际面积，即像元的空间分辨率，n为像元总个数。

本方法数据准备阶段利用MATLAB软件对数据进行了提取、筛选、匹配和矢量转换，以便利用ArcGIS进行后续的处理，使得处理过程更加高效准确。

2.根据权利要求1所述海冰体积的遥感估算方法，其特征在于：所述第一步中，CryoSat-2卫星SAR模式L2I级数据产品为DBL格式，读取DBL格式的原始文件获得数据集，从数据集中提取经纬度、浮冰出水高度、雪深和海冰密集度信息；欧洲气象卫星应用组织海洋海冰应用中心海冰类型数据格式为NetCDF4格式。

3.根据权利要求1所述海冰体积的遥感估算方法，其特征在于：所述第二步中，使用MATLAB软件中shapewrite函数将匹配好的数据转换成带经纬度坐标的矢量点数据，存储为shapefile文件。

4.根据权利要求1所述海冰体积的遥感估算方法，其特征在于：所述第三步具体包括如下步骤：

a.对第二步中获得的矢量点数据在WGS_1984地理坐标系下进行极地立体方位投影；

b.将极地立体方位投影后的矢量点数据进行栅格转换，获得浮冰出水高度栅格数据、积雪深度栅格数据、海冰密集度栅格数据和海冰类型栅格数据；

c.通过空间重采样将上述栅格数据的空间分辨率采样至相同分辨率。

5.根据权利要求1所述海冰体积的遥感估算方法，其特征在于：所述第四步中，海水密度取1023.8kg/m³，一年海冰和多年海冰的密度分别取916.7kg/m³和882.0kg/m³，积雪密度取400.0kg/m³。

6.根据权利要求1所述海冰体积的遥感估算方法，其特征在于：所述第五步中，通过ArcGIS软件中的Raster Calculator工具，根据流体静力学平衡模型反演得到海冰厚度。