0 引言

由全球变暖和人类活动引起的全球海平面上升，已经成为不争的事实.政府间气候变化专门委员会(IPCC)在最新的气候变化评估报告中指出(Pachauri and Meyer, 2014)，海平面在1993—2010年间的上升速度是1901—2010年的2倍，全球66%地区的海平面将比1986—2005年间高出0.29~0.82 m.海平面的持续上升对沿海地区的生态环境和社会经济都造成了不可逆转的破坏(Evadzi et al., 2018; Jayanth et al., 2018).因此，合理的预测海平面的变化，对于人口日益增长、资源不断减少的沿海国家和地区应对全球海平面上升具有一定的科学和现实意义.

目前，海平面的预测方法主要分为两种，气候模型预测(Collins and Allen, 2002; Luo et al., 2011)和数理统计预测(Chambers et al., 2002; Grinsted et al., 2010).气候模型主要是对大尺度的全球海平面变化进行长期预测(Raper and Braithwaite, 2006).但在全球海平面持续上升的大趋势下，人们更为关心的还是自身区域海平面的变化(Yuan et al., 2009; 郭金运等, 2015; Guo et al., 2016).数理统计预测是通过分析海平面变化的时间序列预测海平面变化的，具有资料可靠、针对性强、精度高的优点，是当前区域海平面变化预测的有效方法(Menendez et al., 2009).常规数理统计预测是通过简单的线性拟合，对海平面变化趋势进行求取和外推，预测结果较为粗略，精度不高.随着随机动态模型的建立，将海平面变化分为确定性趋势项、确定性周期项、剩余随机序列和白噪声序列四部分，其中，白噪声序列影响较小，可以忽略不计.通过确定模型中确定项和随机项的具体表达式，实现对原始序列的拟合和外推.其中，趋势项可通过线性拟合或滑动平均提取，但随着对海平面非线性变化趋势认知的加深，研究者们更多的是利用多变量逐步回归、经验正交函数(Church et al., 2004)、灰色模型(Liang, 2009)等方法提取海平面变化的非线性趋势.周期是对海平面变化影响最大的一项，可通过小波分析、谱分析、经验模态分解等方法提取.剩余随机序列可利用自回归(Auto Regression, AR)模型或者ARMA模型进行拟合.根据原始时间序列数据的特点，选择合适的预测方法进行组合，预测效果会更理想(Niedzielski and Kosek, 2009; 段晓峰等, 2014).Rahmstorf(2007)提出一种半经验模型预测21世纪全球海平变化，Grinsted等(2010)利用蒙特卡洛反演方法也对21世纪海平面的变化进行了研究和预测，两者的预测结果均与IPCC AR4报告结果差异较大，引起了社会广泛关注.因此，基于时间序列数据及数理统计模型预测海平面变化领域，还需要海平面变化研究者对其进行更深入的研究.

SSA作为一种数字信号处理技术(Vautard et al., 1992)，不仅能从时间序列中提取观测数据的非线性趋势，且不受正弦波假定的约束，能够稳定识别和强化周期信号(Guo et al., 2018; 周茂盛等, 2018).SSA作为经验正交函数展开的一种扩展形式，能够提取时间序列中主成分的特性，并利用加权相关分析(w-correlation)方法(Hassani, 2007)选择合适阶次对序列进行重构和外推，特别适合分析和预测有周期震荡的时间序列数据(郭金运等, 2018).

ARMA是一种应用广泛、精度较高的线性平稳时间序列预测模型，只要确定有限个参数的值，模型就能被确定(Valipour et al., 2013)，被广泛的应用于金融学、测量学、海洋学等领域相关时间序列数据的预报(Krankowski et al., 2005; Shen et al., 2017, 2018).该模型适用于平稳的时间序列，若要对非平稳时间序列进行分析和预测，需要通过差分运算，对非平稳时间序列进行平稳化处理.

日本是位于西北太平洋的一个岛国，海平面的变化直接影响当地地区的生态环境和经济发展.且因日本处于板块交界地带，地壳运动频繁，很难获取并预测日本近海海平面绝对变化.本文拟沿日本海岸线选取28个验潮站，结合验潮站附近的GPS站，分别利用SSA和SSA+ARMA方法分析1994—2013年20年的验潮数据和2003—2013年10年的GPS时间序列数据，并预测2014—2018年5年的海平面绝对变化，与同时段实测数据比对，证明预测方法的可行性.然后，利用1999—2018年的验潮数据和2003—2018年的GPS数据，基于SSA+ARMA方法预测日本沿岸2019—2023年5年的海平面绝对变化.

1 研究区域和数据 1.1 研究区域

日本位于太平洋板块和亚欧板块的交界地带，地处太平洋西北和亚洲东部，主要由北海道、本州、四国、九州等岛屿组成，且鄂霍次克板块和阿穆尔板块在本州岛中部交汇(Bird, 2003).因此日本岛火山、地震活动频繁，地壳垂直运动较为复杂.

日本海岸线十分复杂，全长约33889 km.日本沿海受多种暖流及寒流影响，其中日本暖流自南向北影响日本东部，里曼寒流和对马暖流在日本西部沿海相遇，千岛寒流影响日本东北部.日本四季分明，全年降水量较多，尤其是6、7月份.日本沿岸的海平面变化由于受到黑潮入侵、日本海内环流、厄尔尼诺及相关动力机制、热力机制的影响，存在长、短不同的周期成分，变化特征较为复杂(Ishii et al., 2006).

1.2 验潮数据及预处理

日本沿岸的平均海平面永久服务(Permanent Service for Mean Sea Level, PSMSL)网站约35个验潮站有1994—2013年的连续月数据，但仅28个验潮站有2003—2018年16年连续的GPS联测数据.考虑到验潮数据与GPS数据在时间和空间上的一致性，本文选择了图 1的28个验潮站进行日本沿岸近海绝对海平面变化的分析和预测.验潮月数据引自PSMSL网站公布的修正地方基准海面(Revised Local Reference, RLR))数据(Holgate et al., 2013).

从PSMSL网站下载的验潮月数据的缺失值皆以极值-99999填充，28个验潮站的数据缺失率为1.547%.本文利用SSA迭代插值补全数据.SSA迭代插值(Shen et al., 2018; 周茂盛等, 2018)的具体方法是：首先用0替换验潮数据中的缺失值，将补全后的时间序列中心化处理；然后对中心化后的时间序列做SSA分解，利用w-correlation确定前n阶为时间序列的主成分，时间序列中缺失数值用第一个重构成分中的值代替，反复重复此过程，当两次插补的数据残差RMS小于0.001 m时，结束此插补过程；继续对第2个重构成分重复上述过程，直至第n个重构成分插补完成.其中，验潮站Ito Ⅱ 2018年5月份的验潮数据被标记为验潮数据集的峰值，比相邻两个月的验潮数据高出1.3 m左右，具体原因未给出说明.为了不影响预测精度的评定，本文将Ito Ⅱ 2018年5月份的验潮数据赋值为-99999，利用SSA进行迭代插值作为此月的验潮数据.

1.3 GPS数据及预处理

本文使用的与验潮站联测的GPS坐标时间序列是由The University of La Rochelle(ULR)利用GAMIT/GLOBK软件解算，GPS时间序列数据(天解)可从SONEL(Systèmed′ Observation duNiveau des Eaux Littorales)网站(https://www.sonel.org)下载.GPS数据的参考椭球为GRS80，参考框架为ITRF08(牛余朋等, 2019)，各测站时间序列的平均解算精度为6.65 mm，且GPS时间序列数据对地震等突发情况做了相应的数据改正(Santamaría-Gómez et al., 2017)，符合本文的研究要求.ULR坐标时间序列的解算策略(Santamaría-Gómez et al., 2012)如表 1所示.

对于时间序列中存在的粗差，以3倍标准差(3σ)准则确定阈值，超过阈值的被认为是粗差并剔除，并对剔除的数据利用上述SSA迭代插值方法补全.因为验潮数据是月均海平面数据，为了保持时间尺度的一致性，我们将GPS每个月的天数据相加取平均作为这个月的GPS月数据，利用GPS月数据分析和预测地壳垂直变化.

2 原理与方法

SSA可以对含有噪声的时间序列数据进行非线性趋势提取、周期判别、降噪处理、数据重构和预测，能够从未知物理本质的数据中提取尽可能多的可靠信息(Elsner, 2002).本文利用SSA分别对验潮数据和GPS数据的主成分进行中长期预报，并将扣除主成分的剩余序列当做随机序列，利用ARMA(p, q)模型对其进行外推预测，SSA主成分预测结果加上ARMA外推结果即为最后SSA+ARMA的预测结果.

2.1 SSA预测模型

将各验潮站1994—2013年20年的验潮月数据作为原始序列，分别对各验潮的验潮数据进行SSA预测.各验潮站时间序列长度N=240.设窗口长度为整数，令K=N-M+1，构建M×K阶的时滞矩阵X

(1)

其中X_i=(x_i, …, x_i+M－1)^T，时滞矩阵X的各行、列都是原始时间序列的子序列，其滞后协方差T_x是一个M×M的对称矩阵

(2)

式中t_i(0≤i≤M－1)为原始序列迟后i的自协方差.T_x的特征值为λ_M，对应的特征向量为E_M.设d=rank(X)，，则时滞矩阵X可表示为

(3)

其中初等矩阵是X_r的奇异值，是前j个X_i的贡献率.

通过对角平均化，将初等矩阵X_r转换为长度为N的多个重建向量(reconstruction component, RC)，原始序列等于所有RC之和.令矩阵Z=X_r，Z经过对角平均化得到的时间序列为z₁, z₂, …, z_N，则对角平均化的公式可表示为

(4)

其中M^*=min(M, K)，K^*=max(M, K).假设RC为Y_k，对应的元素为y₁^k, y₂^k, …, y_N^k，则任意两个重建向量之间的相关系数可表示为

(5)

式中，w_k=min(k, M, N－k).两个RC之间相关系数的绝对值越趋近于1，代表两个RC之间的相关性越强，极大可能代表同一种信号.

利用SSA对原始时间序列预测的本质就是利用SSA迭代插值(Kondrashov and Ghil, 2006; Kondrashov et al., 2010).将验潮时间序列中需要预测的数据位置用0填充，然后对补充后的时间序列做中心化处理.再利用SSA分解出的RC值代替0值，选择合适的迭代次数循环此过程，以达到最好的预测效果.具体步骤如下：

(1) 设置预测长度为m，用0值填充，则重建的时间序列长度为N+m，对其进行中心化处理；

(2) 对重建的时间序列进行SSA分解，预测位置的数据用第一个RC(记为RC₁)中的最后m个数据替代，循环此过程，直到两次预测数据残差的RMS小于0.001 m；

(3) 第一次循环结束后，通过添加第二个RC(RC₂)重建预测段数据，即预测部分数据由RC₁和RC₂通过线性叠加的方式得到，重复步骤(2)的过程，直至RC₁+RC₂序列收敛；

(4) 根据w-correlation和RC_i的贡献率，选择合适数量的RC重复上述过程，直至获得较为理想的预测结果.

2.2 ARMA模型预测

本文将验潮序列分离出主成分后的剩余序列作为一组随机序列，利用ARMA对其进行外推预测.首先利用ADF(Augmented Dickey-Fuller Test)检验方法检验剩余随机序列的平稳性，28个验潮站剩余随机序列的ADF检验结果P_value值均小于0.05，显著拒绝原假设(时间序列非平稳)，说明28个验潮站的剩余随机序列均通过平稳性检验.然后通过扩展自相关函数(Extended Auto Correlation Function，EACF)确定ARMA模型阶数p和q的值(Cryer and Chan, 2008)，通过最大自然估计确定模型各参数.ARMA(p, q)模型可表示为

(6)

式中为自回归(AR)系数，为滑动平均(MA)系数，(p, q)≥0且为整数.{ε_t}为白噪声序列，服从N(0, δ²)正态分布.

预测满足方差最小原则，由条件期望可得

(7)

(8)

(9)

当l>p且l>q时，则(9)式可表示为

(10)

其中表示由t开始进行预测长度为k步的预测.

将SSA主成分预测结果和ARMA外推结果合并即为最终验潮预测序列.以同样的方法处理GPS时间序列数据，从而获得地壳运动的预测值.

3 结果与分析 3.1 相对海平面变化预测

本文分别利用SSA和SSA+ARMA对日本沿岸的28个验潮站1994—2013年20年的验潮月数据进行分析，预测了日本沿岸2014—2018年相对海平面变化，并用相同时间段的实际测量数据对预测精度进行了检验.

利用SSA对时间序列数据进行预测，重点是确定窗口长度M和迭代阶次，一般且为周期的最小公倍数.日本沿岸验潮站的原始序列数据是1994—2013年20年的月数据，N=240，且日本海平变化受季风影响较大，年周期与半年周信号显著，所以窗口长度M=12.迭代阶次通过w-correlation分析结果和各RC的贡献率来确定，其中w-correlation分析能得到各RC之间的相关性，判断分离结果的好坏.本文研究了28个验潮站，由于篇幅限制，在此根据纬度从大到小选择了6个验潮站进行迭代阶次确定的说明.

如图 2所示，本文选择的6个验潮站前7阶RC_S分离的较好，对各重构成分进行快速傅里叶变换(Fast Fourier Transform, FFT)频谱分析后发现，RC₁和RC₂相关性较强，均代表周年信号.RC₃均代表非线性趋势项.各验潮虽然均位于日本沿岸，但各地所处的冷高压位置不同，所受季风、洋流也有所差异，所以各验潮站除周年信号、趋势主要的影响因素外，其他季节信号和年际信号对各验潮站的影响由于地理位置的不同而有所不同(Chen et al., 2000).验潮站Wakkanai、Hakodate I、Kashiwazaki、Aburatsu的RC₄+RC₅均代表半年周期信号，RC₆+RC₇均代表 2年的周期信号.而验潮站Kushimoto和Naha的RC₆+RC₇均代表半年周期信号，RC₄+RC₅分别代表 5年左右的周期信号和2年的周期信号.

由于日本沿岸28个验潮站所处地理位置不同，各周期信号对其海平面变化的影响也不同.通过SSA和w-correlation分析，根据各RC_S的贡献率(C_j>0.95)及相关性(ρ_{i, j}^w>0.9)，选择前7阶RC_S作为原始序列的主成分进行迭代预测.对于扣除前7阶RC_S成分的剩余随机序列可以建立ARMA(p, q)模型进行外推预测.根据EACF准则，ARMA(3, 4)比较适宜于验潮数据剩余随机序列的预测.ARMA模型外推剩余成分结果加上SSA主成分预测结果即为最终SSA+ARMA的预测结果.

如图 3所示，本文利用SSA对日本沿岸28个验潮站的2003—2018年验潮数据进行趋势提取并对其进行线性拟合，获得各验潮站相对海平面变化年速率.日本岛的相对海平面变化速率空间差异明显，以本州岛中部为分界线，分界线以北相对海平面上升速率较小，且部分区域的相对海平面呈下降趋势.分界线以南的相对海平面上升速率较大，验潮站Tajiri附近的相对海平面变化速率更是达到0.0083 m·a^-1.

为了评定各验潮站2014—2018年相对海平面变化的预测精度，将环日本沿岸28个验潮站2014—2018年的实际观测值作为真实值，采用均方根误差(Root Mean Square Error，RMSE)评定各验潮站5年的预测精度.

由图 4所示，对于大部分验潮站，SSA和SSA+ARMA的预测精度均随着预测时间的增加而减小，但在环日本沿岸相对海平面变化第5年的预测中，SSA+ARMA预测序列的均方根误差仍可达到0.0607 m的精度.且在为期5年的预测中，28个验潮站SSA+ARMA的预测精度均始终高于SSA预测精度，改善效果在第1年预测中尤为显著.

对于Oga验潮站，第4年的相对海平面变化的预测精度低于第5年的预测精度，因为在2017年Oga相对海平面高于往年，且当年的海平面的最低值也高于相邻年份的最低值.验潮站Abashiri第2年的预测精度较低，由于该站在2015年相对海平面变化幅度较大，近几年相对海平面高的极大值和极小值都出现在2015年.对于Mera、Naha、Hosojima预测精度不严格单调递减的验潮站，近海相对海平面变化均像Oga和Abashiri一样存在特殊状况，由于预测方法不能够对突发事件给出相应预测，所以该时段的预测精度较低，但整体预测趋势跟真值相差不大，预测结果还是具有参考意义.

3.2 地壳垂直变化预测

海平面的绝对变化包括海平面的相对变化和地壳的垂直运动.海平面的相对变化可以由验潮数据获得.地壳的垂直运动由与验潮站水尺零点联测的GPS测量得到.本文分别利用SSA和SSA+ARMA对28个联测GPS点2003—2013年的GPS月数据进行分析，并预测了各GPS站2014—2018年垂直变化.其中，利用SSA对日本沿岸28个联测GPS站2003—2018年GPS数据进行趋势提取并对其进行线性拟合，获得地壳垂直变化速率，结果如图 5所示.

如图 5所示，地壳垂直上升速率由西南向东北方向逐渐增加，在鄂霍次克板块和太平洋板块交界处地壳运动速率达到最大.且以本州岛中部，即鄂霍次克板块和阿穆尔板块交界处为分界线，分界线以北地壳运动以上升为主，分界线以南地壳运动以下降为主，与相对海平面变化(图 3)恰好相反.当相对海平面上升较快时，当地的地壳下降的也相对较快或上升较慢，两者的相关系数为-0.78.由此可以看出，在全球变暖的大趋势下，地面沉降(上升)已经成为日本岛沿海相对海平面上升(下降)的重要影响因素.

为了评定验潮站2014—2018年地壳垂直变化的预测精度，将同时段28个联测GPS的实际观测值作为真实值进行精度检验，统计结果如表 2所示.

由表 2可以看出，对于GPS时间序列数据，随着预测时间的增加，SSA和SSA+ARMA的预测精度都逐渐降低，但5年内的预测精度均在毫米级，且SSA+ARMA的预测精度均优于SSA预测.对于GPS时间序列数据，SSA+ARMA的预测相比于单纯用SSA预测，精度提高效果没有验潮数据预测明显.因为GPS趋势项、周期与半周期信号已经涵盖原始序列96%左右的信息，SSA分解重构后，剩余随机序列包含信息较少，所以经ARMA预测之后虽然精度有所提高，但是效果不明显.

3.3 绝对海平面变化预测

将28个验潮站2003—2018年的验潮数据及联测GPS数据线性相加，获得绝对海平面变化时间序列.利用SSA对绝对海平面变化时间序列进行趋势提取并对其进行线性拟合，获得日本沿岸2003—2018年海平面绝对变化速率，结果如图 6所示.

图中可见，2003—2018年日本沿岸的绝对海平面总体呈上升趋势，是相对海平面变化和地壳垂直变化共同作用的结果.由于日本沿岸东北部地壳垂直上升速率较快(图 5)，远远大于该区域相对海平面下降(上升)趋势，所以该区域的绝对海平面上升也比较较快.同时，验潮站Oga、Ogi附近海域相对海平面下降速率较快，大于该区域地壳垂直下降(上升)速率，该区域的绝对海平面变化呈下降趋势.

为了获得日本沿岸2014—2018年的近海绝对海平面变化，本文利用SSA和SSA+ARMA分别对1994—2013年的验潮数据和2003—2013年的GPS数据进行了分析和预测，把两者预测结果线性相加，作为最终的绝对海平面变化预测结果.把验潮数据和GPS数据的实际测量值线性相加最为绝对海平面变化的真值，并以此为参照，采用均方根误差评定SSA和SSA+ARMA的预测精度.

如图 7所示，随着预测时间的加长，SSA+ARMA和SSA的RMSE逐渐增加，且从2016年起，增加较为缓慢.在为期5年的近海绝对海平面预测中，SSA+ARMA的预测精度始终略高于SSA预测，但随着预测时间的加长，两者之间的差距逐渐缩小，可见ARMA短期预测效果比较明显.如表 3所示，虽然SSA和SSA+ARMA的预测精度随着时间的增加精度逐渐降低，但在5年的日本岛近海绝对海平面预测中，两者的预测精度都在厘米级，且2018年SSA+ARMA的预测精度为0.0683 m，远远大于卫星测高的近海精度.

为了更好的研究2014—2018年日本沿岸的绝对海平面变化变化速率及进一步检验SSA+ARMA的预测精度，本文分别对2014—2018年绝对海平面变化时间序列的测量值和预测值(时间序列数据经中心化处理)进行趋势提取并对其进行线性拟合，拟合结果如图 8所示.

图 8中2014—2018年日本沿岸绝对海平面变化的真实速率和预测的速率分别为0.0035 m·a^-1和0.0033 m·a^-1.预测速率与真实速率相差较小，说明SSA+ARMA趋势预测结果比较可靠.对于2014—2018年绝对海平面变化时间序列数据，无论是年变化速率还是绝对海平面高的峰值，预测值均略低于真实值.因为在全球变暖的大环境，加快了冰川融解速度及海水膨胀，海平面的变化速率也越来越快.且2014—2016年发生的超强厄尔尼诺事件，导致日本区域降水偏多，海平面上升的速率及幅度也比往年大.近5年环日本沿岸的近海绝对海平面变化年速率略高于全球海平面年变化速率，与Ho等(2005)的研究结果一致，这与环日本岛地壳整体上升和近海海水对温度变化较为敏感有关，但这样区域的上升速率结果对沿海区域的规划才具有针对性的现实意义.

基于以上模型预测精度，本文利用日本沿岸28个验潮站1999—2018年20年的验潮数据及2008—2018年11年的GPS数据(两者时间序列数据均经中心化处理)，通过SSA+ARMA预测模型，预测了环日本沿岸2019—2023年近海绝对海平面变化，预测结果如图 9所示.

用SSA+ARMA预测的2019—2023年日本近海绝对海平面高整体高于往年(图 9)，2019—2023年的平均海面高相对于2014—2018年平均海面高高0.0353 m.由于2019年海面高出现了小幅度的上升，2019—2023年绝对海平面呈略微下降趋势，而2003—2023年海平面的变化率为0.0039 m·a^-1，略高于2014—2018年的海面变化速率0.0035 m·a^-1.对于2019—2023年绝对海平面的预测，无论是19年的小幅度的上升，还是2003—2023年海平面的上升速率，均与Luo等(2011)利用主成分分析(Principal Component Analysis, PCA)预测的结果相吻合，预测结果较为理想.

4 结论

本文从日本沿岸选取了28个验潮站及联测GPS点，利用1994—2013年的验潮月数据和2003—2013年的GPS月数据，通过SSA和SSA+ARMA对日本沿岸2014—2018年的近海海平面变化进行了预测，并用同时段的实际测量值对预测结果进行检验，且在此基础上利用SSA+ARMA预测了日本沿岸2019—2023年的近海绝对海平面变化.实验结果显示，无论是对海平面预测还是地壳变化预测，SSA和SSA+ARMA的预测精度对随着时间的增加而降低，且SSA+ARMA的预测结果始终优于单纯SSA预测.

在2014—2018年间，SSA+ARMA对日本沿海相对海平面的预测精度为0.0357~0.0607 m，地壳垂直运动的预测精度为0.0049~0.0077 m，绝对海平面高的预测精度为0.0433~0.0683 m.且SSA+ARMA预测的2014—2018年的海平面年变化速率为0.0033 m·a^-1，与真实速率0.0035 m·a^-1相差较小.利用SSA+ARMA预测的2019—2023年的平均近海绝对海面高较往年(2014—2018)高0.0353 m，2003—2023年海平面的变化率为0.0039 m·a^-1，预测结果较为理想，精度比较可靠.

由于GPS观测年限限制，本文只用了11年的时间序列数据进行预测，如果能够获取更长的时间序列数据，预测结果可能将会更加精确.对于海平面变化预测，不同的预测方法数据适用性不同，优缺点也不同.因此，海平面变化预测应该集各种方法所长，根据原始序列数据特性选择不同的预测方法进行组合，构建完整优化的海平面预测体系，并在此基础上进行海平面变化的预测和研究.

致谢  感谢匿名审稿专家的意见和建议.感谢SONEL提供的GPS数据和PSMSL提供的验潮数据.