魏思涵,焦鹿怀,黎琪,刘俊材,杨剑轩,张兰

(武汉大学 电子信息学院，湖北 武汉 430072)

0 引 言

随着我国海上交通和旅游业的发展，从事水上工作的人员数量急剧增加，同时也导致海难事故遇难人数的急剧上升，因此海上搜救工作显得极为重要。我国目前海上搜救力量的投入主要侧重于装备的更新，还上搜救应急反应速度慢，在搜救区域的确定上仍显不足[1]。目前我国国家海洋局现有的预报系统中，关于漂流物轨迹预报以及数值预报的模型都已成型，数值的获取主要依靠浮标潜标等传统测量工具进行测量以及星载雷达遥感测量[2]。传统测量工具只能得到点或面上海流数据，无法获得大面积实时信息，精度不够高[3]；星载雷达遥感覆盖面积大，空间分辨率较高，但时间分辨率不足，数据率较低，难以实现连续实时海态观测[4]。利用无人机开展海事搜救和巡航监管与船艇巡航搜救相比，从经济性、实用性、高效性方面具有不可比拟的优势。但采用微波传输技术实现现场图像实时回传受微波本身特性限制，直传距离难以超过100 km,且利用无人机进行常态化监管及海事搜救任务须取得飞行许可[5]。高频地波雷达具有作用距离远、监测范围大、实时性好等特点，利用其优势进行海上人员搜救辅助决策系统的研究与开发具有重要研究意义与实用价值[6]。

国外将高频地波雷达用于流场信息探测领域，且已经达到业务化运行程度，如美国CODAR公司已为世界各国提供数百套地波雷达系统用于海洋环境的监测[7]。国内目前也有地波雷达对于海流及波浪场的较为深入的探测研究，暂时仅应用于海上漂流物体的轨迹跟踪，结合海洋动力学及人体特征参数以进行海上漂流人员轨迹的预测的文献较少[8]。

国内对于利用高频地波雷达进行海洋搜救区域算法的研究较少，近年有一些相关的成果发表。2012年，潘军探讨了利用蒙特卡罗方法确定海上搜救区域的思路并进行了误差估计，但并未结合高频地波雷达的实时海态数据进行仿真和算法验证[9]。李伦在2013年利用高频地波雷达对海洋动力学参数进行了反演研究，但对于海上物体漂移运动轨迹的预测论述较少[10]。2016年王道生等在一种改进的海上搜救区域预测方法及系统中，结合现有气象环境数据、三维流场、改进AP98风压模型以及海上搜救目标漂移轨迹计算模型实现对海上搜救对象的动态漂移轨迹进行预测，但该发明需要收集历史气象数据建立气象数据库，并未结合高频地波雷达的实时海态数据[11]。宋军等在2013年提出一种海上搜救目标漂移路径的预测方法，基于海面风、海流、海浪等海洋动力环境场数值预报技术，构建海上搜救预报模型预报落水人员和海上失事船舶的漂移路径[12]。刘建等在2012年发明一种搜救终端和基于该搜救终端的水上搜救方法，利用传感器向船载基站发送求救信号，可以快速定位落水人员[13]。本文旨在利用高频地波雷达返回的实时海态数据，建立动力学模型，模拟海上物体的漂移轨迹，并结合计算机工具，引入蒙特卡罗方法对算法进行仿真并验证其可行性，预测搜救区域，从而对海上救援行动起到辅助决策作用。

1 落水人员漂移轨迹相关影响因素分析

建立漂移模型需要知道物体运动的速度和方向，一般来说，漂移轨迹的影响因素包括环境因素和人体自身的影响。

1.1 流对漂移的影响

随流漂移是物体受周围水流的影响而产生的漂移。为简化运动过程，假设流的速度变化时，人体的速度也随之变化。在无风的状态下，流的速度视作人的漂移速度。

1.2 风对漂移的影响

考虑到海面上漂移物体露出海面的部分是不规则的，不具有对称性。因此，物体在风的影响下，漂移速度除了下风方向外，还会有一个偏风的分量[14]。如图1所示：

图1 风致漂移分量

研究表明，海平面上10 m风的速度与物体的风致漂移速度基本是线性的。即对于下风分量Ld和偏风分量Lc，满足：

Ld=αdW10d

Lc=αcW10c

(1)

其中，W10d和W10c为海平面上10m处风速在下风方向和偏风方向的分量，αd和αc分别是下风方向和偏风方向的风导因子，其取值具有经验性，一般在0.01～0.05之间，可以通过拟合实验数据得到，一般与漂移物体的类型有关。本文根据Allen建立的LeewayFormula Table[15]，取αd=0.017,αc=0.0136。

1.3 浪等其他因素对漂移的影响

波浪的漂流作用力，就是漂移物体反射波浪所产生的作用力。由于人体的实际尺寸远小于波浪的波长因此波浪对人体的漂移轨迹影响可忽略不计。

海洋中的湍流则是漂移运动不确定性的主要来源。由于湍流的实际大小和方向都具有随机性，因此人体的漂移模型也相应的具有随机性，其运动遵循维纳过程[16]。

1.4 人体自身的影响

落水人员的性别，是否穿戴救生衣等状况会影响其浸没比例，而风致漂移速度与人体的浸没比例有关，如图2[17]。

图2 风致漂移速度与浸没比例关系

根据浮力定律：

ρsgvs=ρpgvp

(2)

其中ρs,ρp分别为海水和人体的密度，vs,vp分别为人体的浸没体积和人体的体积。则人体在海水中的浸没比例：

(3)

用户在系统的输入界面输入落水人员的性别，身高，体重，是否穿戴救生衣等信息，系统将自动对人体的浸没比例进行计算，从而实现对风致漂移速度的估计。

2 落水人员漂移运动模型的建立

对漂移运动的影响因素分析可知，物体的漂移运动是风和流引起的随流运动与湍流引起的随机运动的叠加。

2.1 漂移运动方程的导出

根据牛顿运动学定律，海洋表面漂移物体的受力满足：

(4)

其中m是物体的质量，v是物体的运动速度，∑F是物体所受的合外力。附加质量m′由物体周围水粒子受加速度产生[18]。

海上漂移物在海面的运动方程：

(5)

S为位置矢量，v(t)为漂移速度。

漂流速度受流速风速及湍流等影响

v(t)=vc(t)+αW10(t)+β(t)

(6)

vc(t)为海面流矢量，W10(t)为10 m处风速，α为风因子，β(t)涡流的影响产生的随机漂移速度。

物体受到海水的携带和风力的拖曳作用，其结果是表现为在风力和潮流共同作用下向前漂移。设物体的初始位置矢量为S0，经过一段时间Δt后漂移到新的位置S1，根据拉格朗日追踪法[19]：

(7)

其中ξ(t)是物体受湍流影响产生的随机位移。

当Δt取值足够小时，可以将t0时刻物体的速度vt0作为Δt这段时间内物体的运动速度，这样，我们可以设置一系列微小的时间步长，将原本复杂的曲线运动变为一系列匀速运动的叠加。使用迭代的方法，当物体位移到S1之后再将其作为初始位置S0，用在该处的速度计算到再下一个时间步长的位移，以此类推，从而对漂移物的运动进行预测。

漂移物体的运动速度主要受到海流和风的影响，而随机因素的影响则会使搜救区域存在不确定性。分析漂移轨迹的影响因素并结合拉格朗日追踪法可知，漂移模型建立的重点在于计算物体漂移速度并结合随机因素的影响预测漂移轨迹，其流程如图3所示。

图3 漂移模型建立流程

其中v表示物体漂移速度，vc和vw分别表示随流漂移速度和风致漂移速度。

则最终建立的物体漂移运动方程为：

(8)

2.2 漂流速度的求解

雷达输出的某一点的流场和风场数据的时间分辨率为10 min，对于这10 min内某一时刻流场和风场的数据需要通过插值的方式来进行精确地刻画。考虑到场变化的连续性且10 min内变化不会太大，可以用拉格朗日线性插值公式来进行插值。

(9)

vc(t1),vc(t2)为两个相邻采样时间雷达输出矢量流。

在实际计算过程中，线性插值因为无法反映出海上流场、风场变化的周期性，计算结果与实际估计值误差较大，因此采用傅里叶级数方法[20]进行数据插值。

由于高频地波雷达返回的是一系列经纬度点的风场和流场数据，这些经纬度点可将海平面划分为二维直角网格，人员的落水点就在这些由经纬度点确定的网格区域内。用(i,j)表示每一个数据点的坐标，网格的顶点即为雷达的数据采样点，该处任意时刻海流和风场的速度由傅里叶级数插值法给出。考虑到风速和流速的计算满足矢量加法，故统一将其分解为x和y方向，分别计算物体x和y方向上的位移，从而确定运动轨迹上一系列点的坐标。以海流为例，将海流的速度按x和y方向分解为vcx(i,j)和vcy(i,j)，记网格边界上海流的速度分别为u(i,j)和v(i,j)，如图4所示。

图4 网格边界速度分量

设每一个小矩形网格的x,y方向的宽度分别为：ΔX,ΔY。

那么对于以点(i,j)为左上顶点确定的网格中的落水点，设其坐标为(x0,y0)。定义其在网格中的坐标(iT,jT)为：

(10)

根据线性插值公式，网格边界速度为：

u(i,j)=(vcx(i+1,j)-vcx(i,j))·iT+vcx(i,j)

v(i,j)=(vcy(i,j+1)-vcy(i,j))·jT+vcy(i,j)

(11)

落水点的海流速度即为网格边界速度的线性插值：

Vx(iT,jT)=(1-jT)·u(i,j)+jT·u(i+1,j)

Vy(iT,jT)=(1-iT)·v(i,j)+iT·v(i,j+1)

(12)

其中Vx和Vy分别为落水点速度的x和y分量。

2.3 随机运动的求解

由于雷达返回的数据本身有误差，且漂流过程中物体所受湍流的作用不可知，因此物体的漂移运动势必会表现出随机性，这种随机运动遵循维纳过程。

蒙特卡罗方法的基本思想是将漂移人体定义为单个粒子，粒子的运动轨迹是一个满足一定概率分布的随机过程。大量复制这些粒子，则一段时间后粒子分布的统计区域就是最终确定搜救的概率区域。

同一粒子在时间步长Δt前后位置满足：

(13)

其中ξ(t)由Matlab生成大量满足正态分布的随机数构造。

蒙特卡罗方法运用统计学的相关思想，引进随机过程进行分析。系统具有良好的通用性。漂移模型建立参数相对简单，当需要进行大量的数据模拟时，可以方便的采用分布式运算将模拟过程放置于若干个微机系统独立运行，再将运算结果n(i,j)予以叠加。既可以满足统计的需要，同时也不必担心负荷过重，耗费过多机时，可以快速的取得计算结果。

2.4 系统总体流程图

高频地波雷达返回的实时海流海风数据，经过数据库的分类和处理后，输出可用的海风海流速度大小和方向，进入系统后台。用户在系统前端网页上输入落水人员落水时间，位置，性别，是否穿戴救生衣等相关信息，后台自动开始计算，模拟出漂移人员的运动轨迹，并在前端海图上显示可视的搜救区域。系统总体流程如图5所示。

图5 系统总体流程

3 数据仿真与结果分析

3.1 漂移轨迹的绘制

调用高频地波雷达数据库2013年2月1日6：35—2013年2月1日19：35东经118°0'北纬23°56′至东经118°6′北纬23°43′共13小时98个点的数据，每隔1′建立一个网格，构建二维经纬网络。

用户端输入落水时间8：30，落水位置东经118°2′北纬23°45′，则系统预测的落水人员10 h运动轨迹如图6所示。

图6 系统预测落水人员漂移轨迹

取3 h预测时长，以每30 min为单位，落水人员的位置和距落水点的漂移距离如表1所示。

表1 落水人员位置和距落水点距离

3.2 搜救区域的确定

用户端输入落水时间为7：00，落水位置为东经118°1′北纬23°54′24″，系统每隔30 min输出一次搜救区域，则4 h内可能的搜救区域如图7所示。

图7 4 h内搜救区域

图8 4 h后搜救区域放大图

从图中可以看出，搜救区域整体的偏移运动符合预测出的漂移轨迹，同时由于随机运动的存在，搜救区域的面积随时间增长在不断扩大。放大4小时的搜救区域如图8所示，粒子集中于东经118.06°至118.061°，北纬23.9265°至23.927°确定的矩形网格内，该网格东西向宽度约为101.58 m，南北向宽度约为55.66 m，实际搜救区域面积约为5653.80 m2。同时粒子分布的密集程度也表示了物体在该区域的概率大小。

4 结 语

本文从高频地波雷达应用的实例着手，探讨地波雷达在海上人员搜救方面应用的优越性。根据雷达返回数据的矢量流，通过拉格朗日追踪法建立动力学模型，并构建二维网格插值模型进行求解。在插值方法上，针对传统线性插值法的不足之处，使用傅里叶级数插值方法，极大地减小误差。在搜救区域的确定上，针对传统解析法的计算量大，迭代次数多，精度低的问题，选用蒙特卡罗方法，运用统计学规律，提高计算效率，从而较为精确地确定搜救区域，为搜救工作起到辅助决策的作用。但从仿真结果中发现，由于雷达返回的海风海流数据存在有±15cm/s的误差，用插值方法替代速度随时间连续变化会产生误差，以及风导因子的选择存在经验性和地区适应性，不同取值会影响最终结果。

高频地波雷达返回的矢量流数据在溢油的预测和赤潮的监测等问题中同样可以应用。对于本文的算法和模型，可以将其进一步推广，增加溢油扩散预测和赤潮监测功能。

[1] 马晓雪,石树凯,张成瑶,等.提升我国海上搜救公共服务能力之对策[J].水运管理.2016(05):12-15.

[2] 王辉,万莉颖,秦英豪,等. 中国全球业务化海洋学预报系统的发展和应用[J].地球科学进展.2016(10):1090-1104.

[3] 待茂崇,高郭平,鲍献文. 海洋调查方法导论[M].青岛：中国海洋大学出版社,2008:5-12.

[4] 林幼权, 星载海洋监视雷达系统[J].现代雷达,2012(11):6-10.

[5] 赵金赛,米伟,白树祥. 无人机在海上搜救中的应用探索[J].中国海事, 2014(08):42-44.

[6] 魏国妹,商少平,贺志刚,等. OSMAR-S100便携式高频地波雷达海浪和海面风探测性能分析[J].海洋与湖沼,2016(01):52-60.

[7] LEI R,NASH S, HARTNETT M. Observation and modeling of tide- and wind-induced surface currents in Galway Bay[J].Water Science and Engineering.2015(04):345-352.

[8] 吴雄斌,张兰,柳剑飞. 海洋雷达探测技术综述[J]. 海洋技术学报, 2015(03):8-15.

[9] 潘军.基于 Monte Carlo 方法的确定海上搜救区域模型[J].中国水运,2012(09):43-44.

[10] 李伦. 高频地波雷达海洋动力学参数反演与应用方法研究[D]. 武汉大学,2013:86-88.

[11] 王道生,刘晓东,庄海东.一种改进的海上搜救区域预测方法及系统:中国,CN 105653826 A[P].2016-06-08.

[12] 牟林,宋军,李欢,等.一种海上搜救目标漂移路径的预测方法:中国,CN 103366227 A[P].2013-10-23.

[13] 刘建,沈兵,冯伟,等.一种搜救终端和基于该搜救终端的水上搜救方法:中国,CN 102956089 A[P].2013-03-06.

[14] 姜华林,孙昭晨,李黎,等. 基于蒙特卡洛方法的海上搜寻区域确定模型[J]. 水道港口, 2011(04):285-290.

[15] ALLEN A,PLOURDEJ.Review of leeway:field experiments and implementation[R].USA:US Coast Guard Research and Development Center.1999:8-12.

[16] 苗星妍. 海上运动点连续溢油预测及其三维可视化研究[D]. 集美大学,2015：12-16.

[17] DANIEL P,MARTYF,JOSSE P.Improvement of existing operational oil spill and object drift prediction system in western Mediterranean sea[J]. Proceedings of Enroute to GODAE . 2002:329-330.

[18] 肖方兵,尹勇,金一丞,等. 基于随机粒子仿真的海上搜寻区域确定[J]. 中国航海, 2011(03):34-39.

[19] 翁怡婵,杨金湘,江毓武. 台湾海峡漂移物运动轨迹的数值模拟[J]. 厦门大学学报(自然科学版). 2009(03):446-449.

[20] 王向磊,苏牡丹,丁硕,等. 区域地磁场的连续傅里叶分析建模方法[J]. 测绘科学技术学报, 2011(06):425-428.