基于遥感技术的1988—2018年渤海海冰冰情时空变化特征*

2022-03-14孙静琪陈笑娟高吉喆许映军

灾害学 2022年1期

孙静琪，陈笑娟，李倩，高吉喆，张化，许映军

(1.环境演变与自然灾害教育部重点实验室，北京 100875； 2.北京师范大学地理科学学部，北京 100875; 3.河北省气象灾害防御中心，河北石家庄 050021； 4.应急管理部-教育部减灾与应急管理研究院，北京 100875；5.北京师范大学国家安全与应急管理学院，北京 100875； 6.地表过程与资源生态国家重点实验室，北京 100875)

渤海是北半球纬度最低的结冰海域[1]，渤海及环渤海地区是中国的经济中心之一，已发现胜利油田、辽河油田等丰富的油气资源，且海上运输、渔业养殖等人类经济活动频繁[2-3]。海冰冻结和漂移均会对渤海海运、油气勘探和生产等造成不同程度的影响，形成海冰灾害。截止目前，海冰已导致了多起石油平台坍塌、船舶损坏、航道堵塞等灾害事件发生，造成了严重的经济损失[4-6]。海冰是全球冰雪圈层的重要组成部分，对气候变化较为敏感，研究海冰年际变化、日变化差异，可以精确指示局部区域内的气候变化[7]。海冰淡化因海冰平均盐度远低于海水，有利于降低淡化成本，可以为中国北方地区解决缺水问题提供新的战略水资源[8-10]。海冰资源开发利用的基础之一就是要较准确地估算海冰储量，而海冰面积和厚度时空变化特征研究是海冰资源量估算的基础[11]。自20世纪70年代以来，遥感在海冰监测中发挥着重要作用[12-15]。基于主被动遥感数据可长期、大范围、准实时、动态地监测渤海海冰冰情参数变化[16-17]。因此，基于多源遥感数据研究海冰冰情时空变化特征，对预防和减轻海冰灾害、研究环渤海沿岸地区气候变化、估算海冰资源储量和完善中长期海冰预报统计模型等均有重要意义。但是，遥感卫星由于时空分辨率低、云层遮挡等原因，难获取海冰冰情指标长时序连续日数据，因此，国内外的学者们利用光学遥感例如MODIS和AVHRR提取的海冰数据来建立渤海海冰和气象因子的统计学模型。目前，很少有学者利用GOCI数据来研究长时序的渤海海冰，而GOCI的高时间分辨率和高空间分辨率更适合对渤海海冰进行大范围、近实时监测。此外，有关渤海海冰长达数十年的面积和厚度的日数据集的研究也很少。本文基于AVHRR、MODIS、GOCI三种遥感数据解译出的1988—2018年的渤海海冰面积和厚度数据集分析了渤海海冰冰情指标的变化特征，以及结合气象资料建立统计学模型，补全了1988—2018年渤海海冰面积日数据集。

1 研究区概况

渤海是中国唯一的内海，三面被陆地环绕，也是中纬度典型季风气候区，位于37°～41°N，117°～123°E之间，每年冬季都会结冰，是北半球纬度最低的有冰海域[18]。渤海东北西南向长约300 n mile，东西向长约187 n mile，面积约为77 284 km2，平均水深18 m，最大水深85 m。海冰的形成和分布是区域气象要素和物理海洋要素共同作用的结果，具有明显的区域差异。根据地理位置不同及冰情差异，渤海可分成辽东湾、渤海湾、莱州湾三个海区[19]。每年冬季渤海冰冰期长短、分布范围大小、厚度、表面堆积程度等因素差异较大，年内和年际变化明显，海冰灾害时有发生，造成了一定的区域社会经济等损失。例如2009/2010年冬季，渤海及黄海北部发生了海冰灾害，造成辽宁、河北、天津、山东等沿海三省一市受灾人口6.1万人，船只损毁7 157艘，港口及码头封冻296个，水产养殖受损面积约2.08×105hm2，因灾直接经济损失63.18亿元(灾情数据来源于2010年中国海洋灾害公报)。因此，渤海海冰冰情监测及其时空变化特征研究对于海冰灾害防灾减灾具有重要意义。

图1 研究区地理位置(该图根据审图号为GS(2015)732号的标准地图制作，边界无修改)

2 数据及方法

2.1 数据来源

遥感数据来源于美国国家海洋和大气管理局(NOAA)系列气象卫星上搭载的多光谱通道高分辨率辐射扫描仪(Advanced Very High Resolution Radiometer, AVHRR)1988年1月至2000年3月期间累计298幅，美国宇航局大型空间遥感仪器中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer, MODIS)2000年12月至2015年3月期间累计385幅，以及韩国新一代海洋水色遥感器所涉地球静止海洋彩色成像仪，(Geostationary Ocean Color Imager, GOCI)2015年12月至2018年3月期间累计184幅无云有冰、质量较好的影像。

气象数据来源于中国地面国际交换站气候日数据集，可从中国气象数据服务中心网站(http://data.cma.cn/en)下载，包括渤海沿岸18个气象站点(营口、熊岳、大连、绥中、瓦房店、兴城、黄骅、秦皇岛、乐亭、唐山、天津、塘沽、东营、龙口、惠民、潍坊、长岛、大洼)1988年1月至2018年3月期间平均气温、最低气温和最高气温的观测数据。

2.2 方法

在海冰面积提取和海冰厚度反演之前，需要先对下载下来的遥感图像做预处理，其中包含格式转换、辐射校正、几何校正和大气校正等。

2.2.1 海冰面积提取

对于AVHRR图像，基于分区阈值法将渤海划分成辽东湾、渤海湾、莱州湾、渤海中部四个海洋环境相对较均一的区域，使各子区域内海水反射率小于海冰反射率，根据冰-水反射率分区设定阈值提取海冰面积信息[20-21]。对于MODIS数据和GOCI数据采用基于样本的面向对象分类方法提取海冰面积，分割尺度与海冰最小尺寸相关，设置为20；合并尺度则表示类相似度，设置为70。上述步骤之后，通过目视解译进一步校正，去除明显错分的海冰多边形，最终得到海冰覆盖区域。

2.2.2 海冰厚度估算

海冰面积提取后，采用太阳短波辐射吸收理论的海冰遥感反演模型来反演厚度，并确定不同分区遥感反演模型相关参数，在此基础上得到渤海海冰厚度分布。海冰光学性质的理论模式和实验观测表明一定的厚度范围内，海冰反照率随冰厚的增加而增加[22]，海冰厚度与海冰反照率之间呈现良好指数关系[23]，尤其对于一年冰或冰厚小于1 m的海冰而言，海冰反照率反演模型具有时间和空间上的优势[24-25]。海冰厚度与反照率之间的关系式可表达如下：

α(h)=αmax[1-k·exp(-μαh)]。

(1)

式中：α(h)是与海冰厚度相关的海冰反照率；αmax为无限大冰厚斜对应的反照率；k为与αmax有关的相关系数；μα是关于反照率的衰减系数，h是冰厚(m)。

2.2.3 传感器校正

为了更准确的获取海冰面积和厚度时间序列，本文选择交叉定标进行传感器间校正[26]。利用2000—2011年间的MODIS和AVHRR的数据集建立线性回归模型，然后将获得的模型用于校准AVHRR结果，因为MODIS结果似乎更准确[27]。选择2012—2015年间MODIS和GOCI数据集建立线性回归模型，并对GOCI提取海冰面积进行传感器间校正，选择MODIS数据作为标定基准进而校正GOCI数据，最终将不同分辨率的遥感数据统一成MODIS数据所对应的分辨率。这是由于MODIS数据与AVHRR和GOCI数据序列均有部分时段重叠。通过传感器间校准提高海冰面积和厚度时间序列的一致性和准确性。

2.2.4 气象数据计算

通常将累积冻冰度日(CFDD)定义为一定时间段内每日平均气温和海洋冰冻温度之差的总和(当平均地面气温低于海洋冰冻温度)，累积融冰度日(CTDD)定义为一定时间段内海洋冰冻温度和每日平均气温之差的总和(当平均气温高于海洋冰冻温度)[28]，如式(2)和式(3)所述。

(2)

(3)

式中：Tf是渤海的海洋冰冻温度，参考标准选取-1.8℃[29]；Ta是日平均气温，本文取渤海沿海的18个台站日均气温的平均值；t是冰期天数，ds和de分别为开始和结束的日期。

通常将低于0℃的日平均气温累加值称为负积温[30]，本文将3 d-t℃积温定义为连续的3 d内低于t℃的日平均气温之和。

2.2.5 海冰指标计算

国家海洋局每年以辽东湾、渤海湾和莱州湾三个海区的海冰结冰范围和冰厚作为主要指标，把海冰冰情划分成5个等级：轻冰年(1.0)、常冰年偏轻(2.0)、常冰年(3.0)、常冰年偏重(4.0)、重冰年(5.0)，并在这些等级范围内，细分出常冰年明显偏轻(1.5)、常冰年略偏轻(2.5)、常冰年略偏重(3.5)和常冰年明显偏重(4.5)4个亚等级。1989—2018年间渤海冰情等级来自于国家海洋局每年发布的海洋公报数据；1988—1989年的渤海冰情等级数据来自于文献[31]。

本文将海冰曾结冰频率定义为：

(4)

式中：Pij为第i月第j个栅格的曾结冰频率，ei总为该栅格在1988—2018年i月的遥感影像总数，ei结冰为该栅格在1988—2018年i月的遥感影像中结冰的天数。

冰期指海冰从出现到消融的总天数(一般为12月至次年的3月)。本文以海冰面积达到最大值的日期为分界，将冰期分为冻冰期和融冰期，即冻冰期为出现海冰的第1 d到此分界日的总天数，融冰期为此分界日到海冰消失至无的总天数。

2.2.6 多元非线性回归模型构建

选取与海冰面积相关系数高的四个气象因子建立多元非线性回归模型，利用最小二乘法分别在不同发展阶段拟合出各个阶段的系数，补全遥感影像缺失的日期。

3 结果

3.1 1988-2018渤海海冰时间序列变化特征

3.1.1 年际变化

1988—2018年渤海冰情等级时间序列如图2所示。整体来看渤海冰情等级有明显的年际波动，且2000年之后的年际波动更为明显。刘煜等人[32]的研究表明黄渤海的冰情等级从1950—1980年平均冰情等级的2.9级下降到1980—2010年平均冰情等级的2.3级，降幅为0.6，与全球和东亚地区气候变暖的趋势一致。本文根据资料[16,31]记录的冰情等级，计算了1950—1984年和1985—2018年的平均冰情等级，分别为2.91和2.40，下降了0.51。然而，以冰情等级变化来看，1989—2003年的平均冰情等级比2004—2018年的平均冰情等级上升了0.3，表明2004年后的冰情比之前略有加重且波动明显。1988—2018年间冰情等级最低为2002年的1.0级，属于轻冰年；最高为2001年的4.0级，属于常冰年偏重。1988—2018年间并未出现重冰年，而属于常冰年偏轻(2.0)的年份最多，共有10个年份。

图2 1988—2018年渤海冰情等级时间序列图

图3中可以看到1988—2018年间渤海年平均海冰面积数值主要在4 000～20 000 km2之间，年际变化明显。年均海冰面积均值为8 390.96 km2。其中，最大值出现在2010年，为18 895.59 km2，远大于平均值，这一年发生了1988—2018年间最严重的一次海冰灾害。最小值出现在2015年，为4 533.01 km2。2015年也是冰情等级最低的年份之一，冰情等级为1.0。图3还可以看到1988—2018年间海冰厚度年平均数值主要在5～17 cm之间，有明显的年际变化，年均海冰厚度均值为11.09 cm。其中，年均海冰厚度最大值也出现在2010年，为16.16 cm；最小值出现在2007年，为5.38 cm。2002年冰情等级为1.0，年内平均海冰厚度也非常小，为5.66 cm，接近最小值，也是冰情等级最低的年份之一。

图3 1988—2018年渤海年平均海冰面积和海冰厚度时间序列

3.1.2 年内变化

从图4中可以看出，1988—2018年渤海日平均海冰面积基本呈现为先冻结增加达到年最大值，之后再消融减少，整个过程呈单峰型。日平均海冰面积的峰值出现在1月25日，对应日平均海冰面积为16 278.74 km2。渤海日平均海冰厚度变化也基本呈现为单峰型，即日平均海冰厚度先增加达到峰值，再减小。日平均海冰厚度峰值也出现在1月25日，对应日平均海冰厚度为16.37 cm。从平均趋势来看，海冰冻融过程转折点为1月25日，此前海冰整体呈增加的趋势，此后呈消融趋势，逐渐消融至无。

图4 1988—2018年渤海日平均海冰面积和海冰厚度时间序列

3.2 1988-2018渤海海冰空间分布特征

3.2.1 海冰曾结冰频率空间分布

从图5来看，渤海海冰曾结冰频率空间分布的特征是辽东湾曾结冰的范围最大，其次是渤海湾，最后是莱州湾。1988—2018年冬季渤海曾结冰的范围达到54 544 km2，占渤海总面积的70.58%，其中辽东湾约占结冰范围的63%。冬季3个月中，2月份渤海曾结冰范围最大，达到了49 655 km2，其中辽东湾曾结冰范围约占68%。12月份渤海曾出现海冰面积约占渤海总面积的39.88%，1月份渤海约60.90%的范围以及2月份约64.25%的范围曾出现海冰。

图5 1988—2018年渤海曾结冰频率空间分布

3.2.2 海冰厚度空间分布

从1988—2018年渤海海冰12月、1月、2月和冬季平均厚度空间分布图(图6)来看，渤海海冰冰厚空间分布的特征是离岸越近，厚度越大，辽东湾的冰厚最大，其次是渤海湾，莱州湾的冬季平均海冰厚度是三个海区中最小的。对于辽东湾海区来说，东岸的冬季平均海冰厚度比西岸的更大，冬季平均海冰厚度的最大值出现在营口市和盖州市附近，约为0.24 cm。沿岸的海冰厚度相较于其他区域更厚，且渤海中间区一般没有海冰的分布，海冰厚度平均值为0。从时间上来看，三个海区1月和2月的平均冰厚明显比12月的平均厚度增大，且海冰范围增大、海冰外缘线离岸距离增大。

图6 1988—2018年渤海海冰平均厚度空间分布

3.3 海冰冰情指标对气象因子的响应

海冰的生消与该地区的温度密切相关[33]，本文选取一些相关的气象因子例如：日平均温度、日最低温度、3 d-1.8℃积温、累积冻冰度日等来分析渤海海冰冰情指标与气象因子的关系。本节使用的冻冰期和融冰期的分界点是图4得出的1月25日，将从出现海冰到次年1月25日算做冻冰期，每年的1月25日到海冰消融算做融冰期。

对比表1看出在海冰面积和气象因子的相关性分析中，以12月到次年3月期间海冰面积最大值出现的日期作为分界点，将海冰冻融过程分为冻冰期和融冰期分别与气象因子做相关性分析的相关程度高于直接将冰期与气象因子做相关分析的相关程度。尤其是冻冰期海冰面积与CFDD的相关性为高度相关(r=0.821)，明显高于总冰期海冰面积与CFDD的相关性(r=0.524)；融冰期海冰面积与CTDD的相关性(r=-0.521)也明显高于总冰期海冰面积与CTDD的相关性(r=-0.132)。在冻冰期，与海冰面积相关系数最高的气象因子是CFDD；在融冰期，与海冰面积相关系数最高的气象因子是3 d-1.8℃积温。

表1 海冰面积厚度与气象因子相关系数

在海冰厚度和气象因子的相关性分析中，将海冰冻融过程分为冻冰期和融冰期分别与气象因子做相关性分析的相关程度也高于直接将冰期与气象因子做相关分析的相关程度。尤其是冻冰期海冰厚度与CFDD的相关系数为0.465，明显高于总冰期海冰面积与CFDD的相关系数0.257。在冻冰期，与海冰厚度相关系数最高的气象因子是CFDD；在融冰期，与海冰厚度相关系数最高的气象因子是CTDD。

3.4 海冰面积时间序列数据集补全

多元非线性回归方程：

(5)

式中：Tmin为日最低气温，T为3 d-1.8℃积温，CFDD为累积冻冰度日，CTDD为累积融冰度日。冻冰期系数：a1=-386.164,a2=-132.346,a3=-3.254,a4=135.72,a5=-70.41,a6=2.932,a7=1.072,a8=-0.685,a9=-3.824,a10=0.124,a11=-0.003,a12=-0.004,a13=-683.058。融冰期系数：a1=-207.436,a2=-114.502,a3=-58.967,a4=99.492,a5=109.589,a6=13.155,a7=-0.685a8=-1.387,a9=5.137,a10=0.32,a11=-0.002,a12=0.005,a13=1 541.771。

冻冰期的多元非线性回归方程的R2=0.803；融冰期的多元非线性回归方程的R2=0.689。

从图7中可以看出，2016年2月2日之后的几天和2018年2月20日前后的几天，真实值比预测值偏大的较为明显。2016年和2018年的冰情等级分别是3.0和2.5，比用来回归拟合的26个年份的平均冰情等级2.4大，属于偏重的冰情年份，这或许导致了预测值在此时间段出现低估的现象。2015年1月21日之后的几天，真实值比预测值偏小的较为明显，而2015年的冰情等级为1.0，属于轻冰年，这或许导致了预测值在此时间段出现高估的现象。由于轻冰年和重冰年的海冰面积差异较大，此回归模型在预测重冰年的海冰面积时，可能会一定程度上的低估，在预测轻冰年的时候可能出现一定程度上的高估。如果遥感提取样本数足够多，可以按轻冰年、常冰年和重冰年分别建立回归模型，这会使得预测的结果更接近真实值。目前，遥感提取海冰面积与回归模型计算的海冰面积整体趋势一致，且回归模型具有较好的拟合优度，具有一定的补全和预报价值。