基于逐旬滚动主成分回归分析的渤海海冰预测方法研究

周须文，史印山，井元元，吴国明

（1.河北省气候中心，河北石家庄050021；2.河北省气象科学研究所，河北石家庄050021）

基于逐旬滚动主成分回归分析的渤海海冰预测方法研究

周须文1，史印山1，井元元1，吴国明2

（1.河北省气候中心，河北石家庄050021；2.河北省气象科学研究所，河北石家庄050021）

利用1960—2013年渤海冰情等级资料和500 hPa、1000 hPa NCEP高度场再分析格点资料，采用滑动平均、相关分析、主成分分析等统计方法，找出了与渤海冰情密切相关的大气环流因子，探讨了采用主成分分析对渤海冰情进行定性、定量预测方法。结果表明：在主成分得分散点图上不同冰情年份的落区有一定的规律，根据年份的落点位置可以对冰情等级做出定性预估；采用主成分回归（PCR）方法建立的定量预测模型，具有较好的历史拟合。对2005—2014年的独立样本试报试验，相关系数为0.76，通过了α=0.01的显著性检验。

环流异常；环流因子；主成分分析；海冰预测

1前言

海冰是渤海的主要海洋灾害之一。近年来随着环渤海经济的快速发展，海冰的危害日益凸显，严重的冰情会给涉海企业带来毁灭性的破坏。因此，研究渤海海冰异常成因、提高海冰预测能力，对于保障海上交通运输及各类工程作业的安全，具有非常重要的意义。

20世纪60年代，高纬沿海国家相继开展了海冰观测和研究工作，发布海冰预报。我国自70年代起开展了海冰预报业务。多年来，国内外学者对海冰的预报进行了不懈的探索。其中，美国海军的海冰预报系统[1]主要基于模式和资料的结合，制作5—7 d的短期海冰预报。我国国家海洋环境预报中心根据海洋、水文和气象要素场与冰情之间的统计关系，建立了业务化的海冰预报系统[2]。邓冰等[3]使用了隶属函数、最优气候均态、最优子集回归、逐步回归等方法对渤、黄海海冰的预报方法进行了分析研究；国家海洋环境预报中心研制开发了渤海海冰数值模式[4-5]，建立了海冰数值预报系统。但受海冰观测资料、数值气象场预报精度和时效等因素限制，只能作1—5 d预报。本文通过渤海冰情与大气环流关系的分析，找出影响渤海冰情的异常关键区因子，建立因子与渤海冰情之间的关系，为渤海冰情的中长期预测提供技术支撑。

2资料与方法

海冰资料为1960—2013年冰情等级资料，来自国家海洋环境预报中心、文献[6]和国家海洋局1991—2013年海洋灾害公报。冰情等级的划分是：1.0—1.5为轻冰年，2.0—2.5为偏轻冰年，3.0为常冰年，3.5—4.0为偏重冰年，4.5—5.0为重冰年。海冰资料中，如2010年的冰情等级是指2009—2010年冬季的冰情。

高度场资料为北半球1959—2014a NCEP/ NCAR再分析500 hPa和1000 hPa高度场逐日资料，水平格距为2.5°×2.5°。气温资料为国家气候中心整编的1959—2013年12月—翌年3月逐日平均气温资料，站点选取环渤海沿岸的绥中、营口、乐亭、黄骅4个站。

采用滑动平均、相关分析和主成分分析等统计

方法进行研究。渤海海冰的状况是对中短期气候变化的响应，采用月资料进行分析，为了消除自然月的时间长度与天气过程的时间长度不完全一致，对高度场和气温资料采取了逐旬滚动的月平均处理，这里逐旬滚动月平均高度场（月平均气温）是指对连续三个旬的高度场（逐日平均气温）求平均，如12月3旬的月平均高度场（月平均气温）为12月21日到翌年1月20日高度场（日平均气温）的平均。预测方法流程见图1。

图1 预测方法流程图

3环流异常特征及因子的选取

3.1 环流异常关键区的时空分布

为了找出大气环流对渤海海冰异常有影响的时段和环流因子，对7月1旬到翌年3月1旬逐旬滚动的500 hPa（1000 hPa）月平均高度场与渤海冰情进行相关分析，发现较好相关区的空间分布具有相似性，时间上具有连续性，且具有明确的物理意义，这样的相关区称为异常关键区。

图2逐旬滚动月高度场与渤海冰情等级相关的空间分布

图2 给出了12月1旬—翌年2月2旬的月平均高度场与渤海冰情等级相关的空间分布（图中不同颜色的闭合线代表不同时段的相关区，闭合线中各点的相关系数≥0.35，显著水平α=0.01）。由图2a可以清楚看出，在500 hPa上存在6个异常关键区，分别为高纬欧洲北部和亚洲北部的正相关区，它们与北半球极涡的位置、强度密切相关，极涡代表北方的冷空气，极涡位置、强弱决定了影响渤海区域的冷空气路径和强度；中纬度环流带上北美大陆西部、欧洲大陆和亚洲蒙古一带负相关区，反映了中纬度锋区的位置及西风带槽脊活动对渤海区域气温的影响[7—8]；低纬度太平洋上空负相关区，该区域与西北太平副热带高压相对应，多位学者研究表明[9—14]：西北太平副热带高压指数与渤海海冰关系密切。

由图2b可以看出在1000 hPa有两个异常关键区，一个是欧亚北部正相关区，一个是青藏高原东北部的负相关区。显然欧亚北部的冷空气东移南下影响着渤海海冰的冰情；青藏高原东北部的负相关区是青藏高原大地形对高层冷空气响应的结果，它对渤海冰情也具有明显的指示意义。

为了找出对渤海冰情有影响的大气环流时段，计算了异常关键区的面积和相关程度。异常关键区面积由异常区相关系数绝对值≥0.35的格点数表示，相关程度由相关系数绝对值最大值表示。图3可以看出，500 hPa异常关键区与1000 hPa异常关键区的变化趋势基本一致，异常关键区面积与其相关程度成正相关。大气环流与渤海冰情关系密切的

时段分为两部分，一是海冰出现前期10月2旬到11月1旬，二是海冰同期的12月1旬到翌年2月2旬。同期高度场相关性明显好于前期。

3.2 环流因子的选取

为了简洁、定量、正确的描写影响渤海冰情的大气环流异常状态，把环流异常关键区的状况用环流因子来定量表示，具体做法：找出异常关键区中相关系数绝对值≥0.35格点，计算每个关键区所选格点高度场的平均值，以1981—2010年高度场平均值作为常年值，求其距平，并进行标准化值处理，即形成8个环流因子。表1给出了环流因子对应的异常关键区及所在位置。a1—a6为500 hPa上的环流因子，其中，a5、a6代表极涡的位置和强度变化，a1—a3代表中纬度上槽脊状态，a4代表西北太平洋副热带高压的变化；b1、b2为1000 hPa上的环流因子，b1代表西伯里亚高压的强弱，b2代表高层冷空气在亚洲东移南下的状况。

图3 相关区的面积和最大相关系数时间分布

表1 因子对应的异常关键区及所在位置

4主成分分析

4.1 因子间的相关分析

对所选因子以及海冰等级的相关分析（表2），可以看出：海冰等级与各因子间都有较好的相关性，说明这些因子能够较好的描述海冰异常年的大气环流异常状态。a1与海冰等级的相关性最好，并且a1与其它因子都具有较好的相关，说明a1（500hPa蒙古异常区的环流）是影响渤海冰情最直接、最重要的因子。a5与a6有较好的正相关，表明欧洲北部和亚洲北部的环流对极涡的响应具有较好的一致性；a5、a6与a1、a2、a3、b1有较好的相关，表明极涡的变化直接影响着中纬度槽脊，地面上西伯利亚高压也会有明显的表现。b2只与a1有较好的相关，表明1000 hPa青藏高原东北部负相关区是对亚洲中纬度高空冷空气的响应。a4与各因子的相关性较差，表明低纬太平洋区域的环流受中高纬环流的影响较小。

表2 8个因子及海冰等级的相关性

4.2 主成分分析

由于选取的8个因子间大多存在一定的线性相关性，因此利用主成分分析从8个因子的观测样本矩阵中提取主成分，它们是原变量的线性组合且相互正交。这样可以保留大部分原有信息，消除复共线性。以12月3旬的月高度场为例，选取1960—2013年8个因子的资料作为观测样本矩阵，对其进行主成分分析，结果为：第一主成分的贡献率为53%，第二主成分的贡献率为15%，第三主成分的贡献率为12%。下面给出了前2个主成分的表达式：

由式（1）看出：第一主成分在每个因子上有相近的载荷，说明每个因子对第一主成分的得分（y1）重要性相差不多，并且正相关因子a5、a6和b1是负载荷，负相关因子a1、a2、a3、a4和b2是正载荷，因子与它们对应载荷的乘积具有一致性，因此第一主成分的贡献率最大；式（2）看到，第二主成分各因子的载荷相差较大，并且负相关因子a1、b2是负载荷，与其他因子不一致，第二主成分的贡献率明显小于第一主成分。

5渤海冰情逐旬滚动预测方法

5.1 渤海冰情的定性预估

图4a是12月3旬的主成分得分散点图（y1为第一主成分得分，y2为第二主成分得分），可以看出：在y1＜0，y2＜0.5的区域中多为重冰年，12 a重冰年中10 a出现在该区域，频率为71.4%；0≤y1＜1的区域中，多为常冰年，出现的频率为80%；在1≤y1的区域中多为轻冰年。图4b是2月2旬的主成分得分散点图，可以看出：在斜线的左下部多为重冰年，右上部多为轻冰年。显然，前两个主成分包含了大部分的原始信息，根据年份的落点位置可以给出冰情的定性预估。

根据渤海沿岸气象站逐旬滚动的月平均气温变化（图5），可以给出图4a中特殊年份落点位置的解释。1964年为重冰年，出现在轻冰年（1≤y1）区域中，主要原因是渤海沿岸1964年冬季是典型的前冬暖、后冬冷的年份，冰情较重是因后冬气温明显偏低造成的。1961年为轻冰年，出现在重冰年区域中，1961年渤海冰情不重是由于前冬气温偏低，后冬气温偏高，负积温小，不利于海冰的发展。2010年冬季12月和2月气温偏低，负积温时间长，冰情较重，但由于1月气温偏高，因而2010年偏离了重冰年落区。1969年是冰情最严重的年份，整个冬季气温明显偏低，落点在重冰年区且远离轻冰年。这些特殊年份在图4b中的落区都是正确的。以上分析可知，海冰生成和发展在时段上有较大差别，海冰的冰情是一段时间累积结果，年份点的落区和位置对冰情具有一定的指示意义，只要根据逐旬滚动的主成分得分散点图仔细分析，就可以更准确地制作出渤海冰情的定性预估。

5.2 渤海冰情的定量预测

建立逐旬滚动的渤海冰情的定量预测模型，采

用的是主成分回归（PCR）预测方法。按累计方差贡献率达80%左右作为选取主成分的标准[15]，选取前几项方差较大的主成分，略去方差较小的一些主成分，建立多元回归方程。以12月3旬的主成分分析为例，前3个主成分的累计方差贡献率为80%，因此选取前3个主成分进行多元回归计算，回归方程为

式中：y为海冰等级预报值，y1、y2、y3分别为第一主成分得分、第二主成分得分、第三主成分得分。经F检验，显著水平α=0.001时，回归方程是显著的，并具有较好的历史拟合。

采用NCEP实况资料和derf 1.0模式（预报未来4个旬）资料分别进行了独立样本的试报检验，对2005—2014年进行的独立样本试报检验结果为，NCEP实况资料的预测结果与实况的相关系数为0.76，显著水平为α=0.01；derf 1.0模式资料的预测结果与实况的相关系数为0.6，显著水平为α=0.05。

在实际业务中为了提高预测准确率，一是采用NCEP实况资料，应用前期的回归方程和主成分散点图作出预测，二是采用数值预报产品统计释用的PP方法，应用同期的回归方程和主成分散点图作出预测，最后对两种资料的预测结果进行综合分析，得出预测结论；并随着时间的推移、资料的更新，逐旬滚动作出预测，及时修订预测结论。

图4 主成分得分散点图（蓝色*代表重冰年，绿色*代表常冰年，红色•代表轻冰年）

图5 特殊年份逐旬滚动月气温距平分布

6结论与讨论

（1）渤海海冰与前期（10月2旬—11月1旬）和同期（12月1旬—2月1旬）的月高度场有较好的相关性。500 hPa和1000 hPa上的环流异常关键区具有一致的变化趋势；

（2）海冰等级与环流因子具有较好的相关性，这些因子可以较好描述海冰异常年大气环流的异常状态；

（3）在主成分得分散点图上，不同冰情的年份落点位置有一定的规律，对定性预测冰情有指示意义；

（4）采用主成分回归建立的冰情定量预测模

型，对中长期冰情的定量预测提供了理论依据；

（5）该方法在预测准确率以及预测效果稳定性等方面还会存在一定缺陷，因此在实际业务应用中采用实况资料和模式资料同时应用，预测结论综合考虑，并逐旬滚动作出预测，及时修订预测结论。

[1]高文洋,郭婷婷,王长甫,等.美国海军海冰预报系统[J].海洋预报,2009,26(2):103-109.

[2]李春花,刘煜,白珊,等.海冰预报综合系统简介[J].冰川冻土, 2003,25(增刊2):349-351.

[3]邓冰,佟凯,张学宏,等.渤、黄海海冰预报方法研究[J].海洋预报,2004,21(3):15-21.

[4]吴辉碇.海冰客观分型、四维同化和数值预报产品的研制(1991—1995)[R].北京:国家海洋局国家海洋环境预报研究中心,1995.

[5]白珊,吴辉碇.渤海的海冰数值预报[J].气象学报,1998,56(2): 139-153.

[6]杨华庭,田素珍,叶琳,等.中国海洋灾害四十年资料汇编(1949—1990)[M].北京:海洋出版社,1993:171-181.

[7]毛睿,龚道溢,房巧敏.冬季东亚中纬度西风急流对我国气候的影响[J].应用气象学报,2007,18(2):137-146.

[8]杨小波,王永光,梁潇云.11月气候异常型及前期环流信号[J].应用气象学报,2011,22(3):275-282.

[9]耿淑琴,王旭.2001—2002年冬季黄渤海天气气候特征及对渤海海冰的影响[J].海洋预报,2002,19(4):38-47.

[10]耿淑琴,王咏亮,白珊.2000—2001年冬季北半球大气环流变化对渤海气候及海冰的影响[J].海洋预报,2001,18(3):25-33.

[11]李彦青,苏洁,汪洋,等.渤海海冰外缘线候平均离岸距离的变化及其关键影响因子[J].中国海洋大学学报,2013,43(7):7-16.

[12]白珊,刘钦政,吴辉碇,等.渤海、北黄海海冰与气候变化的关系[J].海洋学报,2001,23(5):33-41.

[13]刘钦政,黄嘉佑,白珊,等.渤海冬季海冰气候变异的成因分析[J].海洋学报,2004,26(2):11-19.

[14]李剑,黄嘉佑,刘钦政.黄、渤海海冰长期变化特征分析[J].海洋预报,2005,22(2):22-32.

[15]丁德文.工程海冰学概论[M].北京:海洋出版社,1999.

Research on predictionmethod of the sea ice in the BohaiSea based on ten-days rolling principal component regression analysis

ZHOU Xu-wen1,SHIYin-shan1,JING Yuan-yuan1,WU Guo-m ing2
(1.HebeiClimate Center,Shijiazhuang 050021China;2.Hebei Institute ofMeteorological Science,Shijiazhuang 050021 China)

Based on the sea-ice grade data in the Bohai Sea from 1960 to 2013 and the NCEP/NCAR reanalysis data at 500 hPa and 1000 hPa,the atmospheric circulation factors associating closely to sea-ice grade were selected by using the statisticalmethods ofmoving average,the correlation analysis and the principal component analysis(PCA),and the qualitative and quantitative prediction method for sea-ice conditions in the Bohai Sea was investigated.The results show that there is some rule of falling areas in different sea-ice grades at the PCA scatter chart,thequalitative prediction of sea-icegrades in the BohaiSea can beestimated according to the falling areas of the year.The quantitative forecastmodel is established by the principal component regression method (PCR)and has a good historical fitting rate.The correlation coefficient is 0.76 w ith significant at the confidence levelof0.01 for the independentsamples test from 2005 to 2014.

circulation anomalies;circulation factors;PCA;sea ice prediction

P731.15

A

1003-0239（2015）06-0074-06

2015-02-28

河北省科技厅渤海海冰发生规律及预测技术研究（14273501D）

周须文（1963-），女，高级工程师，大学本科，从事短期气候预测工作。E-mail：hbszxwgood@163.com

10.11737/j.issn.1003-0239.2015.06.010