CN108805350B - 基于多维蒙特卡洛理论的搜救范围预测方法 - Google Patents

Abstract

本发明公开一种基于多维蒙特卡洛理论的搜救范围预测方法，包括以下步骤：考虑风、浪、流的共同作用，建立遇险目标的漂移轨迹预测模型；将遇险目标作为满足马尔可夫过程的单个随机粒子，在时间、空间和风漂系数以及预报风浪流三个维度上，采用Monte‑Carlo方法利用所述漂移轨迹预测模型计算每个随机粒子的最终漂移位置；对大样本容量的随机粒子轨迹进行结果统计分析，确定遇险目标概率最高的搜救范围。

Description

基于多维蒙特卡洛理论的搜救范围预测方法

技术领域

本发明涉及海上搜救技术领域，具体而言，涉及一种针对海上遇险目标的、基于多维蒙特卡洛理论的搜救范围预测方法。

背景技术

近年来，随着海洋开发规模的扩大，我国海上对外贸易、旅游业、水产养殖业、资源开采勘探、科学研究等一系列涉海活动的空前发展，涉海人员数量也迅速增多。日益频繁的海上活动极大地增加了海上事故发生的频率，海难事故更为人们所关注。海上搜救行动能否在事故发生后迅速有效地搜寻到搜救目标并实施救助，对保障人员生命财产安全具有重要的意义。

在整个救援行动中，最关键的就是确定搜寻失事目标的范围，该范围应该包含所有可能存在幸存者的位置。准确的搜索区域不仅要能以最大概率包含搜救对象、不遗漏搜救对象，并且要求搜救区域确定得尽量详细，尽量小，使搜救力量可以集中在最短期内搜索可能性最高的区域。对于搜救计划的决策者来说，海上遇险海难事故发生后，发生的具体位置、具体时间以及搜救的对象类型和样本数都是不确定的，例如申请号为201610136114.8的专利中将这几种关键的条件看作是确定的，直接输入其所发明的“一种改进的海上搜救区域预测方法及系统”中，这显然是不太合理的，因为在海难事故发生时，所有关键的信息中大多都不是确定的，需要综合考虑这些不确定的变量，整合所有已知的信息将变量整合优化，这样才能满足实际运用的需要，更好地服务于整个搜救活动。

发明内容

本发明提供一种基于多维蒙特卡洛理论的搜救范围预测方法，针对现有漂移模型的不足，考虑了时间、空间、漂浮物的风漂系数及预报的风、浪、流场等数据的波动性，在多个维度上采用蒙特卡洛估计，对模型进行优化处理，使其具有更好的实际意义。

本发明提供一种基于多维蒙特卡洛理论的搜救范围预测方法，包括以下步骤：

S1：根据历史气象、海洋数据建立气象、海洋数据库；

S2：通过海上失事遇险第一时间已知信息，获取遇险事故发生地点坐标或地点坐标范围，发生时间或时间范围，搜救对象类型、现场实时观测的有可能残缺或者有一定误差的气象、海洋数据信息；

S3：根据所获取的已知失事信息，并从气象、海洋数据库中调取失事区域附近的气象、海洋数据，将两者数据进行耦合校正并对残缺信息进行初步处理，从而建立初始数据库；根据所创建的基于三维Monte-Carlo理论的搜救范围预测模型矩阵分别考虑海洋环境动力数据、失事时间、失事地点的范围，在保证失事目标最大概率覆盖的情况下，通过Monte-Carlo理论对海上失事目标的三维数据，利用统计实验的随机模拟，得到基于观测数据和历史数据的补充数据，从而建立完整数据库；

S4：根据海洋动力模型、失事时间、失事地点的相互关系，利用所述基于三维Monte-Carlo理论的搜救范围预测模型矩阵对S3所得到的完整数据库进行优化处理；

S5：将S4所优化的数据库作为模拟失事目标的三维数据样本，在失事目标漂移轨迹满足马尔科夫模型理论的基础上，将失事遇险目标抽象优化为单个粒子，在模拟时作为搜救范围预测模型矩阵的信息载体,对优化的数据库中的样本粒子，运用以四阶的Runge-Kutta迭代算法为核心的漂移轨迹预测计算模型，对S4所优化的数据库中的粒子进行漂移轨迹模拟；

S6：利用所提出的非结构化网格数据映射快速寻址技术，遵从低精度结构网格定位、高精度非结构网格计算的主要流程，基于海洋、气象动力模型模拟计算遇险目标漂移轨迹；

S7：确定一定的时间步长，利用S6所述的新技术，将所得的每一条漂移轨迹的最终停留位置在海洋、气象动力模型中进行可视化直观反映；

S8：将S7中得到的可视化的所有粒子运用凸包算法模拟计算出包含所有粒子的搜救最佳范围；

S9：对S8所得到的基于凸包算法的凸包区域进行几何图形上的优化，最终得到更加优化、符合实际搜救情景的搜救范围。

本发明实现了以下有益效果：

本发明基于多维蒙特卡洛理论对遇险目标的搜救范围进行预测，针对现有漂移模型的不足，考虑了时间、空间、漂浮物的风漂系数及预报的风、浪、流场等数据的波动性，在多个维度上采用蒙特卡洛估计，对模型进行优化处理，使其具有更好的实际意义。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为基于多维蒙特卡洛理论的搜救范围预测方法流程图；

图2为基于三维Monte-Carlo理论的搜救范围预测模型矩阵示意图；

图3为根据Monte-Carlo理论对失事目标的三维数据处理过程示意图；

图4为FVCOM模型网格与SIFOM模型网格嵌套覆盖示意图；

图5为FVCOM三角网格中定义SIFOM模型中的垂向坐标示意图；

图6为某凸包图形优化为四边形搜救区域过程示意图；

图7为遇险目标粒子搜救范围计算示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明主要针对现有漂移搜救模型的不足，考虑了在海上遇险事故发生后的短时间内时间、空间、漂浮物的风漂系数及预报的风、浪、流场等数据的波动性，在多个维度上采用蒙特卡洛估计方法，考虑更加全面的各种不确定因素，对模型进行优化处理，使其具有更好的实际意义，为搜救初期建立最大限度的失事目标覆盖区域提供支持。

具体方案如下：一种基于多维蒙特卡洛理论的搜救范围预测方法，包括以下步骤：

一种基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，包括以下步骤：

S1：根据历史气象、海洋数据建立气象、海洋数据库；

S2：通过海上失事遇险第一时间已知信息，获取遇险事故发生地点坐标或地点坐标范围，发生时间或时间范围，搜救对象类型、现场实时观测的有可能残缺或者有一定误差的气象、海洋数据信息；

S3：根据所获取的已知失事信息，并从气象、海洋数据库中调取失事区域附近的气象、海洋数据，将两者数据进行耦合校正并对残缺信息进行初步处理，从而建立初始数据库；根据所创建的基于三维Monte-Carlo理论的搜救范围预测模型矩阵分别考虑海洋环境动力数据、失事时间、失事地点的范围，在保证失事目标最大概率覆盖的情况下，通过Monte-Carlo理论对海上失事目标的三维数据，利用统计实验的随机模拟，得到基于观测数据和历史数据的补充数据，从而建立完整数据库；

S4：根据海洋动力模型、失事时间、失事地点的相互关系，利用所述基于三维Monte-Carlo理论的搜救范围预测模型矩阵对S3所得到的完整数据库进行优化处理；

S5：将S4所优化的数据库作为模拟失事目标的三维数据样本，在失事目标漂移轨迹满足马尔科夫模型理论的基础上，将失事遇险目标抽象优化为单个粒子，在模拟时作为搜救范围预测模型矩阵的信息载体,对优化的数据库中的样本粒子，运用以四阶的Runge-Kutta迭代算法为核心的漂移轨迹预测计算模型，对S4所优化的数据库中的粒子进行漂移轨迹模拟；

S6：利用所提出的非结构化网格数据映射快速寻址技术，遵从低精度结构网格定位、高精度非结构网格计算的主要流程，基于海洋、气象动力模型模拟计算遇险目标漂移轨迹；

S7：确定一定的时间步长，利用S6所述的新技术，将所得的每一条漂移轨迹的最终停留位置在海洋、气象动力模型中进行可视化直观反映；

S8：将S7中得到的可视化的所有粒子运用凸包算法模拟计算出包含所有粒子的搜救最佳范围；

S9：对S8所得到的基于凸包算法的凸包区域进行几何图形上的优化，最终得到更加优化、符合实际搜救情景的搜救范围。

上述方法中，S3所述的创新性的建立基于三维Monte-Carlo理论的搜救范围预测模型矩阵具体如下：

在搜救实例中，假定当满足以下条件(1)、(2)、(3)时，对于搜救目标的确切位置是可以确定的。

(1)确定失事目标失事位置；

(2)确定失事目标失事时间；

(3)确定失事目标失事时的环境数据库，环境数据库包括海洋风、浪、流场情况及漂移物的风漂系数等参数；

根据相关文献资料以及国内外搜救实例，当以上三个条件满足时，在一定的精度内搜救目标的位置是可以确定的。当然这其中排除一些人为的因素，以及环境不确定因素的影响。以上三个参数也是在分析大量搜救案例的基础上，对影响目标区域的多种因素中整合、总结出来的，具有相对独立性和代表性的有效变量。

Monte-Carlo法具有逼真地描述具有随机性质的事物以及物理实验过程的特点，而海上漂浮物的漂移具有很大的随机性，故Monte-Carlo理论非常适合该类问题的求解。针对以上假设，创建基于上述三个因素的三维Monte-Carlo理论的搜救范围预测模型矩阵。

现提出，以海上目标遇险时间T(t)，海面二维坐标L(l)，海洋动力条件H(h)为三个广义的维度建立三维数组Ri,j,k，称为搜救范围预测模型矩阵，在海上事故发生后的短时间内，能够最大限度地覆盖海上遇险事故的各种情景。

设有三维数组R_i,j,k＝[t,l,h]，如简图2所示。

其中t＝{海上目标遇险时间}，t为一维时间信息数据；

l＝{海上目标遇险时位置坐标}，l为二维位置坐标信息数据；

h＝{海上目标遇险时间附近的海洋动力条件}，海洋动力条件h为多维集合的环境数据信息。利用此基于三维Monte-Carlo理论的搜救范围预测模型矩阵可以归纳出几乎所有的海上失事遇险类型：

(1)在已知遇险时间、地点的情况下，基于海上风、海流、海浪等海洋动力环境数值预报技术，构建海上搜救预报模型，这其中的相关海洋动力参数通过Monte-Carlo理论进行选择率定，简称基于海洋环境动力信息的蒙特卡洛模型，即为第一维度的蒙特卡洛模型；

(2)较为常见的搜救类型是在已知航线但未知遇险时间的情况下，此时基于Monte-Carlo理论对失事时间进行预估，简称基于失事时间信息的蒙特卡洛模型，此为第二维度的蒙特卡洛模型；

(3)当接收到搜救信息，但信息中未知失事的位置范围或者位置范围较大，通过Monte-Carlo理论对搜救范围进行预测，简称基于失事位置范围信息的蒙特卡洛模型，此为第三维度的蒙特卡洛模型。

当然实际运用中情况极其复杂，海洋环境动力信息、失事时间、失事地点之间有着千丝万缕的关系，将上述的三维乃至多维的Monte-Carlo理论联系起来综合分析出较为合理的海上搜救模型；通过上述方式，本发明能够为海洋搜救提供较为合理且切实可行的技术支持。

上述方法中，S3、S4中利用Monte-Carlo方法的对失事目标粒子信息的处理过程(如图3所示)，具体如下：

一、获取失事目标信息

海上遇险事故发生后，根据第一时间得到的失事目标海量数据信息，并对该信息进行归纳整理、剔除对于海洋模型要素无用的信息；遴选出具有相对独立性和代表性的有效变量，即S3所述的创新性的建立基于三维Monte-Carlo理论的搜救范围预测模型矩阵中的变量：目标遇险时间T(t)，海面二维坐标L(l)，海洋动力条件H(h)。根据以上三维数据建立搜救范围预测模型矩阵：

R_i,j,k＝[t,l,h]

目标遇险时间T(t)是狭义的一维数据，T(t)∈{t|t＝t₀+dε}；dε满足高斯分布；

海面二维坐标L(l)为狭义的二维数据：此二维数据根据不同的情景，满足相应的基于笛卡尔坐标系的某种解析几何分布，在蒙特卡洛搜救范围预测模型矩阵里为定义为广义的第二维度的数据；

海洋环境动力条件H(h)为多维环境数据的耦合，包括失事区域附近的气象、海洋数据，具体为遇险海域的风浪流数据、风漂系数、漂移目标的几何特征等数据，在蒙特卡洛搜救范围预测模型矩阵里为定义为广义的第三维度的数据。

根据海上遇险事故具体情况的不同，目标失事遇险的第一时间所获取的三维数据有较多的残缺，且每一类数据均有一定的波动范围，即所获取的数据精度无法确定。

二、补全遇险信息

根据第一步所建立的残缺的、不完整的搜救范围预测模型矩阵，对相关数据进行基于Monte-Carlo理论方法的数据扰动：即对残缺的数据进行补全，对较大的波动范围进行优化处理、使其满足特定的数理特征，得到具有完整信息的三维搜救范围预测模型矩阵：

R_i,j,k＝[t_n,l_n,h_n]

该扰动根据不同的信息类别特点，所依据的数学模型将有所差异，例如一维的时间数据扰动，可以依据高斯分布；二维的位置数据可以根据不同的遇险失事类型，依据相应的解析几何图像分布进行扰动；多维的环境数据耦合则需要根据不同的环境数据分量，进行针对性的扰动。此处Monte-Carlo理论的具体数学方法及步骤具体如下：

(1)将所述漂移轨迹预测模型描述成概率过程

设y＝g(ξ₁,ξ₂,...ξ_n,)是n个相互独立的随机变量ξ_i的函数，其概率分布密度为并设定y的数学期望θ＝E(g(ξ₁,ξ₂,...,ξ_n)|(ξ₁,ξ₂,...,ξ_n))为所述漂移轨迹预测模型的真解；

(2)实现从已知概率分布抽样

对所述概率过程进行数字模拟，从f_i(ξ_i)(i＝1,2,…,n)中随机抽样，产生随机变量ξ_i的一个值ξ_i′,i＝1,2,…,n，计算出随机变量y的一个样本值y＝g(ξ₁′,ξ₂′,...ξ_n′)，将此过程重复计算N次，就得到y的N个样本值y₁,y₂,...,y_N；

(3)给出所述漂移轨迹预测模型的解

用θ(y₁),θ(y₂),...,θ(y_N)的均值

作为θ的解。

三、优化失事遇险目标信息

由于第二步所建立的完整的搜救范围预测模型矩阵在补全三个矩阵分量的过程中仅仅局限于各自分量的信息，没有考虑三个分量之间的约束性信息，造成的大量冗杂数据且存在部分不符合实际应用的情况，故需根据三维数据之间的数理关系对相应数据进行选择性删除而进行优化。得到优化的搜救范围预测模型矩阵：

R_i,j,k＝[t_n1,l_n1,h_n1]

此完整的优化的搜救范围预测模型矩阵即可作为能模拟失事目标的三维数据样本，将失事遇险目标抽象优化为单个粒子，在模拟时作为搜救范围预测模型矩阵的信息载体。

上述方法中，S5所述的运用四阶的Runge-Kutta迭代算法预测其漂移轨迹具体如下：

对于遇险目标漂移轨迹预测模型是在风漂作用的基础上，充分考虑风、浪、流的共同作用，建立漂移轨迹预测模型。漂移轨迹预测模型表达式为：

式中：为目标位置，

为风漂速度，为水体流速，包括海面流速和波致Stokes漂流流速，为因为湍流影响发生的随机走动距离。在获取遇险目标初始位置的前提下，给定风场、流场以及波浪场，可以计算遇险目标的漂移轨迹。为提高轨迹计算精度与速度，在进行轨迹模拟计算时，用Runge-Kutta时间积分法对漂移模型[2.1]进行数值求解，四阶Runge–Kutta法如下所示：

其中，f(x)为待求非线性微分方程，x_n为原函数自变量；下一个值(y_n+1)由现在的值(y_n)加上时间间隔(h)和一个估算的斜率的乘积所决定。k₁是时间段开始时的斜率；k₂是时间段中点的斜率，通过欧拉法采用斜率k₁来决定y在点

的值；k₃也是中点的斜率，但是这次采用斜率k₂决定y值；k₄是时间段终点的斜率，其y值用k₃决定。

上述方法中，S5所述的马尔科夫模型(Hidden Markov Model，HMM)，其在海洋应急搜救领域的运用具体如下：

本发明方法中的所用到的模型中失事遇险目标可以优化为单个粒子，这些单个粒子的海上漂移过程满足马尔科夫过程，即目标未来位置的条件概率分布仅依赖于当前状态，与到达当前状态的历史路径无关(在对漂移模型[2.1]进行数值求解时所需遵守的)，其概率分布为：

p(x_n+1|x_n，x_n-1，x_n-2，...，x₁)＝p(x_n+1|x_n) (1)

其中，x_i(i＝1,2,3,...n...)表示第i时刻目标的位置，p表示条件概率。单个粒子的随机漂移过程可表示为：

Δx＝V_o(n)ΔT+dε (2)

x_n+1＝x_n+Δx (3)

其中，V_o(n)为第n时刻目标的漂移速度，ΔT为离散时间步长，dε为位置的随机扰动项，该扰动项符合高斯分布；x_n表示第n时刻粒子的位置。

上述方法中，S5中对S4所优化的数据库中的粒子进行漂移轨迹模拟，具体如下：

根据Monte-Carlo理论方法的基本思想，在实际搜救过程中，把漂浮物认定为单个的粒子，粒子的运动看作是一个马尔可夫随机过程。对漂浮物在某个位置向下一位置的漂移概率进行建模，在同一规则下假设这些粒子的运动是相互独立的，并大量复制这些满足条件的粒子，依据这些粒子在一定时间步长Δt时间后的分布的统计属性，即可分析出漂浮物漂移运动随时间变化的概率区域，计算过程中增加模拟次数，使得漂浮物运行的概率区域更准确。此处步骤具体如下：

1)生成n个初始位置满足遇险事故位置分布的随机粒子，用以表示符合漂流试验的样本。

2)对每个样本按照漂流模型进行计算，得到该样本的最终漂流位置。

(1)确定当前样本粒子的初始位置，该位置由Monte-Carlo理论方法获得，满足遇险事故位置概率分布。

(2)在第一步仿真步长时段，该样本受到所在区域总流压的作用，沿流压矢量方向漂移，流压矢量生成由环境数据库提供，同时受到区域风力作用，产生风压效果，计算风压漂移距离，确定该步长漂移位置。

(3)在未达到预测时间之前，重复上一步，每次使用当前所在位置海区时段环境数据更新计算数据。

(4)到达预测时间，输出该样本最终位置。

3)对包含所有试验结果的区域划分成若干区域，统计该漂浮物在该区域的包含概率。

在整个区域内，所包含的粒子数密集程度代表该区域包含密度大小。

上述方法中，S6所述研发的非结构化网格数据映射快速寻址技术，具体如下：

相对于结构矩形网格，非结构三角网格模拟效果较好，尤其是非结构三角网格的局部加密方案可以进一步提高加密区域的有效模拟精度，但是由于非结构三角网格模型较为复杂，计算速度、效率较低，在目前的所存在的技术中，高精度的非结构三角网格模型在实际的大范围应用中效果欠佳。研发非结构化网格数据映射快速寻址技术，使遇险目标漂移轨迹基于风浪流数值场的计算效率可以满足一定的应用需求，即在实际搜救业务中实现快速预报。具体寻址模型如下：

第一步——寻址模型构建：

对目标海域进行低精度的结构化网格划分。为了增加基于有限差分格式建立的模型的计算效率，对目标海域Ω大尺度地划分为多个区块，并将所划分的区块按顺序进行编号1,2,3,…,n；显然满足

同时将每个结构化网格单元作为一个数据存储单元，记作K₁,K₂,K₃,…,K_n，此时的数据存储单元已经具备了前述的搜救范围预测模型矩阵R_i,j,k中第三维海洋环境的部分数据。在每个大尺度的结构化网格单元内根据单元内不同的海洋、气象条件，对敏感性的区域的计算网格进行非结构化网格的局部加密处理。

第二步——寻址定位：

当事故发生后，根据已知的信息确定失事位置为以r₁为半径的圆形范围Θ₁，并根据历史海洋、气象数据预估范围的半径精度r₂，可以得到目标区域Θ₂，且

根据Θ₂的位置信息，对上述结构化网格单元遍历，即可快速定位在某些区块，将定位区域附近作空间上的正态分布概率分析：确定以预估地点为圆心、以一定的精度长度为半径的覆盖所有目标区块的圆形区域。

第三步——精确计算：

在大尺度网格存储单元内生成的非结构化三角网格，可以根据该网格内特殊的海洋气象条件、固体边界条件对三角网格进行常规的局部的加密处理。本发明提出运用网格嵌套技术，在此运用中该网格嵌套系统由一个全三维流体力学模型(SIFOM)和一个地球物理流体力学模型(FVCOM)组成。

通过嵌套网格的方法，可以将SIFOM网格与FVCOM网格任意嵌套，如图4所示，图中填充圆是SIFOM和FVCOM的接口节点，开圆(元素质心)、三角形(网格节点)和正方形(元素垂直边的中间点)是FVCOM的接口节点。在混合的网格嵌套系统中，SIFOM完全解决了三维小尺度的局部流现象，而FVCOM捕获了大量的背景流，通过区域分解实现了这两个模型的集成耦合，其中两个网格之间的数据的互联是此类网格嵌套的难点所在，本发明提出以下两步解决两种网格的互联：

(1)在水平平面上的三角形网格中搜索，如果一个三角形满足以下条件，则三角与矩形单元格相关联：

β₃₁×β₃₄＞0,β₁₂×β₁₄＞0,β₂₃×β₂₄＞0

其中β是连接其下标所表示的两个点的向量。例如，β₂₁是从2起点以1终点的向量。

(2)在垂直方向上搜索如图5所示，将确定接口节点M在不同σ层(如4和4′、d和d′)上的投影坐标。显然，这些投影的x和y坐标将与M的坐标相同，它们的z坐标将通过插值表示。例如4点的z坐标计算为

A表示为垂直投影方向对应的三角形面积，z为其竖坐标。

以上方法在充分发挥非结构化网格数学模型优势的基础上，通过采用网格嵌套技术，能够很大程度上减小计算域中网格的数量，提高了计算速度和效率。同时，可以布置更加精细的计算网格，以适应复杂的地形、曲折的岸线边界和大范围的海流和风的数据值模拟，提高计算结果的准确性。

上述方法中，S8中将S7所得到的所有单粒子进行可视化展示，运用凸包算法找出搜救的最佳范围，具体如下：

上述模型运行一定时间后，以大量的粒子最终停留的位置为基点，此处引入拓扑几何学中的凸包算法，运用Graham扫描法包围所有粒子的位置点，形成的位置凸包即为搜寻区域。

上述方法中，S9中对S8所得到的基于凸包算法的凸包区域进行几何图形上的优化，可以得到更加合适的、优化的搜救范围，具体如下：

在实际的搜救事件中，当失事遇险目标较多、或者目标确定的参数较少时，由Monte-Carlo理论对三维参数扰动模拟后产生的数据量巨大，这样完全依照凸包算法得到的搜寻区域形状及其不规则，不利于快速搜寻。鉴于此，对于最终确定的搜寻区域，并不是简单地将所有粒子的凸包范围作为搜寻区域，而是将包围所有粒子凸包的最小四边形所围区域定为搜寻区域，如图6所示，即将所得到的粒子凸包优化为四边形，这样在实际搜救中更有可操作性。

实施例1

结合图1，表示出了本发明所提出的基于多维蒙特卡洛理论对遇险目标的搜救范围预测方法流程，具体方案如下：

1、通过海上失事遇险第一时间已知信息，获取遇险事故发生地点坐标或地点坐标范围，发生时间或时间范围，搜救对象类型、现场实时观测的有可能残缺或者有一定误差的气象、海洋数据等信息；

2、根据所获取的已知失事信息，对已得的大量信息进行归纳整理并对残缺信息进行智能化处理；继而从气象、海洋数据库中调取失事区域附近的气象、海洋数据与之耦合校正，从而建立初始数据库。

建立初始数据库，现以最常见的海上遇险情况为例，说明对获取的第一手资料进行智能化处理。

“MH370”事件是在特定航线上航行(由吉隆坡国际机场飞往北京首都国际机场)，事件中对“MH370”失事遇险的最终位置并没有确定。其中在构建模型的预处理时，满足以下条件：

l＝f(t)........................................(1)

l_n～N(0，δ²)................................(2)

h＝g(t，l).....................................(3)

(1)式表示失事位置与最后失联时间之间的关系，这个关系可以根据航空公司提供的时间、航线信息确定；(2)式为根据(Monte-Carlo参数扰动)估计失事位置，即在特定段的航线两边产生多个满足条件的粒子，其最终失事位置概率分布满足：以行经的航线为对称轴的正态分布；(3)式根据蒙特卡洛理论得出多个满足条件的粒子的位置以及所匹配的时间，根据查阅公开的海洋环境动力数据。

3、分别考虑初始数据库中海洋环境动力数据、失事时间、失事地点的范围，在保证失事目标最大概率覆盖的情况下，创建的基于三维Monte-Carlo理论的搜救范围预测模型矩阵R_i，j，k＝[t，l，h]，如图2所示；

4、通过Monte-Carlo理论对海上失事目标的三维数据，利用统计实验的随机模拟，得到基于观测数据和历史数据的补充数据，从而建立完整数据库；此处Monte-Carlo理论的具体数学方法及步骤具体如下：

(1)获取失事目标信息

海上遇险事故发生后，根据第一时间得到的失事目标海量数据信息，并对该信息进行归纳整理、剔除对于海洋模型要素无用的信息；遴选出具有相对独立性和代表性的有效变量，即建立基于三维Monte-Carlo理论的搜救范围预测模型矩阵中的变量：目标遇险时间T(t)，海面二维坐标L(l)，海洋动力条件H(h)。根据以上三维数据建立搜救范围预测模型矩阵：

R_i,j,k＝[t,l,h]

目标遇险时间T(t)是狭义的一维数据，T(t)∈{t|t＝t₀+dε}；dε满足高斯分布；

海面二维坐标L(l)为狭义的二维数据：此二维数据根据不同的情景，满足相应的基于笛卡尔坐标系的某种解析几何分布，在蒙特卡洛搜救范围预测模型矩阵里为定义为广义的第二维度的数据；

海洋环境动力条件H(h)为多维环境数据的耦合，包括失事区域附近的气象、海洋数据，具体为遇险海域的风浪流数据、风漂系数、漂移目标的几何特征等数据，在蒙特卡洛搜救范围预测模型矩阵里为定义为广义的第三维度的数据。根据海上遇险事故具体情况的不同，目标失事遇险的第一时间所获取的三维数据有较多的残缺，且每一类数据均有一定的波动范围，即所获取的数据精度无法确定。

(2)补全遇险信息

根据第一步所建立的残缺的、不完整的搜救范围预测模型矩阵，对相关数据进行基于Monte-Carlo理论方法的数据扰动，即对残缺的数据进行补全，对较大的波动范围进行优化处理、使其满足特定的数理特征，得到具有完整信息的三维搜救范围预测模型矩阵：

R_i,j,k＝[t_n,l_n,h_n]

该扰动根据不同的信息类别特点，所依据的数学模型将有所差异，例如一维的时间数据扰动，可以依据高斯分布；二维的位置数据可以根据不同的遇险失事类型，依据相应的解析几何图像分布进行扰动；多维的环境数据耦合则需要根据不同的环境数据分量，进行针对性的扰动。此处Monte-Carlo理论的具体数学方法及步骤具体如下：

1)将所述漂移轨迹预测模型描述成概率过程

设y＝g(ξ₁,ξ₂,…ξ_n,)是n个相互独立的随机变量ξ_i的函数，其概率分布密度为f_i(ξ_i)(i＝1,2,…,N)，并设定y的数学期望θ＝E(g(ξ₁,ξ₂,...,ξ_n)|(ξ₁,ξ₂,...,ξ_n))为所述漂移轨迹预测模型的真解；

2)实现从已知概率分布抽样

对所述概率过程进行数字模拟，从f_i(ξ_i)(i＝1,2,…,n)中随机抽样，产生随机变量ξ_i的一个值ξ_i′,i＝1,2,…,n，计算出随机变量y的一个样本值y＝g(ξ₁′,ξ₂′,…ξ_n′)，将此过程重复计算N次，就得到y的N个样本值y₁,y₂,…,y_N；

3)给出所述漂移轨迹预测模型的解

用θ(y₁),θ(y₂),…,θ(y_N)的均值作为θ的解。

(3)优化失事遇险目标信息

由于第二步所建立的完整的搜救范围预测模型矩阵在补全三个矩阵分量的过程中仅仅局限于各自分量的信息，没有考虑三个分量之间的约束性信息，造成的大量冗杂数据且存在部分不符合实际应用情况，故需根据三维数据之间的数理关系对相应数据进行选择性删除而进行优化。得到优化的搜救范围预测模型矩阵：

R_i,j,k＝[t_n1,l_n1,h_n1]

此完整的优化的搜救范围预测模型矩阵即可作为能模拟失事目标的三维数据样本，将失事遇险目标抽象优化为单个粒子，在模拟时作为搜救范围预测模型矩阵信息载体。具体过程如图3所示。

5、对完整数据库中的大量粒子运用以四阶Runge-Kutta漂移轨迹预测模型方法，在设置一定的时间步长的情况下，可以计算出其最终漂移位置；在实际运用中需要建立漂移轨迹预测模型。对于遇险目标漂移轨迹预测模型是在风漂作用的基础上，充分考虑风、浪、流的共同作用，建立漂移轨迹预测模型。漂移轨迹预测模型表达式为：

式中：

为目标位置，

为风漂速度，

为水体流速，包括海面流速和波致Stokes漂流流速，

为因为湍流影响发生的随机走动距离。在获取遇险目标初始位置的前提下，给定风场、流场以及波浪场，可以计算遇险目标的漂移轨迹。为提高轨迹计算精度与速度，在进行轨迹模拟计算时，用Runge-Kutta时间积分法对漂移模型[2.1]进行数值求解，四阶Runge–Kutta法如下所示：

其中，f(x)为待求非线性微分方程，x_n为原函数自变量；下一个值(y_n+1)由现在的值(y_n)加上时间间隔(h)和一个估算的斜率的乘积所决定。k₁是时间段开始时的斜率；k₂是时间段中点的斜率，通过欧拉法采用斜率k₁来决定y在点

的值；k₃也是中点的斜率，但是这次采用斜率k₂决定y值；k₄是时间段终点的斜率，其y值用k₃决定。

然而计算过程中存在着很多的不确定因素和误差，将导致计算的漂移轨迹与真实漂移轨迹随着时间推移不断偏离，因此在这里继续引入Monte-Carlo算法。在漂移模型中，其基本思想是将漂浮物定义为单个粒子，粒子的移动轨迹是一个满足一定概率分布的随机过程，在同一规则下，大量复制这些粒子，并设这些粒子的运动互不相关，那么这些粒子在一段时间后的统计区域即是漂浮物漂移的概率区域。将在海面上的失事目标看做满足马尔科夫模型，其马尔科夫模型即目标未来位置的条件概率分布仅依赖于当前状态，与到达当前状态的历史路径无关，其概率分布为：

p(x_n+1|x_n，x_n-1，x_n-2，...，x₁)＝p(x_n+1|x_n) (1)

其中，x_i(i＝1,2,3,…n…)表示第i时刻目标的位置，p表示条件概率。单个粒子的随机漂移过程可表示为：

Δx＝V_o(n)ΔT+dε (2)

x_n+1＝x_n+Δx (3)

其中，V_o(n)为第n时刻目标的漂移速度，ΔT为离散时间步长，dε为位置的随机扰动项，该扰动项符合高斯分布；x_n表示第n时刻粒子的位置。

6、利用所研发的非结构化网格数据映射快速寻址技术，以快速得到海上搜救目标漂移轨迹，提高对大量粒子进行预测漂移轨迹预测运算的工作效率；具体操作如下：

第一步——寻址模型构建：

对目标海域进行低精度的结构化网格划分。为了增加基于有限差分格式建立的模型的计算效率，对目标海域Ω大尺度地划分为多个区块，并将所划分的区块按顺序进行编号1,2,3,…,n；显然满足

同时将每个结构化网格单元作为一个数据存储单元，记作K₁,K₂,K₃,…,K_n，此时的数据存储单元已经具备了前述的搜救范围预测模型矩阵R_i,j,k中第三维海洋环境的部分数据。在每个大尺度的结构化网格单元内根据单元内不同的海洋、气象条件，对敏感性的区域的计算网格进行非结构化网格的局部加密处理。

第二步——寻址定位：

当事故发生后，根据已知的信息确定失事位置为以r₁为半径的圆形范围Θ₁，并根据历史海洋、气象数据预估范围的半径精度r₂，可以得到目标区域Θ₂，且S_Θ2＝π(r₁+r₂)²。根据Θ₂的位置信息，对上述结构化网格单元遍历，即可快速定位在某些区块，将定位区域附近作空间上的正态分布概率分析：确定以预估地点为圆心、以一定的精度长度为半径的覆盖所有目标区块的圆形区域。

第三步——精确计算：

在大尺度网格存储单元内生成的非结构化三角网格，可以根据该网格内特殊的海洋气象条件、固体边界条件对三角网格进行常规的局部的加密处理。本发明提出运用网格嵌套技术，在此运用中该网格嵌套系统由一个全三维流体力学模型(SIFOM)和一个地球物理流体力学模型(FVCOM)组成。

通过嵌套网格的方法，可以将SIFOM网格与FVCOM网格任意嵌套，如图4所示，图中填充圆是SIFOM和FVCOM的接口节点，开圆(元素质心)、三角形(网格节点)和正方形(元素垂直边的中间点)是FVCOM的接口节点。在混合的网格嵌套系统中，SIFOM完全解决了三维小尺度的局部流现象，而FVCOM捕获了大量的背景流，通过区域分解实现了这两个模型的集成耦合。其中两个网格之间的数据的互联是此类网格嵌套的难点所在，本发明专利提出以下两步解决两种网格的互联：

(1)在水平平面上的三角形网格中搜索，如果一个三角形满足以下条件，则三角与矩形单元格相关联：

β₃₁×β₃₄＞0,β₁₂×β₁₄＞0,β₂₃×β₂₄＞0

其中β是连接其下标所表示的两个点的向量。例如，β₂₁是从2起点以1终点的向量。

(2)在垂直方向上搜索如图5所示，将确定接口节点M在不同σ层(如4和4′、d和d′)上的投影坐标。显然，这些投影的x和y坐标将与M的坐标相同，它们的z坐标将通过插值表示。例如4点的z坐标计算为

A表示为垂直投影方向对应的三角形面积，z为其竖坐标。

根据以上方法，利用“MH370”在特定航线上失事遇险的时间的波动性(Monte-Carlo参数扰动)可以估计失事位置，进而估计其漂移路径上的海洋动力条件，从而对其漂移轨迹和搜救范围进行预测。

6、上述模型运行一定时间后，以大量的粒子最终停留的位置为基点，此处引入拓扑几何学中的凸包算法，运用Graham扫描法包围所有粒子的位置基点，形成的位置凸包即为搜寻区域。

7、对于最终确定的搜寻区域，并不是简单地将所有粒子的凸包范围作为搜寻区域，而是将包围所有粒子凸包的最小四边形所围区域定为搜寻区域，如图6所示，即将所得到的粒子凸包优化为四边形，这样在实际搜救中更有可操作性。对于一个已经获得所有三维信息的单个粒子，其模型整体运行示意图如图7所示。

与现有技术相比，本发明实现了以下有益效果：

本发明基于多维蒙特卡洛理论对遇险目标的搜救范围进行预测，针对现有漂移模型的不足，考虑了时间、空间、漂浮物的风漂系数及预报的风、浪、流场等海洋、气象数据的波动性(不同于以往仅仅考虑风浪流数据的波动性)，在多个维度上采用蒙特卡洛估计，考虑了更加全面的影响因素，对模型进行优化处理，使其具有更好的实际意义，为搜救初期建立最大限度的失事目标覆盖区域提供支持。

以上说明对本发明而言只是说明性的，而非限制性的，本领域普通技术人员理解，在不脱离以下所附权利要求所限定的精神和范围的情况下，可做出许多修改、变化或等效，但都将落入本发明的保护范围内。

Claims (8)

1.一种基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，包括以下步骤：

S1：根据历史气象、海洋数据建立气象、海洋数据库；

S2：通过海上失事遇险第一时间已知信息，获取遇险事故发生地点坐标或地点坐标范围，发生时间或时间范围，搜救对象类型、现场实时观测的有可能残缺或者有一定误差的气象、海洋数据信息；

S3：建立完整数据库，具体处理过程如下：

一、获取失事目标信息

海上遇险事故发生后，根据第一时间得到的失事目标海量数据信息，并对该信息进行归纳整理、剔除对于海洋模型要素无用的信息；遴选出具有相对独立性和代表性的有效变量，即：目标遇险时间T(t)，海面二维坐标L(l)，海洋动力条件H(h)；根据以上三维数据建立搜救范围预测模型矩阵：

R_i,j,k＝[t,l,h]

目标遇险时间T(t)是狭义的一维数据，T(t)∈{t|t＝t₀+dε}；dε满足高斯分布；

海面二维坐标L(l)为狭义的二维数据：此二维数据根据不同的情景，满足相应的基于笛卡尔坐标系的解析几何分布，在蒙特卡洛搜救范围预测模型矩阵里为定义为广义的第二维度的数据；

海洋环境动力条件H(h)为多维环境数据的耦合，包括失事区域附近的气象、海洋数据，具体为遇险海域的风浪流数据、风漂系数、漂移目标的几何特征数据，在蒙特卡洛搜救范围预测模型矩阵里为定义为广义的第三维度的数据；

二、补全遇险信息

根据第一步所建立的残缺的、不完整的搜救范围预测模型矩阵，对相关数据进行基于Monte-Carlo理论方法的数据扰动：即对残缺的数据进行补全，对较大的波动范围进行优化处理、使其满足特定的数理特征，得到具有完整信息的三维搜救范围预测模型矩阵：

R_i,j,k＝[t_n,l_n,h_n]

该扰动根据不同的信息类别特点，所依据的数学模型将有所差异，一维的时间数据扰动，依据高斯分布；二维的位置数据根据不同的遇险失事类型，依据相应的解析几何图像分布进行扰动；多维的环境数据耦合则需要根据不同的环境数据分量，进行针对性的扰动；

S4：根据海洋动力模型、失事时间、失事地点的相互关系，利用所述基于三维Monte-Carlo理论的搜救范围预测模型矩阵对S3所得到的完整数据库进行优化处理；

S5：将S4所优化的数据库作为模拟失事目标的三维数据样本，在失事目标漂移轨迹满足马尔科夫模型理论的基础上，将失事遇险目标抽象优化为单个粒子，在模拟时作为搜救范围预测模型矩阵的信息载体,对优化的数据库中的样本粒子，运用以四阶的Runge-Kutta迭代算法为核心的漂移轨迹预测计算模型，对S4所优化的数据库中的粒子进行漂移轨迹模拟；

S6：利用所提出的非结构化网格数据映射快速寻址技术，遵从低精度结构网格定位、高精度非结构网格计算的主要流程，基于海洋、气象动力模型模拟计算遇险目标漂移轨迹；

S7：确定一定的时间步长，利用S6所述的非结构化网格数据映射快速寻址技术，将所得的每一条漂移轨迹的最终停留位置进行可视化直观反映；

S8：将S7中得到的可视化的所有粒子运用凸包算法模拟计算出包含所有粒子的搜救最佳范围；

S9：对S8所得到的基于凸包算法的凸包区域进行几何图形上的优化，最终得到更加优化、符合实际搜救情景的搜救范围。

2.根据权利要求1所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，步骤S5中所述的运用四阶的Runge-Kutta迭代算法预测其漂移轨迹具体如下：

在风漂作用的基础上，考虑风、浪、流的共同作用，建立漂移轨迹预测模型，所述漂移轨迹预测模型表达式为：

式中：

为目标位置，

为风漂速度，为水体流速，包括海面流速和波致Stokes漂流流速，

为因为湍流影响发生的随机走动距离；为提高轨迹计算精度与速度，在进行轨迹模拟计算时，用Runge-Kutta时间积分法对漂移模型[2.1]进行数值求解，四阶Runge–Kutta法如下所示：

其中，f(x)为待求非线性微分方程，x_n为原函数自变量；下一个值(y_n+1)由现在的值(y_n)加上时间间隔(h)和一个估算的斜率的乘积所决定；k₁是时间段开始时的斜率；k₂是时间段中点的斜率，通过欧拉法采用斜率k₁来决定y在点

的值；k₃也是中点的斜率，但是这次采用斜率k₂决定y值；k₄是时间段终点的斜率，其y值用k₃决定。

3.根据权利要求2所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，步骤S5中利用马尔科夫模型对S4所优化的数据库中的粒子进行漂移轨迹模拟具体如下：

假定粒子的漂移满足马尔科夫过程，即目标未来位置的条件概率分布仅依赖于当前状态，与到达当前状态的历史路径无关，其概率分布为：

p(x_n+1|x_n，x_n-1，x_n-2，...，x₁)＝p(x_n+1|x_n) (1)

其中，x_i(i＝1,2,3,...n...)表示第i时刻目标的位置，n为自然数，p表示条件概率；单个粒子的随机漂移过程可表示为：

Δx＝V_o(n)ΔT+dε (2)

x_n+1＝x_n+Δx (3)

其中，V_o(n)为第n时刻目标的漂移速度，ΔT为离散时间步长，dε为位置的随机扰动项，该扰动项符合高斯分布；x_n表示第n时刻粒子的位置。

4.根据权利要求3所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，步骤S6中利用所提出的非结构化网格数据映射快速寻址技术，遵从低精度结构网格定位、高精度非结构网格计算的流程，基于海洋、气象数据模型计算遇险目标漂移轨迹包括：

第一步——寻址模型构建：

对目标海域进行低精度的结构化网格划分，为了增加基于有限差分格式建立的模型的计算效率，对目标海域Ω大尺度地划分为多个区块，并将所划分的区块按顺序进行编号1,2,3,…,n；显然目标海域面积满足同时将每个结构化网格单元作为一个数据存储单元，记作K₁,K₂,K₃,…,K_n，此时的数据存储单元已经具备了前述的搜救范围预测模型矩阵R_i,j,k中第三维海洋环境的部分数据，其中，所述第三维海洋环境是指海上目标遇险时间附近的海洋动力条件；在每个大尺度的结构化网格单元内根据单元内不同的海洋、气象条件，对敏感性的区域的计算网格进行非结构化网格的局部加密处理；

第二步——寻址定位：

当事故发生后，根据已知的信息确定失事位置为以r₁为半径的圆形范围Θ₁，并根据历史海洋、气象数据预估范围的半径精度r₂，可以得到目标区域Θ₂，且

根据Θ₂的位置信息，对上述结构化网格单元遍历，即可快速定位在某些区块，将定位区域附近作空间上的正态分布概率分析：确定以预估地点为圆心、以一定的精度长度为半径的覆盖所有目标区块的圆形区域；

第三步——精确计算：

在大尺度网格存储单元内生成的非结构化三角网格，根据该网格内特定的海洋气象条件、固体边界条件对非结构三角网格进行常规的局部的加密处理，并利用全三维流体力学模型(SIFOM)和一个地球物理流体力学模型(FVCOM)两种网格相互嵌套提高计算结果的准确性。

5.根据权利要求1所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，所述将遇险目标作为单个随机粒子，在时间、空间和包括有风漂系数、预报风浪流在内的海洋环境数据三个维度上，采用Monte-Carlo方法利用所述漂移轨迹预测模型计算每个随机粒子的最终漂移位置步骤包括：

(1)将所述漂移轨迹预测模型描述成概率过程：

设y＝g(ξ₁,ξ₂,...ξ_n,)是n个相互独立的随机变量ξ_i的函数，其概率分布密度为f_i(ξ_i)(i＝1,2,…,N)，并设定y的数学期望θ＝E(g(ξ₁,ξ₂,...,ξ_n)|(ξ₁,ξ₂,...,ξ_n))为所述漂移轨迹预测模型的真解；

(2)实现从已知概率分布抽样

对所述概率过程进行数字模拟，从f_i(ξ_i)(i＝1,2,…,n)中随机抽样，产生随机变量ξ_i的一个值ξ_i′,i＝1,2,...,n，计算出随机变量y的一个样本值y＝g(ξ₁′,ξ₂′,...ξ_n′)，将此过程重复计算N次，就得到y的N个样本值y₁,y₂,...,y_N；

(3)给出所述漂移轨迹预测模型的解

用θ(y₁),θ(y₂),...,θ(y_N)的均值作为θ的解。

6.根据权利要求2所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，所述搜救范围通过以下方式确定：

通过所述漂移轨迹预测模型对每个随机粒子的漂移轨迹进行预测，得到每个随机粒子的当前位置；模型运行一定时间后，以大量的粒子最终停留的位置为基点，引入拓扑几何学中的凸包算法，运用Graham扫描法包围所有粒子的位置基点，形成的位置凸包即为搜寻区域。

7.根据权利要求6所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，将包围所有粒子凸包的最小四边形所围区域确定为搜寻区域。

8.根据权利要求6所述的基于多维蒙特卡洛理论的搜救范围预测方法，其特征在于，所述随机粒子的漂移满足的马尔科夫过程用公式描述为：

p(x_n+1|x_n，x_n-1，x_n-2，...，x₁)＝p(x_n+1|x_n) (1)

其中，x_i(i＝1,2,3,...n...)表示第i时刻遇险目标的位置，p表示条件概率；

单个粒子的随机漂移过程可表示为：

Δx＝V_o(n)ΔT+dε (2)

x_n+1＝x_n+Δx (3)

其中，V_o(n)为第n时刻目标的漂移速度，ΔT为离散时间步长，dε为位置的随机扰动项，该扰动项符合高斯分布；x_n表示第n时刻粒子的位置。

Images