“21世纪海上丝绸之路”旨在通过海上互联互通推动中国与沿线国家的经贸合作，是实现沿线国家经贸互联互通的重要途径，其有两大走向：一是从中国沿海港口过南海，经过马六甲海峡到印度洋，延伸至欧洲；二是从中国沿海港口过南海，向南太平洋延伸。“21世纪海上丝绸之路”沿线国家众多，海上贸易联系密切，航线交通流量大，但沿线海峡众多，自然灾害频发导致主航线船舶交通流结构复杂，对船舶航行、海洋生态环境等造成严重挑战。

由于“21世纪海上丝绸之路”是实现沿线国家多元、可持续发展的关键所在，学者们已经对其进行了诸多方面的研究。贸易方面：张桂梅等^[1]基于拓展引力模型对海上丝绸之路西线国家做了国家贸易潜力研究，发现了对贸易量影响较大的正相关因素，并对各个国家进行了潜力型划分；王凤婷等^[2]利用随机前沿模型对海上丝绸之路沿线国家农产品出口贸易进行研究，分析了现状和影响因素，估算了贸易效率和出口潜力。港口航线方面：郑崇伟等^[3-4]搜集梳理了海上丝绸之路的重要航线、节点和港口的地理水文特征、气候特征、人文背景和商业价值等信息，集中展现了大量背景知识；牟乃夏等^[5]从海陆联运视角出发，构建港口区位优势评价模型，对海上丝绸之路沿线港口进行了区位优势度分析。网络方面：李振福等^[6]利用复杂网络分析方法，对海上丝绸之路沿线34个集装箱港口的各类中心性进行了测度，结果显示了香港港、深圳港、大连港、新加坡港和上海港在该集装箱海运网络中具有较高的地位。于安琪等^[7]分析了海上丝绸之路沿线港口的分布特点和网络特征，并利用班轮连通性指数测度了集装箱海运网络的连通性，指出了海上丝绸之路沿线国家连通性指数普遍较低且增长缓慢的现状。朱小林等^[8]针对21个海上丝绸之路沿线国家，建立区域整体贸易运输网络，量化了网络枢纽间规模经济效应和各个国家的物流绩效，通过遗传算法选择出更符合海上丝绸之路沿线区域联动发展的贸易海运网络。文化法律方面：孟夏韵^[9]探究了“一带一路”背景下中国文化在拉美的传播路径；杨泽伟^[10]研究了建设海上丝绸之路重要节点面临的政治风险和地区法律问题。安全方面：孙海泳^[11]研究了“一带一路”背景下，中非海上互通的安全风险，并且提出了防范风险的对策；宫晓婞等^[12]利用改进松弛变量模型对海上丝绸之路关键海峡、运河的安全效率进行了排序，指出海盗和船舶事故是影响安全效率的重要因素。海上丝绸之路海上交通流量预测对于提升船舶航行安全性，完善海上交通安全管理制度具有重要意义。但是，目前对海上丝绸之路交通流的研究相对较少。

交通流理论是运用物理和数学的定律来描述交通特性的一门边缘学科^[13]，广泛应用于交通工程的诸多研究领域。交通流研究最初集中于陆上交通，但随着全球卫星定位、船舶自动识别系统(Automatic Identification System，AIS)等技术的不断成熟，交通流理论逐渐被应用于海上船舶交通流特征统计及交通流量预测等方面。例如，甘浪雄等^[14]利用场方法统计了单位时间内某区域的船舶交通流密度；熊振南等^[15]利用点统计方法来判断船舶分布类型；金兴赋等^[16]、刘敬贤等^[17]针对特定水域，利用AIS数据分析船舶交通流特征；文哲等^[18]基于交通波理论，分析了枯水期长江中游的船舶流特征。近些年，研究人员在统计分析基础上，开始着手对交通流量预测进行研究，比如，惠子刚等^[19]建立自适应模糊推理系统对船舶交通流数据进行了仿真预测；钮浩东等^[20]运用优化的广义回归神经网络对东海大桥船舶流量进行了仿真预测；张晓雷等^[21]综合季节因素，建立灰色马尔可夫模型对船舶流进行预测；肖进丽等^[22]针对船舶交通流时间序列非平稳的特点，使用差分进化算法、BP神经网络等进行了较高精度的组合预测。上述研究均可以较好地预测船舶交通流，但是对交通流规律的解释性不够，不能反映当前船舶交通流会受哪些历史数据的影响。谢振东等^[23-24]指出交通流预测可以缓解交通拥堵的问题，Yip^[25]指出研究交通流可以更好地进行船舶交通安全管理，所以刻画交通流规律对通航管理具有一定的应用价值。

本文从交通流角度出发，基于2018年海上丝绸之路AIS数据，利用ARIMA模型(Auto-regressive Intergrated Moving Average Model，差分自回归移动平均模型)，分别对货船、油轮和货船−油轮3种情形下的船舶交通流进行拟合和预测，探究海上丝绸之路船舶交通流规律。

1 数据和模型描述 1.1 AIS数据描述

本文数据来源于《基于卫星组网的“21世纪海上丝绸之路”通道支点检测与保障研究》项目所提供的2018年海上丝绸之路AIS数据，包含以下2类信息：

1) 船舶静态信息。

船舶静态信息包括：MMSI码(Maritime Mobile Service Identify，水上移动通信业务标识码)、船名、船舶类型、IMO编号(International Maritime Organization，国际海事组织，内贸船无该编号)、船长、船宽、船旗国7个字段数据。

2) 船舶动态信息。

船舶动态信息包括：MMSI码、时间、经度、纬度、航向、航速、航行状态7个字段数据。

上述信息中，MMSI码用来识别属于同一船舶的信息，船舶类型用来区分不同类型的船舶，经纬度用来判断船舶的活动位置，时间、航速和航向可以反映船舶活动规律并用来研究交通流。本文根据研究需要，将提取船舶的MMSI码、船舶类型和船舶动态等信息。

1.2 AIS数据筛选

2018年全年AIS数据量巨大，并且动态信息记录存在诸多问题，例如，动态信息时间间隔不规范，船舶无有效动态记录(即航速始终为0)，船舶航速记录异常等。因此，本文采取以下3个筛选原则：

(1) 筛除无有效动态信息的船舶，即航速始终为0的船舶；

(2) 筛除船舶异常航速信息，即航速明显超出该类型船舶的正常速度范围的船舶记录；

(3) 由于动态数据时间间隔并不规范，所以本文只抽取每小时内的第一条船舶动态数据记录。

此外，船舶静态信息也存在一些问题，例如，船长、船宽和船舶类型等静态信息记录存在错误或缺失。借助船讯网( http://www.shipxy.com/)和船队在线网( http://www.hifleet.com/)补全更正船舶静态信息。

经过以上处理，得到31 554条船舶的有效动态信息，其中包含货船、油轮和捕捞船等23个类型的船舶，具体船舶数量见图1，本文重点研究货船和油轮。

1.3 时间序列预测模型构建

船舶交通流时间序列是指在一定时间范围内，以固定时间间隔排列的船舶数量序列。虽然构成该序列的任一序列值具有不确定性，但序列值之间存在一定关联，因此可通过构建数学模型对后续序列值进行预测。常用的时间序列预测模型包括AR模型(Autoregressive Model，自回归模型)、MA模型(Moving Average Model，移动平均模型)、ARMA模型(Autoregressive Moving Average Model，自回归移动平均模型)和ARIMA模型，其中ARMA模型是AR模型和MA模型的组合，ARIMA模型是在ARMA模型的基础上增加了差分处理的改进。

1.3.1 模型简介

ARIMA模型由美国统计学家Box和英国统计学家Jenkins提出，是一类预测精度较高的时间序列模型，广泛应用于各类时间序列的预测^[26-30]，其特点是利用历史数据就可以预测未来响应的时间序列。该模型将具有明显趋势的非平稳数据，经过差分处理使其转化为平稳序列，建立合适的模型，利用历史数据预测下一阶段的数据。ARIMA模型由3部分组成，包括自回归过程AR、差分过程I和移动平均过程MA。

自回归过程的公式： ${y_t} = \mu + \displaystyle \sum\limits_{i = 1}^p {{\gamma _i}{y_{t - i}}} + {\varepsilon _t}$ ，其中 ${y_t}$ 是当前值， $\;\mu$ 是常数项， $p$ 表示阶数(表示当前值与其前 $p$ 个序列值有关)， ${\gamma _i}$ 是自相关系数， ${\varepsilon _t}$ 是误差。自回归模型利用自身的数据进行预测，所以必须具有自相关性，一般要求自相关系数不小于0.2。

移动平均过程关注的是自回归模型中误差项的累加，AR与MA结合得到自回归移动平均模型ARMA，其表达式为 ${y_t} = \mu + \displaystyle \sum\limits_{i = 1}^p {{\gamma _i}{y_{t - i}}} + {\varepsilon _t} + \displaystyle \sum\limits_{i = 1}^q {{\theta _i}{\varepsilon _{t - i}}}$ ，其中 $q$ 是移动平均的阶数， ${\theta _i}$ 是移动平均参数。ARMA模型中 $p$ ， $q$ 的值需要由自相关函数和偏相关函数来确定。

差分过程I表示对原始时间序列进行差分处理，一般用 $d$ 来表示 $d$ 阶差分。若原始数据已经是平稳序列，那么 $d = 0$ ，一般情况下 $d = 1$ 。

综上，差分自回归移动平均模型一般用 ${\rm{ARIMA}} (p,d,q)$ 来表示，其中 $p$ 是自回归模型阶数， $d$ 是差分的阶数， $q$ 是移动平均的阶数。

1.3.2 建模步骤

本文利用Python实现ARIMA模型建模过程，调用了Pandas，Numpy和StatsModels库，选取从2018年1月1日起至2018年12月24日止的货船和油轮动态AIS数据，记录每隔7天的2类船舶活动数量，每类船舶有52条数据记录，设置85%的训练集构建模型，15%的测试集来验证模型的准确性。ARIMA模型建模步骤如图2所示。

2 结果与讨论

本节以货船为例，实现ARIMA预测模型。

2.1 差分平稳化

货船交通流原始时间序列数据如图3所示。

从图3观察可知，该时间序列是非平稳的，首先对该原始时间序列进行一阶差分处理，得到一阶差分时间序列，如图4所示。

由图4可知，一阶差分时间序列比原始时间序列更趋于平稳。分别对其进行平稳性检验，检验序列中是否存在单位根，若存在单位根就是非平稳时间序列。对比结果见表1。

表1中ADF Test值(T值)是对时间序列的检验结果，显著性水平值表示不同显著性水平下拒绝原假设的统计值，ADF检验的原假设是存在单位根，若T值小于1%、5%、10%这3个显著性水平值，即说明非常好地拒绝该假设。Prob(P值)表示T值的概率值，P值要求小于给定的显著性水平，若小于1%就可以极显著地拒绝原假设，表明数据平稳。由表1可知，一阶差分时间序列不存在单位根，为平稳序列。

2.2 模型选择

取前85%的原始时间序列值作为训练集，进行一阶差分处理后得到数据时间序列，分别计算其自相关函数(Autocorrelation Function，ACF)和偏相关函数(Partial Autocorrelation Function，PACF)。

ACF计算公式为 $\rho = \displaystyle\sum\limits_{i = 1}^{n - h} {\dfrac{{\left( {{x_i} - \overline \mu } \right)\left( {{x_{i + h}} - \overline \mu } \right)}}{{\displaystyle\sum\nolimits_{i = 1}^n {{{\left( {{x_i} - \overline \mu } \right)}^2}} }}}$ ，其中 ${x_i}$ 表示第 $i$ 个序列值， ${x_{i + h}}$ 表示滞后阶数为 $h$ 的序列值， $\overline \mu $ 表示 $n$ 个序列值的平均值。ACF反映的是一组数据与其滞后数据间的相关性，即 ${x_i}$ 与 ${x_{i + 1}}, \cdots ,{x_{i + h}}$ 之间关系。而PACF则反映剔除 ${x_{i + 1}}, \cdots ,{x_{i + h - 1}}$ 干扰后 ${x_i}$ 对 ${x_{i + h}}$ 的影响。时间序列数据对应的ACF和PACF图，如图5所示。

图5中，横坐标表示滞后阶数，阴影部分表示置信区间，当滞后阶数为2时，ACF和PACF超过了置信区间。 $p$ 和 $q$ 的选择原则见表2。

从图5中ACF和PACF图的拖尾性来看， $p = 1$ ， $q = 2$ ，可以考虑使用ARMA(1,0)模型，ARMA(0,2)模型和ARMA(1,2)模型。借助计算最小化信息量AIC(Akaike Information Criterion)和贝叶斯信息量BIC(Bayesian Information Criterion)来选择合适的模型，计算公式分别为 ${\rm{AIC}} = 2k - 2\ln L$ ， ${\rm{BIC}} = k\ln n - 2\ln L$ ，其中 $k$ 是模型参数个数， $n$ 是样本数量， $L$ 是似然函数。利用AIC和BIC可以在有限的阶数范围内寻找相对最优拟合模型，有效弥补根据ACF图和PACF图定阶的不足。计算 $p,q \leqslant 4$ 所有组合下AIC、BIC的值，取其信息量值最小的阶数组合，计算结果见表3和表4。

由表3和表4可知，AIC值最小的是ARMA(1,2)模型，BIC值最小的是ARMA(1,0)模型，分别求其残差对应的ACF和PACF，见图6。

图6中显示ARMA(1,2)模型下残差的ACF和PACF均在置信区间内，不存在相关关系。说明在该模型下，数据序列的有效信息已经被全部提取，拟合效果显更优，所以选择ARMA(1,2)模型。

2.3 模型检验

对选定的ARMA(1,2)模型进行检验，分为2个步骤：(1) 选用Durbin-Watson检验(D-W检验)来检查模型残差的自相关性，D-W值接于0或4，说明存在自相关性，而D-W值接近于2时，说明不存在相关性；(2) 选用Ljung-Box Q检验(LBQ检验)来判断残差是否为高斯白噪声。在ARIMA模型中残差被假定为高斯白噪声序列，如果检验概率小于给定显著性水平(5%或10%)，那么残差就不是高斯白噪声，说明该模型并不是一个合适的模型。

ARMA(1,2)模型下，残差的D-W值为1.8261，接近于2，说明残差不存在自相关性。残差的LBQ检验结果见表5，其中各阶滞后的检验概率均大于5%，说明ARMA(1,2)模型的残差是高斯白噪声序列，即该模型是一个合适的模型。

2.4 模型预测

利用ARMA(1,2)模型进行预测，模型表达式为

$\Delta {y_t} = 9.63 + 0.45\Delta {y_{t - 1}} + {\varepsilon _t} + 0.43{\varepsilon _{t - 1}} + 0.56{\varepsilon _{t - 2}}$

将表达式转换后，得到

$ {y_t} = 9.63 + 1.45{y_{t - 1}} - 0.45{y_{t - 2}} + {\varepsilon _t} + 0.43{\varepsilon _{t - 1}} + 0.56{\varepsilon _{t - 2}} $

由上述表达式可以看出，当前船舶交通流会受到之前2个时间周期内的船舶交通流的影响，预测结果、真实值及误差见表6。

常用的回归预测评价方法有均方误差 ${\rm{MSE}}$ 、均方根误差 ${\rm{RMSE}}$ 、平均绝对误差 ${\rm{MAE}}$ 和决定系数R²。R²不存在量纲问题，取值范围为[0,1]，R²越接近1说明模型拟合效果越好。

$ {R^{\rm{2}}} = 1{\rm{ - }}{{\Bigg( {\displaystyle\sum\limits_{i = 1}^m {{{({y_i} - y_i')}^2}} } \Bigg)} \Big/ {\Bigg( {\displaystyle\sum\limits_{i = 1}^m {{{(\overline y - y_i')}^2}} } \Bigg)}} $

其中， $m$ 表示真实值和预测值的数量， ${y_i}$ 表示真实值， ${y'}$ 表示预测值， $\overline y $ 表示真实值的平均值。由于MSE、RMSE和MAE均有量纲存在且没有一定取值范围，可读性不强。

根据表6结果，计算ARMA(1,2)模型的决定系数 ${R^{\rm{2}}} = 0.8{\rm{09\;19}}$ ，表示预测效果比较好。由图7可知，预测曲线趋势与实际值趋势基本一致。

利用ARIMA模型对油轮和货船−油轮这两类情形的交通流进行预测，建模过程与研究货船时一致，各情形的模型选择和R²值，见表7。

从表7可知，考虑货船、油轮和货船−油轮这3种情形时，最优模型均为ARIMA(1,1,2)，表明当前交通流会受到近两周内的交通流的影响。对比不同情形下的R²值，可以发现ARIMA模型对货船和油轮交通流预测效果明显优于对货船−油轮交通流预测效果，表明ARIMA模型更适合用来研究单一船型的交通流。

3 结论

本文基于2018年“21世纪海上丝绸之路”AIS数据，针对货船、油轮、货船−油轮3种情形，构建了ARIMA模型拟合船舶交通流并做出预测。结果显示，(1) ARIMA是一个合适的模型，对于不同的情形都可给出较好的拟合预测；(2) 由ARIMA(1,1,2)模型得出的拟合公式可知，当前船舶交通流会受之前2个时间周期内的交通流影响，这为下一周期的船舶水路规划和通航安全管理提供充裕的时间保障；(3) 由本文所求得的R²值可知，ARIMA(1,1,2)模型对单一船型的交通流预测效果明显优于货船−油轮组合船舶。综上表明ARIMA预测模型对“21世纪海上丝绸之路”海上交通流研究具有一定实用价值。建立“21世纪海上丝绸之路”是中国加强与沿线国家贸易联系和打造区域经济一体化的重要战略。沿线的中国港口可以基于船舶到港口数据或AIS数据，利用ARIMA模型对船舶交通流进行预测，确保港口航线通航安全，针对特定船型，通过交通流预测可以预估贸易量，了解其他国家货物需求，这对我国港口发展建设及港口腹地经济研究具有重要意义。

本文不足之处：(1) 项目只包含了2018年的AIS数据，数据量偏少，因此导致预测的精度下降；(2) 影响海上交通流的因素众多，本文利用ARIMA模型研究海上丝绸之路交通流，仅使用了AIS的历史数据作为研究的时间序列，不能反映不确定性因素。