1 引言

海平面变化是气候变化的重要参数之一，受大气圈、水圈、岩石圈以及人类活动等多要素的综合作用，时空特性鲜明，变化机理复杂^{[1, 2, 3]}，海平面变化预测是海平面变化研究乃至全球研究中的热点也是难点问题^{[4, 5, 6, 7]}。目前国际上关于器测时期海平面变化预测可分为数理统计预测和动力模型预测两种方法。对于大尺度如全球性海平面变化，可利用动力模型根据影响海平面变化的驱动因子（如全球平均气温）及作用机制进行过程清晰的长期预测^[8]。然而对于区域影响要素复杂多样、物理过程及作用机制均不够清晰的区域性海平面变化研究，动力模型预测精度有待提高。数理统计预测是基于海平面过去时间变化序列进行建模和外推，具有资料可靠，精度高的优点，是现阶段区域海平面变化预测的有效途径^[9]。

早期统计预测通常假定海平面过去变化特征在未来一段时期保持一致，线性拟合其趋势加以外推，是最简单可行的方法^[10]，但预测结果也相对较粗略。随着傅利叶变换、小波分析等时频分析方法被逐渐应用于海平面变化研究，验潮站序列被分解为确定性趋势项、周期项以及剩余随机序列^{[11, 12]}。预测时将趋势项及周期项外推，并对剩余随机序列进行线性/非线性拟合预测^[13]。近年，针对傅里叶变化基于固定时频窗以及小波分析的基函数选择难度大等问题，Huang等^[14]提出新的时频分析方法——经验模态分解（Empirical Mode Decomposition，简称EMD），依据各周期信号振幅频率特征不同将原始数据分解为一系列具有不同特征尺度的序列，即本征模函数（Intrinsic Mode Function，简称IMF）。EMD方法提取的信号具有很高的信噪比，目前已被广泛应用于地学、水文、气象等领域^{[15, 16, 17]}。王卫强等^[18]利用EMD方法对中国沿岸海平面变化进行研究，并与灰色系统方法和随机动态预测模型比较，认为EMD方法不但能够得到较精确的年变率，而且能够分离出各种天文周期，较好地解决了资料序列的非平稳和非线性问题。集合经验模态分解（Ensemble Empirical Mode Decomposition，简称EEMD）是基于EMD方法的改进，利用噪声辅助减弱模态混叠现象的影响从而改善分解效果^[19]。然而傅利叶变换和EEMD等时频解析方法始终是对现状海平面变化的分解研究，以此直接外推生成未来时间变化序列其精度有待提高。

另一方面，随着时间序列预测模型日益发展成熟，神经网络作为一门新兴边缘学科，因能较好解决非平稳序列间的非线性映射被广泛应用于地学、海洋、大气等领域^{[20, 21, 22, 23, 24]}。BP（Back Propagation）神经网络是其中应用最为广泛、发展最为成熟的模型之一^{[25, 26, 27, 28]}。但神经网络基于纯数据驱动预测，建模过程中无法考虑时间序列本身可能存在的物理意义及规律，无法解析海平面变化的过程及特征。因此，建立一种新的思路，打破海平面时间序列的时频分析方法及纯数据驱动预测方法的局限性，在科学解析过去海平面变化过程及特征的基础上，准确预估其变化规律是当前国际关于区域海平面变化研究需要解决的重要科学问题^{[29, 30]}。

珠江三角洲位于中国广东沿海，区域内人口、财富密集，受海平面变化威胁严重^{[31, 32, 33]}。三角洲内海平面在全球及区域尺度要素等共同影响下表现出相当的复杂性和不确定性，国内外已开展大量研究并取得一定成果^{[34, 35, 36, 37]}。黄镇国等^[10]基于香港及澳门站验潮序列预测2030年珠江三角洲海平面将上升30cm； 冯伟等^[38]利用卫星测高技术研究发现1993～2009 年南海海平面变化为 5.5±0.7mm/a，明显高于同期全球平均海平面变化趋势（3.3±0.4mm/a）。然而以往关于海平面变化预测的研究通常利用现状时间序列进行线性拟合外推或者基于时频分析方法进行周期外推，研究精度有待进一步提高^[39]。因此为明晰海平面变化过程并科学预测其趋势，本研究拟采用EEMD结合BP神经网络建立EEMDBP预测模型，以研究区内标准海洋站闸坡站验潮数据为例，对月均海平面高度序列进行拟合预测。 一方面分析海平面变化的规律，另一方面探讨提高海平面变化预测精度的思路，为区域海平面长期变化预测研究提供参考借鉴。

本研究收集全球海平面监测网络（Global Sea Level Observing System，简称GLOSS）标准站闸坡站1959～2011年海平面月均变化序列资料作为实证研究对象，对海平面长期变化进行预测研究。闸坡站位于中国南部广东沿海（21°35′N，111°50′E； 图1），是中国南海验潮站中历史最长、记录最完整的验潮站序列^[40]。本文采用三次样条插值方法补齐序列中的缺测数据（1～2个月）^[41]。

图1(Fig.1)

图1 闸坡站位置图 Fig.1 Location of Zhapo tide gauge

采用EEMD-BP模型进行海平面变化趋势预测主要包括时间序列的分解和预测两部分。 首先对原始海平面时间序列进行EEMD分解，提取海平面变化过程具有物理含义的信号并降低噪音影响，得到不同尺度周期或趋势(IMFs)，实现序列的平稳化； 继而采用BP神经网络对分解后的各IMFs分别建立预测模型研究其未来变化趋势； 最后将各IMFs的预测结果重构，得到原始序列的预测值。

采用EEMD提取海平面变化序列中的周期项和趋势项，并经显著性检验，确定海平面变化隐含的周期及最终趋势，具体方法如下^{[42, 43]}: 首先在原始序列中K次加入均值为0，幅值标准差为常数的高斯白噪声n_i(i=1，2，…，K)得到时间序列f_i(t)。然后对f_i(t)进行EMD分解: 找出序列中所有极大值和极小值点，分别用三次样条函数合成上下包络线，并得到平均包络线m_i1，将f_i(t)减去m_i1得到去掉低频趋势的h_i1，多次重复上述过程，使平均包络线趋于零从而得到第一个IMF分量c_i1(t)，代表原始序列中最高频的分量，有:

对剩余序列r_i1(t)重复上述分解过程，直到j次分解后所得剩余分量满足给定的结束标准^[14]，分解结束。

则原始序列可表示为:

根据不相关随机序列的统计值为0的原理^[44]，将上述步骤对应的IMF进行总体平均运算，以消除多次加入高斯白噪声对真实IMF的影响。

BP神经网络的预测思路是认为未来值与现状值之间隐含某种函数关系，通过建立神经网络来拟合这种关系并推导出未来值^[45]。典型的BP网络包括输入层、隐含层(中间层)和输出层(见 图2)，建模过程由信息的正向传播和误差的反向传播两部分组成。输入层将接收到的输入信息传递到隐含层，将信息处理变换后传递到输出层输出处理结果。当实际输出与预期输出不符时，进入误差的反向传播阶段。误差通过输出层，按误差梯度下降的方式修正各层权值，向隐含层、输入层逐层反传。神经网络通过反复学习训练，进行多次信息正向传播和误差反向传播，不断调整各层权值，直至输出误差减小到预设值，或者达到预设的学习次数。有关BP神经网络的基本原理已在众多文献中介绍^{[46, 47]}，本文不再赘述。

图2(Fig.2)

图2 神经网络模型示意图 Fig.2 Structure of BP neural network

闸坡站验潮资料经过EEMD分解生成9个IMFs，如下 图3所示，IMF1～IMF9依次对应8个周期项和1个长期趋势项。周期项分别包括3个月、7个月、13个月、26个月、49个月、101个月、220个月和314个月周期，基本对应常用理论天文周期^[48]。其中，前3个周期分别对应季节，半年和1年的变化周期。26个月变化周期为准两年大气振荡，49个月周期与3～7年ENSO周期变化相符。101个月周期为准9年的半交点周期，与黄镇国等^[49]关于广东沿岸珠江三角洲验潮记录的周期分析一致； 而220个月周期可能表现的是潮汐变化周期，与天文学上黄白交点运行的18.61年周期^[50]相近。

图3(Fig.3)

图3 闸坡站验潮资料的EEMD分解结果 Fig.3 Results of the EEMD of Zhapo tide gauge records

为检验IMFs的信息的有效性，对分解结果进行白噪音检验，如下 图4所示，IMF1～3、IMF5～7均通过了置信度为95 ％ 的显著性检验，说明这6个IMFs具有统计学意义，它们分别对应的季节变化、半年变化、年变化、ENSO事件、准9年半交点周期及18.61天文年变化等6个周期项变化显著。未通过检验的其余2个IMFs则被认为是噪音序列，在实际预测中应予以剔除。从各IMF的方差贡献率来看，季节、半年和年周期分别占海平面变化的30.71 ％ ，42.88 ％ 以及20.05 ％ ，说明年内变化对中国南海平面变化影响较大，与前人^[51]研究关于中国海平面变化以年内周期为主一致。同时根据趋势项IMF9的线性拟合显示，1959～2011中国南海近53年年变化率为1.6mm/a，与陈特固等^[52]关于中国南海海域珠江三角洲海平面变化研究较一致。结果显示EEMD能较完整分解出海平面变化的常规周期项，提取的长期变化趋势合理，分解效果较好。

图4(Fig.4)

图4 EEMD分解结果显著性检验 Fig.4 Significance test of IMFs from the EEMD of Zhapo tide gauge records

考虑到海平面变化具有强烈的季节变化特征，本研究采用对应月时间序列组织数据建模，并通过数据序列滑动平移构建样本。原始资料01/1959～12/2011共636个月海平面高度资料被分为两部分，其中前576个资料(01/1959～12/2006)作为训练样本，后60个资料(01/2007～12/2011)作为检验样本^[53]。经过反复调试，建立隐含层为1，输入层神经元为30，输出层为1(如输入层为1959～1988年每年1月海平面高度，输出层则为第31年1989年的1月海平面高度)的神经网络预测模型^[54]。训练过程中逐步调整隐含层的节点数目，得到最佳隐含层节点数为9。在模型的隐含层中选用S型正切函数(tansig)作为输入层至隐含层的传递函数； 选用线性函数(purelin)作为隐含层与输出层间的传递函数； learngdm函数作为学习函数； trainlm函数作为训练函数以及mes函数作为性能函数。研究将学习率确定为0.02，并采用0.003作为误差边界值，即当误差值小于0.003时，认为学习完成，停止计算，输出结果，其他参数采用系数默认值。关于BP神经网络的建模预测过程在Matlab(R 2012a)中完成。

EEMD-BP建模预测即将分解所得IMF1～9分别按上述步骤进行建模预测。IMF1～9预测值与样本资料对比如下 图5所示。从 图5中可以看出，高频分量IMF1( 图5a)的拟合效果最差，但序列整体振幅及走向与样本变化较一致。由IMF2至IMF9，随着各IMF变化趋于平稳，预测效果逐渐变好。说明神经网络虽然适用于非平稳非线性序列预测，但当序列趋向平稳时，预测精度也随之上升。考虑到研究海平面长期变化趋势时，通常应先去除季节等年内高频振荡信号^[55]，因此，IMF1拟合效果不甚理想等对于长期如年代际海平面变化预测研究影响不大。

图5(Fig.5)

图5 EEMD分解各IMF预测值与样本对比 Fig.5 Comparison between the prediction and samples of IMFs from EEMD

为全面评价EEMD-BP模型预测精度，本文同时采用BP神经网络直接建模方法进行对比预测，即不经过EEMD分解，直接利用BP神经网络对闸坡站原始序列进行预测。为避免因降噪处理等造成与原始观测数据比较时产生偏差，进行预测精度对比时将全部IMF1～9预测结果叠加重构，形成EEMD-BP模型预测的最终结果。BP神经网络与EEMD-BP模型的预测结果见 图6，其中蓝色曲线代表观测值，红色曲线代表EEMD-BP模型预测结果，绿色曲线代表BP神经网络预测结果。 图6中可以看出3条曲线的整体趋势较一致，两种预测方法的拟合精度都较高。其中绿色曲线更易出现极值偏差，如在2007年9月，BP神经网络的预测值明显高于观测样本，而EEMD-BP预测结果则与观测样本较接近。另外在2009年11月及2010年11月，绿色曲线都明显偏离红蓝曲线。一方面，红色与蓝色曲线的拟合较绿色与蓝色曲线拟合的整体效果更好，即采用EEMD-BP模型预测比单独使用BP神经网络预测的效果好； 另一方面，BP神经网络预测过程中与原始观测值偏离较大处，即预测误差偏大的点较易出现在时间序列变化的极值点上，如图6中2008年8～10月及2011年8～10月。由此可知，尽管建模过程已通过学习过去变化规律对未来趋势进行模拟预测，但因海平面高度复杂多变，存在大量突变及极值点，相对于长期变化的趋势预测研究，对极值的预测处理相对要差一些，可另开展极值频率、重现期等预测研究。

图6(Fig.6)

图6 BP神经网络预测及EEMD-BP模型预测结果比较 Fig.6 Comparison of the predictions from BP and EEMD-BP models

为进一步分析两种方法的拟合精度及预测效果，本文采用相关系数(R)、平均误差(Mean Error，简称ME)和均方根误差(Root Mean Square Error，简称RMSE)对预测结果进行评价，见 表1。可以看出，两组预测值与检验样本均有较高的相关系数(R>0.75)，较低的均方根误差(RMSE<37mm)及平均误差(ME<4mm)。Cheng等^[56]利用验潮站数据及卫星资料预测英国沿海海平面变化，预测值与观测值最小RMSE为43.9mm，本文整体预测效果较理想。EEMD-BP模型预测结果的相关系数高于BP神经网络模型，而均方根误差及平均误差绝对值则低于BP神经网络模型，总体表现优于BP神经网络。因此，采用EEMD-BP模型预测，首先对数据进行平稳化及降噪处理，依据信号的不同变化特征进行模拟预测优于直接使用BP神经网络对原始数据进行预测。

表1(Table 1)

表1 BP神经网络预测与EEMD-BP神经网络预测误差统计比较Table 1 Error statistics of the predictions from BP and EEMD-BP models R RMSE/mm ME/mm BP 0.76 36.74 -3.46 EEMD-BP 0.89 28.16 2.31

表1 BP神经网络预测与EEMD-BP神经网络预测误差统计比较 Table 1 Error statistics of the predictions from BP and EEMD-BP models

考虑数据的有效性，本研究采用EEMD-BP模型，选用IMF2～3、IMF5～7及趋势项IMF9对珠江三角洲未来2050年海平面变化进行预测，结果如 图7所示。其中，蓝色为闸坡站1959～2011年海平面高度观测值，线性拟合结果显示该时期珠江三角洲海平面变化呈现上升趋势，速率约为2.22mm/a。 图7中红色曲线为预测值，线性拟合显示至2050年，上升速率将达2.72mm/a。

图7(Fig.7)

图7 基于EEMD-BP模型的珠江三角洲1959～2050年海平面变化预测 Fig.7 Past sea level observation and sea level prediction from 1959 to 2050 in the Pearl River delta by EEMD-BP

根据《2010年中国海平面公报》^[57]显示，近30年广东沿海海平面以2.5mm/a速率上升，并预测未来30年将上升75～145mm。本文以闸坡站为代表站点研究珠江三角洲海平面变化，与以上研究相近，预测结果可靠。

本文利用珠江三角洲闸坡站1959～2011年月均海平面高度观测资料，结合EEMD及BP神经网络建立EEMD-BP预测模型进行海平面长期变化趋势预测研究，并将预测结果与直接使用BP神经网络预测对比，分析研究模型预测精度及适用性，在此基础上进行了珠江三角洲2050年海平面变化预测。研究结果表明:

(1)BP神经网络技术因其无需先验知识等优势，可以从纯数据驱动的出发，较好地对海平面高度时间序列进行学习训练、预测，适应海平面变化这类多因素、非平稳和非线性时间序列预测研究。但对于极值模拟预测能力有限，可另行开展极值、重现期等研究。

(2)EEMD-BP模型首先对原序列进行EEMD分解，得到具有不同时频特征的周期项及趋势项之后，再利用BP神经网络对各IMFs建立预测模型，预测精度高于直接使用BP神经网络预测。对各IMFs预测结果可看出，序列越趋于平稳，BP神经网络预测效果越好。EEMD-BP在对数据进行平稳化及降噪处理的基础上完成预测，模型构建合理可行，模型精度显著提高。

(3)闸坡站1959～2011年海平面高度时间序列经EEMD分解得到季节变化、半年变化、年变化、ENSO事件、准9年半交点周期及18.61天文年变化等6项周期项。其中，前3个周期的方差贡献率总和达93.64 ％ ，说明珠江三角洲海平面变化以季节、半年及年等年内周期变化为主。

(4)闸坡站1959～2011年海平面高度资料线性拟合结果显示该时段珠江三角洲以2.22mm/a速率上升，EEMD-BP模型预测至2050海平面上升速率增大至2.72mm/a，有加速上升趋势。