全球气候变化、北极海冰大量融化，使得北极水域船舶商业化运营成为了可能。2016年，我国“永盛轮”、“夏之远”、“天禧轮”等多艘船舶穿越北极东北航道顺利抵达了欧洲港口。然而北极水域作为一个船舶航行高风险、环境脆弱敏感的冰区水域，船舶航行面临复杂海冰、地磁干扰、大风、低能见度等复杂环境条件的考验。为了提升极地水域船舶航行安全水平，开展风险因素识别研究十分必要。

北极水域船舶航行安全研究受到了俄罗斯^[1]、加拿大^[2-3]、挪威^[4]等北极周边国家学者的广泛关注。加拿大交通部门构建了北极冰域航运系统^[2]；KHAN等提出了北极航运系统风险评价框架^[3]；KUM等分析了北极水域船舶碰撞、搁浅事故致因^[4]；MARKEN等分析了东北航道船舶延误的风险^[5]；极地规则^[6]也对不同海冰下各种冰级船舶的航行速度进行了分析。吕宝刚等分析了北极东北航道的船舶通航环境特征^[7-8]；FU等建立了北极水域船舶冰困风险的多因素耦合模型^[9]。

这些研究主要以通航环境和事故风险评价研究为主，运用北极水域船舶通航环境影响因素，构建了碰撞、搁浅等北极水域典型事故的风险评价模型，但对于风险因素的辨识缺乏系统性研究，冰流、磁暴等部分环境风险因素没有在模型中予以考虑。鉴于此，本文针对北极水域船舶航行所涉及的环境因素进行系统性的分析和识别。

1 环境风险因素的层次模型

北极水域是一个涉及复杂气象、海况和冰况的动态环境。本研究通过文献分析和专家调研，识别船舶在北极水域航行面临风、能见度、地磁、气温、磁暴、海冰密集度、冰流等环境风险影响因素^{[6, 9-11]}。针对北极水域船舶航行可能面临的这些环境风险因素，运用层次分析法建立北极水域船舶航行环境风险因素的层次模型，如图 1所示。

该模型以识别北极水域船舶航行环境风险因素为目标，将其划分为气象环境和航道条件2个子系统作为指标层，气象环境子系统被进一步划分为风、能见度、地磁、低温、磁暴5个底层分指标，航道条件子系统被划分为海况、冰况2个分指标层。其中，海况分指标又被进一步划分为海水温度和浪2个底层分指标，冰况分指标层又被进一步划分为海冰密集度、海冰厚度和冰流3个底层分指标。

2 风险因素识别方法

针对北极水域船舶航行环境风险因素中面临的不确定性问题，本研究在图 1所构建的北极水域船舶航行环境风险因素的层次分析模型基础上，提出了一种基于蒙特卡洛仿真的模糊层次分析法(Monte Carlo-based fuzzy-AHP, MC-FAHP)。通过北极水域船舶航行环境风险影响因素层次模型的因素权重计算，以识别影响北极水域船舶航行安全的重要环境风险因素。具体包含以下4个研究步骤：

1) 模糊判断矩阵构建。建立专家判断术语与三角模糊数之间的映射关系，搜集专家对于层次模型各层级指标的判断意见，针对每个专家意见分别构建成对比较的模糊判断矩阵。

2) 权向量的概率分布计算。根据模糊层次分析法计算规则，计算模糊判断矩阵的权向量，结合可能性分析方法，将模糊判断矩阵中涉及的权向量从模糊隶属度函数转化为概率分布形式。

3) 专家权重计算。针对参与问卷调研专家的职位、工作经验、教育背景等方面因素，构建专家权重评价指标体系，评价问卷调研专家的权重系数，并根据参与问卷调研专家的权重系数，运用加权平均方法，对权向量的概率分布函数进行数据融合。

4) 合成权重计算。根据专家意见融合后各指标的权向量，计算底层指标对目标层(北极水域船舶航行环境风险因素识别)合成权重的概率分布函数，并根据概率分布函数众数的大小对其进行排序，识别关键风险因素(指标)。

2.1 模糊判断矩阵的构建

根据层次分析法常用的九级标度法^[12]，结合模糊层次分析法三角模糊集合构建方式^[13]，定义成对比较的九标度模糊判断术语来搜集专家意见，如表 1所示。

根据表 1所述的模糊判断术语，搜集每个专家的判断意见，对各层级指标构建模糊判断矩阵。假设层次模型中某层级具有j个因素，用 ${{\tilde a}_{ij}}$ 表示第i个因素对于第j个因素的重要性比较结果，该层级的判断矩阵 ${\mathit{\boldsymbol{\tilde A}}}$ 可以表示为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde A}} = {{\left( {{{\tilde a}_{ij}}} \right)}_{n \times m}} = \left[ {\begin{array}{*{20}{c}} {\left( {1,1,1} \right)}&{\left( {{l_{12}},{m_{12}},{u_{12}}} \right)}& \cdots &{\left( {{l_{1n}},{m_{1n}},{u_{1n}}} \right)}\\ {\left( {1/{u_{12}},1/{m_{12}},1/{l_{12}}} \right)}&{\left( {1,1,1} \right)}& \cdots &{\left( {{l_{2n}},{m_{2n}},{u_{2n}}} \right)}\\ \vdots&\vdots&\vdots&\vdots \\ {\left( {1/{u_{1n}},1/{m_{1n}},1/{l_{1n}}} \right)}&{\left( {1/{u_{2n}},1/{m_{2n}},1/{l_{2n}}} \right)}& \cdots &{\left( {1,1,1} \right)} \end{array}} \right],}\\ {i,j = 1,2, \cdots ,n;i < j} \end{array} $

(1)

2.2 权向量的概率分布计算

根据模糊层次分析法^[13]和Wang等^[14]提出的一个修正公式，计算模糊判断矩阵 ${\mathit{\boldsymbol{\tilde A}}}$ 列向量之和：

$ \begin{array}{*{20}{c}} {R{\mathit{\boldsymbol{S}}_i} = \sum\limits_{j = 1}^n {{{\tilde a}_{ij}}} = \left( {\sum\limits_{j = 1}^n {{l_{ij}}} ,\sum\limits_{j = 1}^n {{m_{ij}}} ,\sum\limits_{j = 1}^n {{u_{ij}}} } \right),}\\ {i = 1,2, \cdots ,n} \end{array} $

(2)

然后对模糊判断矩阵 ${\mathit{\boldsymbol{\tilde A}}}$ 进行归一化计算，得到各层级权向量的三角模糊数：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde S}}}_i} = \frac{{R{\mathit{\boldsymbol{S}}_i}}}{{\sum\nolimits_{j = 1}^n {R{\mathit{\boldsymbol{S}}_j}} }} = }\\ {\left[ {\frac{{\sum\limits_{j = 1}^n {{l_{ij}}} }}{{\sum\limits_{j = 1}^n {{l_{ij}}} + \sum\limits_{k = 1,k \ne i}^n {\sum\limits_{j = 1}^n {{u_{kj}}} } }},\frac{{\sum\limits_{j = 1}^n {{m_{ij}}} }}{{\sum\limits_{k = 1}^n {\sum\limits_{j = 1}^n {{m_{kj}}} } }},} \right.}\\ {\left. {\frac{{\sum\limits_{j = 1}^n {{u_{ij}}} }}{{\sum\limits_{j = 1}^n {{u_{ij}}} + \sum\limits_{k = 1,k \ne i}^n {\sum\limits_{j = 1}^n {{l_{kj}}} } }}} \right],}\\ {i = 1,2, \cdots ,n} \end{array} $

(3)

式中： ${{\mathit{\boldsymbol{\tilde S}}}_i}$ 表示判断矩阵 ${\mathit{\boldsymbol{\tilde A}}}$ 第i个指标权重的三角模糊数， ${\left( {{{\mathit{\boldsymbol{\tilde S}}}_1},{{\mathit{\boldsymbol{\tilde S}}}_2}, \cdots ,{{\mathit{\boldsymbol{\tilde S}}}_n}} \right)^{\rm{T}}}$ 组成模糊判断矩阵 ${\mathit{\boldsymbol{\tilde A}}}$ 的特征向量。

对模糊集合 ${{\mathit{\boldsymbol{\tilde S}}}_i}$ 取α截集，获得置信度1-α下的模糊截集 $\tilde S_i^\alpha $ ：

$ \begin{array}{l} \mathit{\boldsymbol{\tilde S}}_i^\alpha = \left( {L_i^\alpha ,U_i^\alpha } \right) = \left( {\alpha \left( {{M_i} - {L_i}} \right) + {L_i}} \right.,\\ \;\;\;\;\;\;\;\;\left. {\alpha \left( {{U_i} - {M_i}} \right) + {M_i}} \right),i = 1,2, \cdots ,n \end{array} $

(4)

当α=0时， $\mathit{\boldsymbol{\tilde S}}_i^\alpha $ 表示100%置信区间；当α=5%时， $\mathit{\boldsymbol{\tilde S}}_i^\alpha $ 表示95%置信区间；当α=95%时， $\mathit{\boldsymbol{\tilde S}}_i^\alpha $ 表示5%置信区间。

在蒙特卡洛仿真的过程中，通过生成符合一定概率分布的随机变量，用统计方法估计模型的数字特征，从而得到实际问题的数值解。由于专家判断的输入数据存在一定的随机不确定性，拟采用三角形分布函数来对这种不确定性进行数值仿真。三角形分布是以下限为a，众数为b，上限为c的连续概率分布：

$ f\left( {x\left| {a,b,c} \right.} \right) = \left\{ \begin{array}{l} \frac{{2\left( {x - a} \right)}}{{\left( {b - a} \right)\left( {c - a} \right)}},a \le x \le b\\ \frac{{2\left( {b - x} \right)}}{{\left( {b - a} \right)\left( {b - c} \right)}},b < x \le c \end{array} \right. $

(5)

按照下式将模糊判断矩阵指标权重转化为三角形分布函数：

$ \begin{array}{*{20}{c}} {\left\{ \begin{array}{l} {a_i} = \frac{{\sum\limits_{j = 1}^n {{l_{ij}}} }}{{\sum\limits_{j = 1}^n {{l_{ij}}} + \sum\limits_{k = 1,k \ne i}^n {\sum\limits_{j = 1}^n {{u_{kj}}} } }}\\ {b_i} = \frac{{\sum\limits_{j = 1}^n {{m_{ij}}} }}{{\sum\limits_{k = 1}^n {\sum\limits_{j = 1}^n {{m_{kj}}} } }}\\ {c_i} = \frac{{\sum\limits_{j = 1}^n {{u_{ij}}} }}{{\sum\limits_{j = 1}^n {{u_{ij}}} + \sum\limits_{k = 1,k \ne i}^n {\sum\limits_{j = 1}^n {{l_{kj}}} } }} \end{array} \right.}\\ {i = 1,2, \cdots ,n} \end{array} $

(6)

$ \begin{array}{*{20}{c}} {{f_i}\left( {{x_i}\left| {{a_i},{b_i},{c_i}} \right.} \right) = }\\ {f\left( {{x_i}\left| {\frac{{\sum\limits_{j = 1}^n {{l_{ij}}} }}{{\sum\limits_{j = 1}^n {{l_{ij}}} + \sum\limits_{k = 1,k \ne i}^n {\sum\limits_{j = 1}^n {{u_{kj}}} } }},\frac{{\sum\limits_{j = 1}^n {{m_{ij}}} }}{{\sum\limits_{k = 1}^n {\sum\limits_{j = 1}^n {{m_{kj}}} } }}} \right.} \right.}\\ {\left. {\frac{{\sum\limits_{j = 1}^n {{u_{ij}}} }}{{\sum\limits_{j = 1}^n {{u_{ij}}} + \sum\limits_{k = 1,k \ne i}^n {\sum\limits_{j = 1}^n {{l_{kj}}} } }}} \right)} \end{array} $

(7)

然后通过随机生成1 000组样本数据^[15]来仿真三角形分布函数的随机性，来估算各层级权向量的众数和置信区间。

2.3 专家权重计算

拟根据调研人员职位(e₁)、航运系统工作经验(e₂)、北极航运工作经验(e₃)、教育背景(e₄)等方面因素构建人员工作经验评价指标体系，如表 2所示，对调研专家的资历进行综合评价，从而获得各位专家在模型中的权重。

假设有m位专家参与调研工作，各位专家综合得分S^(k)、权重λ^(k)按照下式进行计算：

$ {S^{\left( k \right)}} = S_{{e_1}}^{\left( k \right)} + S_{{e_2}}^{\left( k \right)} + S_{{e_3}}^{\left( k \right)} + S_{{e_4}}^{\left( k \right)},q = 1,2, \cdots ,m $

(8)

$ {\lambda ^{\left( k \right)}} = \frac{{{S^{\left( k \right)}}}}{{\sum\limits_{k = 1}^m {{S^{\left( k \right)}}} }},k = 1,2, \cdots ,m,\sum\limits_{k = 1}^m {{\lambda ^{\left( k \right)}}} = 1 $

(9)

2.4 合成权重计算

根据各位专家的权重，按照加权平均的方法对权向量的三角形分布函数进行数据融合：

$ \begin{array}{*{20}{c}} {\bar f\left( {{x_i}\left| {{a_i},{b_i},{c_i}} \right.} \right) = \sum\limits_{e = 1}^m {{\lambda ^{\left( e \right)}}{f^{\left( e \right)}}\left( {{x_i}\left| {{a_i},{b_i},{c_i}} \right.} \right)} ,}\\ {i = 1,2, \cdots ,n} \end{array} $

(10)

将目标层下的指标层定义为第1层，各分指标层分别表示为第2~(n-1)层，对各层级的权重进行合成，计算各层级相对于目标层合成权重的分布函数：

$ \begin{array}{*{20}{c}} {\bar W_i^{\left( k \right)} = \prod\nolimits_{k = 1}^{n - 1} {\bar w_i^{\left( k \right)}} = \prod\nolimits_{k = 1}^{n - 1} {\bar f_i^{\left( k \right)}\left( {{x_i}\left| {{a_i},{b_i},{c_i}} \right.} \right)} ,}\\ {i = 1,2, \cdots ,n} \end{array} $

(11)

式中：w_i^(k)表示第k层指标的权重，W_i^(k)表示第k层指标对于目标层的合成权重，f_i^(k)(x_i|a_i, b_i, c_i)表示第k层指标合成权向量的分布函数。

3 环境风险因素的识别研究 3.1 模糊判断矩阵构建

根据2.1节所述的模糊判断术语，通过问卷调研搜集到6名专家对于北极水域航行船舶环境风险因素的判断意见。表 3为各位专家对于气象环境分指标(C₁₁~C₁₅)的判断意见。

将每位专家对于各层级指标的两两比较判断意见组成各层级模糊判断矩阵。例如A专家对于环境风险因素的判断矩阵可以表示为

$ \mathit{\boldsymbol{\bar A}}_{{C_1}}^{\left( A \right)} = {\left[ {\begin{array}{*{20}{c}} {\left( {1,1,1} \right)}&{\left( {1/6,1/5,1/4} \right)}&{\left( {2,3,4} \right)}&{\left( {1,2,3} \right)}&{\left( {4,5,6} \right)}\\ {\left( {4,5,6} \right)}&{\left( {1,1,1} \right)}&{\left( {6,7,8} \right)}&{\left( {4,5,6} \right)}&{\left( {8,9,9} \right)}\\ {\left( {1/4,1/3,1/2} \right)}&{\left( {1/8,1/7,1/6} \right)}&{\left( {1,1,1} \right)}&{\left( {1/3,1/2,1} \right)}&{\left( {1,2,3} \right)}\\ {\left( {1/3,1/2,1} \right)}&{\left( {1/6,1/5,1/4} \right)}&{\left( {1,2,3} \right)}&{\left( {1,1,1} \right)}&{\left( {2,3,4} \right)}\\ {\left( {1/6,1/5,1/4} \right)}&{\left( {1/9,1/9,1/8} \right)}&{\left( {1/3,1/2,1} \right)}&{\left( {1/4,1/3,1/2} \right)}&{\left( {1,1,1} \right)} \end{array}} \right]^{\left( A \right)}} $

3.2 权向量的概率分布计算

按照式(3)计算模糊判断矩阵的权重，可得到 $\mathit{\boldsymbol{\tilde A}}_{{C_1}}^{\left( A \right)}$ 矩阵的权向量：

$ \mathit{\boldsymbol{\tilde S}}_{{C_1}}^{\left( A \right)} = {\left[ {\begin{array}{*{20}{c}} {0.1459}&{0.2195}&{0.3076}\\ {0.4179}&{0.5292}&{0.6351}\\ {0.0458}&{0.7790}&{0.1312}\\ {0.0785}&{0.1313}&{0.2056}\\ {0.0305}&{0.0420}&{0.0697} \end{array}} \right]^{{{\rm{T}}^{\left( A \right)}}}} $

按照式(6)、(7)将三角模糊向量转化为三角形分布函数：

$ \begin{array}{*{20}{c}} {f_{{C_{11}}}^{\left( A \right)}\left( {{x_{{C_{11}}}}\left| {{a_{{C_{11}}}},{b_{{C_{11}}}},{c_{{C_{11}}}}} \right.} \right) = }\\ {f\left( {{x_{{C_{11}}}}\left| {0.1459,0.2195,0.3076} \right.} \right)} \end{array} $

3.3 专家权重计算

按照式(8)、(9)计算得到6位专家(A~F)的权重λ_k分别为(0.238 8, 0.090 9, 0.155 8, 0.194 8, 0.142 9, 0.181 8)^T。按照式(10)对6位专家数据进行融合，得到环境风险因素C₁₁指标的概率密度函数和累积分布函数，如图 2所示。

3.4 合成权重计算

根据公式(11)计算各层级相对于目标层合成权重的分布函数，表 4为计算得出的底层指标相对于目标层合成权重，分别用95%置信程度的区间[f_Ci, f_Ci]、众数和百分比这3个统计量进行表示。

由表 4可见，风、能见度和海冰密集度是影响北极水域船舶航行安全的首要因素，它们众数的权重(在95%置信度情况下)占总权重的13%以上；气温、海冰厚度对北极船舶航行安全的影响次之，它们众数的权重占总权重的9%~10%；冰流、地磁和磁暴是影响北极船舶航行安全的一般因素，它们众数的权重占总权重的6%~7%；海水温度和浪对于北极水域船舶航行安全的影响较小，它们众数的权重约在5%以下。

3.5 分析与验证

从数据和模型这2个方面，对本研究北极水域船舶航行环境风险识别研究的有效性进行分析与验证。

3.5.1 专家判断数据

通过将模糊矩阵中的三角模糊数进行去模糊化计算转换成常规矩阵，可以应用层次分析法的一致性检验方法，对每位专家数据的进行检验。

下面以 $\mathit{\boldsymbol{\tilde A}}_{{C_1}}^{\left( A \right)}$ 为例对计算过程进行介绍。取 $\mathit{\boldsymbol{\tilde A}}_{{C_1}}^{\left( A \right)}$ 矩阵中每个三角模糊数的中间值构成非模糊的AHP判断矩阵A_C1^(A)，并对其进行一致性检验^[13]，计算结果如下:

$ \mathit{\boldsymbol{A}}_{{C_1}}^{\left( A \right)} = {\left[ {\begin{array}{*{20}{c}} 1&{1/5}&3&2&5\\ 5&1&7&5&9\\ {1/3}&{1/7}&1&{1/2}&2\\ {1/2}&{1/5}&2&1&3\\ {1/5}&{1/9}&{1/2}&{1/3}&1 \end{array}} \right]^{\left( A \right)}}, $

$ \begin{array}{*{20}{c}} {w_{{C_1}}^{\left( A \right)} = \left( {0.1923,0.5700,0.0729,0.1207,0.0441} \right),}\\ {{\lambda _{\max }} = 5.1191,{\rm{CI}} = 0.0298,{\rm{CR}} = 0.0266} \end{array} $

根据表 2所述的一致性判断标准，A_C1^(A)矩阵的CI值为0.029 8，小于0.1，认为该判断矩阵的一致性可以接受。

类似地，计算得出每位专家判断矩阵均能通过一致性检验。

3.5.2 模型效度验证

1) 研究方法对比验证

根据CHANG^[13]和WANG^[14]等所述的模糊层次分析法，对6名专家的判断意见进行计算，得到模糊层次分析法层次模型底层指标相对于目标层合成权重，如表 5所示。模糊层次分析法的计算结果用权向量( ${{\mathit{\boldsymbol{\tilde S}}}_i}$ )、模糊区间间隔和去模糊化数值这3个指标表示，以便跟基于蒙特卡洛仿真的模糊层次分析法(MC-FAHP)的计算结果进行比较。其中，模糊层次分析法的去模糊化，按照以下公式进行计算：

$ \begin{array}{*{20}{c}} {{V_{{\rm{defuzzification}}}} = }\\ {\frac{1}{3}\left( {\frac{{\sum\nolimits_{j = 1}^n {{l_{ij}}} }}{{\sum\nolimits_{j = 1}^n {{l_{ij}}} + \sum\nolimits_{k = 1,k \ne i}^n {\sum\nolimits_{j = 1}^n {{u_{kj}}} } }} + \frac{{\sum\nolimits_{j = 1}^n {{m_{ij}}} }}{{\sum\nolimits_{k = 1}^n {\sum\nolimits_{j = 1}^n {{m_{kj}}} } }} + } \right.}\\ {\left. {\frac{{\sum\nolimits_{j = 1}^n {{u_{ij}}} }}{{\sum\nolimits_{j = 1}^n {{u_{ij}}} + \sum\nolimits_{k = 1,k \ne i}^n {\sum\nolimits_{j = 1}^n {{l_{kj}}} } }}} \right)} \end{array} $

从模糊层次分析法方法的去模糊化数值来看，模糊层次分析法方法的计算结果与MC-FAHP的结果呈现出高度的一致性，验证了基于蒙特卡洛仿真的模糊层次分析法的有效性。

2) 研究结果验证

从国际相关研究结果来看：风和能见度在国际上一直被认为是影响船舶航行安全的重要环境风险因素^[16]，海冰密集度作为极地水域典型环境风险因素也在IMO极地规则^[6]中得到了重点分析。本研究基于蒙特卡洛仿真的模糊层次分析法(MC-FAHP)计算得出的结果与国际研究结论相一致，验证了研究结论的合理性。

4 结论

1) 本研究提出了一种基于蒙特卡洛仿真的模型层次分析法(MC-FAHP)，识别出风、能见度、海冰密集度和气温是影响北极水域船舶航行安全的重要环境风险因素，并从数据和模型两方面对本研究北极水域船舶航行环境风险因素识别研究的合理性进行了分析与验证。

2) 本研究识别出来的环境风险影响因素可以为北极水域船舶事故风险建模研究中的模型参数的选取提供支撑，也可运用于北极水域船舶航行风险评价问题的权重计算。