水下“寻宝”猎人：落水物品快速定位与打捞

Underwater “Treasure” Hunter: Quickly Locate and Retrieve Items That Have Fallen into the Water

摘要: 水上游玩时，常容易出现贵重物品掉落(如手机、身份证等)的突发情况，如何高效打捞成为了一大难题。本文对掉落物品物理性质、掉落流域特性的相关数据进行分析，建立了基于纳维–斯托克斯方程的手机落水位置分析模型，针对京杭大运河拱宸桥段流动水域特点，分析了手机落水后的受力情况，包括流水推力、阻力、重力和浮力。通过建立笛卡尔坐标系，运用N-S方程和库埃特流动理论，预测了手机在水平和竖直方向上的运动轨迹。结合水平位置预测模型和深度预测模型，构建了完整的时空预测模型，并通过模拟测试验证了模型的通用性。针对轻质量物品分析了居民身份证物理特性，并对前述模型进行调整，考虑了风速对身份证的影响，优化了流体力学模型和纳维–斯托克斯方程的位置分析模型。并基于遗传算法，分析了最优搜索路径并生成了可视化结果。

Abstract: When playing on water, it is common to encounter unexpected situations where valuable items (such as mobile phones, ID cards, etc.) fall off, making efficient salvage a major challenge. This article analyzes the relevant data on the physical properties and watershed characteristics of falling objects, establishes a model for analyzing the location of mobile phone falling into the water based on the Navier-Stokes equation, and analyzes the force situation of mobile phone falling into the water, including water flow thrust, resistance, gravity, and buoyancy, according to the characteristics of the flowing water area in the Gongchen Bridge section of the Beijing Hangzhou Grand Canal. By establishing a Cartesian coordinate system and applying N-S equations and Couette’s flow theory, the motion trajectory of the mobile phone in both horizontal and vertical directions was predicted. A complete spatiotemporal prediction model was constructed by combining a horizontal position prediction model and a depth prediction model, and the universality of the model was verified through simulation testing. We analyzed the physical characteristics of resident identity cards for lightweight items and adjusted the aforementioned model, taking into account the influence of wind speed on identity cards. We optimized the fluid dynamics model and the Navier-Stokes equation position analysis model. Based on genetic algorithm, we analyzed the optimal search path and generated visual results.

1. 符号说明及数据来源

为更好地解决实际问题、建立数学模型，我们将现实场景中所需运用到的数据进行符号表示，全文中所需的变量、常量如(表1)所示：

Table 1. Model symbol explanation

表1. 模型符号说明

符号	含义	单位
$M_1$	华为Mate60pro的质量	kg
$V_1$	华为Mate60pro的体积	m^3
$S_1$	华为Mate60pro较大面的面积	m^3
$M_2$	身份证的质量	kg
$V_2$	身份证的体积	m^3
$S_2$	身份证较大面的面积	m^3
$h$	物品掉落位置与水面的距离	m
$v_0$	游船行驶速度	km/h
$v_1$	华为Mate60pro/身份证落入水面的初速度	km/h
$v_t$	华为Mate60pro/身份证经过t时间后的速度	km/h
$v_r$	京杭大运河拱宸桥的河流流速	km/h

为更好体现水体变化中的时时性、数据动态分析，我们于earth.nullschool.net中找到了地表风速及洋流流速变化的数据，便于进行模拟检验。同时我们从华为官网中收集到了Mate60pro规格(如表2)所示。

Table 2. Data sources

表2. 数据来源

符号	来源
华为手机Mate60pro规格	https://consumer.huawei.com/cn/phones/mate60-pro/specs/
风速及洋流流速	https://earth.nullschool.net/#current/wind/surface/level/orthographic=-246.80,29.70,887/loc=111.492,23.143

2. 基本假设

• 假设收集的数据均为真实测量数值，不存在记录失误等原因造成的数据错误。

• 不考虑西湖和京杭大运河中除水流以外的其他因素(如水中生物、其他船只经过影响的水流改变等)对手机以及居民身份证的位置影响。

• 假设西湖和京杭大运河拱宸桥附近的天气稳定，不考虑暴雨强风等突发变化。

• 引入水流随时间变化的函数。收集拱宸桥段不同位置和不同时刻的水流速度数据，利用插值法或拟合函数来描述其变化规律，比如通过多项式拟合得到水流速度$v\left(x,y,t\right)$ 关于空间坐标$\left(x,y\right)$ 和时间t的函数表达式，然后将其代入纳维–斯托克斯方程和其他相关的模型方程中，从而使模型能反映实际水流速度的动态变化。

3. 基于纳维–斯托克斯方程的位置分析

3.1. 受力分析与建立笛卡尔坐标系

Figure 1. Force analysis of mobile phones in the Grand Canal between Beijing and Hangzhou

图1. 京杭大运河中手机的受力分析

拱宸桥段于京杭大运河最南端，地处杭嘉湖平原，为长江三角洲的一部分。假设在京杭大运河拱宸桥附近，从游船上掉落到水里一款华为Mate60pro手机。与西湖不同的是，京杭大运河是流动水域，且在江南段的流向为自北向南。因此在手机掉落至运河中，流水对手机会产生推力F，同时为了阻碍手机的相对运动，手机会受到运河水的阻力f。在竖直方向上，手机受到重力与浮力的作用。经推算，$G\approx 2.25\text{N}$ ，$F_{浮}\approx 1.05\text{N}$ 。以手机抛出点为原点，建立笛卡尔坐标系，以水流方向为x轴，垂直水流方向为y轴，z轴竖直向下。得到受力分析空间图如图1。

3.2. 水平位置预测模型

3.2.1. 基于纳维–斯托克斯方程的物体运动轨迹分析

基于问题一的前期建设，我们可以得到手机落水的初速度为$v_1$ 。基于笛卡尔坐标系，以手机抛出点为圆心，对运河水空间进行网格化分割。根据现场监测，可得水流速度为$v_r$ (外界因素对于水流速度的影响较小，可忽略不计)，我们将继续模拟流体动力对其轨迹的影响。在水平面上，当手机移动到坐标点$\left(x,y\right)$ 时，我们有矢量表达式：$\left(\text{d}x,\text{d}y\right)=v\text{d}t$ 。

浮力与重力只与手机竖直方向上的运动有关，与水平点坐标$\left(x,y\right)$ 无关。由此根据受力分析图，我们可知只有阻力f与运河水流动产生的推力F决定水平运动。根据牛顿第二定律，得矢量表达式：

$F+f=Ma=M\frac{\text{d}v}{\text{d}t}=M\frac{{\text{d}}^2\left(x,y\right)}{\text{d}{t}^2};{\text{d}}^2\left(x,y\right)=\frac{F+f}{M}\text{d}{t}^2$

F表示水流冲击力，其可用冲量与动量的微分形式表示：$F=\rho Q{v}^{2}S$ 。

纳维–斯托克斯方程(Navier-Stokes equation)是描述粘性不可压缩流体动量守恒的运动方程，简称N-S方程。手机在$\left(x,y\right)$ 处的速度分量为$\left(u,w\right)$ ；因此可以用来描述运河流水对手机的阻力效应：

$\{\begin{array}{l}\nabla v=0\\ f=\rho \cdot \left(\frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)\cdot v\right)+\nabla p-\mu {\nabla }^{2}v\end{array}$

ρ为河水的密度，t为时间，v为手机相对于运河水的速度矢量，P为手机在水中的压力，µ为动态粘度系数，f₀为阻力，$\nabla$ 是偏微分算子。

${\text{d}}^2\left(x,y\right)=\frac{F+f}{M}\text{d}{t}^2=\frac{\rho Q{v}^{2}S+\rho \cdot \left(\frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)\cdot v\right)+\nabla p-\mu {\nabla }^{2}v}{M}\text{d}{t}^2$

$\{\begin{array}{l}{\text{d}}^{2}x=\frac{\rho Q{u}^{2}S+\rho \cdot \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+w\frac{\partial u}{\partial y}\right)+\frac{\partial p}{\partial x}-\mu \left(u\frac{{\partial }^{2}u}{\partial x^2}+w\frac{{\partial }^{2}u}{\partial y^2}\right)}{M}\text{d}{t}^2\\ {\text{d}}^{2}y=\frac{\rho Q{w}^{2}S+\rho \cdot \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+w\frac{\partial w}{\partial y}\right)+\frac{\partial p}{\partial y}-\mu \left(u\frac{{\partial }^{2}w}{\partial x^2}+w\frac{{\partial }^{2}w}{\partial y^2}\right)}{M}\text{d}{t}^2\end{array}$

它是一个非线性偏微分方程，求解非常困难和复杂，在求解思路或技术没有进一步发展和突破前只有在某些十分简单的特例流动问题上才能求得其精确解；但在部分情况下，可以简化方程而得到近似解。

3.2.2. 基于库埃特流动的偏微分方程求解

N-S方程在偏微分方程中的研究是出了名的难以得到精准的解值，因此我们考虑到了拱宸桥水体的空间特征，对偏微分方程的项进行弱化。基于前期的文献检索，我们在刘益，付小莉关于《Matlab求解库埃特流在计算流体力学课程中的应用》中得到了启发，计算流体力学主要研究如何通过数值方法来求解流体力学的控制方程。本文以不可压缩库埃特流为例,在有理论解析解的基础上建立显式、隐式方法的有限差分法控制方程，并提出了基于压力修正法的第三种方法，采用Matlab编程求解流场分布[1]。希望利用这个技术对方程的解进行估计。

可求得精确解的最简单情况是平行流动，而京杭大运河在江南段自北向南流动恰好可满足平行流动的近似解。查得资料，拱宸桥段中心处最深处为9.8 m，两岸侧边平均水深为5.4 m，而拱宸桥流段，跨度达101.2 m。设岸边最浅点至岸中最深点连线与水平面成$\beta$ 角，$\beta =\arctan \left(4.4/50.6\right)\approx 4.97˚$ ，说明运河底部曲面曲率几乎为0，可近似为平面。于是可将运河截面视为矩形，运河载体为无限延伸的空间长方体(图2)。

Figure 2. Canal section of Gongchen Bridge section

图2. 拱宸桥段运河截面

Figure 3. Motion trajectory of xOy plane (horizontal plane)

图3. xOy平面(水平面)运动轨迹

库埃特流动是指粘性流体在相对运动着的两平行平板之间的层流流动(图3)。在流体动力学中，库埃特流动是两个表面之间的空间中粘性流体的流动，其中一个流体相对于另一个表面正切向移动。

下面考虑具有深度为h与运河长度l的无限长的矩形通道，经受顶壁以恒定速度$v_r$ 移动的条件。无任何施加的压力梯度，Navier-Stokes方程减少到：

$\frac{{\partial }^2{v}_r}{\partial x^2}+\frac{{\partial }^2{v}_r}{\partial y^2}=0$

算法优化：手机掉落的动水搜索模型

① 测量并输入$v_r、h、v_1、\mu$ ，记录手机抛出点的经纬度；

② 由Navier-Stokes方程计算出手机在水中受到的阻力；

③ 微分方程处理：由

${\text{d}}^2\left(x,y\right)=\frac{F+f}{M}\text{d}{t}^2=\frac{\rho Q{v}^{2}S+\rho \cdot \left(\frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)\cdot v\right)+\nabla p-\mu {\nabla }^{2}v}{M}\text{d}{t}^2$

可得该方程为二阶非线性偏微分方程。可根据流体的不可压缩性与化简Navier-Stokes方程为：

$\frac{{\partial }^2{v}_r}{\partial x^2}+\frac{{\partial }^2{v}_r}{\partial y^2}=0$

④ 算法的迭代求解：设置步长h = 1 s，迭代10个小时。标出在xOy平面上的离散点列。

3.3. 深度预测模型

竖直方向上的运动分析

Figure 4. Motion trajectory of xOz plane (vertical direction)

图4. xOz平面(竖直方向)运动轨迹

在竖直方向上仅有重力、浮力与阻力在z轴上的分力$f_z$ 会对手机竖直方向上的速度构成影响。显然在长时间的运动下，阻力与水平面成的夹角会逐渐变小，由于$G>F_{浮}$ ，手机若没有及时打捞，最终会为沉底状态(图4)。

在竖直方向上的运动轨迹，与问题一相同，不再赘述。

$\text{d}z=h+\frac{\left(G-\rho Vg\right)t^2-2kt\sin \alpha }{2M}$

3.4. 航迹线搜索方式

3.4.1. 时空模型预测结果

综合上述水平位置预测模型与深度预测模型，将两个模型在三维空间中叠加张开，我们得到了完整的时空预测模型。输入必要参数后，即可获得手机在京杭大运河落水后的空间轨迹图。

接下来，为了测试模型的通用性，我们选择了初始环境进行模拟。取定京杭大运河平均流速为$v_r=0.5\text{km}/\text{h}$ ，手机抛出经纬坐标为(120˚19'02'', 30˚08'06'')，$h=0.5\text{m}$ ，$v_1=v\left(3,4,5\right)$ ，$\mu =2.98\times {10}^{-3}\text{Pa}\cdot \text{s}$ 。

由于手机掉落至水中，无法准确获得潜水艇的实时位置。面对这种情况，最佳采用的搜索方式为航迹线搜寻。由于我们已经得出潜水艇丧失动力t秒后的坐标为$\left(\text{d}x,\text{d}y,\text{d}z\right)$ ，当Δt取得很小，连续取n个点(n尽可能大)，可视为潜水艇失事的航迹线。在此我们取Δt = 1 s，n = 10,000，画出一条航迹线如下图5所示。

Figure 5. The trajectory of a dropped mobile phone

图5. 掉落手机航迹线

3.4.2. 基于蒙特卡洛概率模型的航迹线搜索方式

运河底部仍然拥有较多我们不能确定的不稳定因素，所以计算出的表达式不能精准的得出手机的“航迹线”，而是以计算出的点为中心，更新手机位置的概率分布。根据每秒接受到的数据，考虑扰动后，估计手机的位置概率。

在这里我们引入蒙特卡洛概率模型。孙少华，尚雅玲，李国林[2]在《基于蒙特卡洛法的制导子弹药命中概率模型》中提出运用蒙特卡洛方法建立了制导子弹药命中概率计算模型，并对不同捕获概率下对给定目标的命中概率进行了分析。因此我们将其模型迁移运用于手机掉落的位置概率点。在流体力学的应用方面，蒙特卡洛方法可以用于模拟流体的输运和扩散过程。流体的扩散过程近似于潜水艇受外界干扰，而偏离指定航迹。通过随机抽样和迭代计算，我们可以模拟手机的掉落过程，并分析流体中的输运性质和混合效应。

运用蒙特卡洛算法计算流体扩散到边界的概率，我们首先要选择一个适当的边界形状和位置。结合实际，我们选定t时刻手机所处位置M为中心点，扩散概率的边界为以M为中心的半径为r的同心圆。

选定完边界之后我们开始模拟扩散情况，设置Δt = 1 s，针对每一个时间步长进行布朗运动模拟，计算潜水艇位移并不断更新位置。在单位时间内检测潜水艇(拟为流体粒子)是否与边界发生碰撞。若碰撞到相应边界，则标记为偏移位置的一次可能性。

进行n次模拟后，对模拟结果进行统计分析，记录每次模拟扩散到对应边界的粒子数目，最终得到计算边界的概率结果。将概率结果可视化后我们得到概率密度图，为方便观测，我们给画出单点的概率密度切面形式如下(图6)：

Figure 6. Probability density section of a single point

图6. 单点的概率密度切面

根据蒙特卡洛算法我们可以得到找到潜水艇的概率。并通过matlab软件，模拟手机掉落5 h与10 h的手机于xOy平面可能掉落的位置，进行对比预测(图7，图8)。

Figure 7. 30 minute probability prediction of mobile phone position xOy plane

图7. 手机位置xOy平面30 min概率预测

Figure 8. 1-hour probability prediction of mobile phone position xOy plane

图8. 手机位置xOy平面1 h概率预测

结果表明，手机刚刚落水30分钟内，手机掉落位置概率稿的密集范围基本呈条带状；落水时间越长，手机掉落的高概率位置分布越广，概率图的无序度扩大，规律性消失，呈散点状搜寻难度大大增加。

4. 基于时空预测模型下的身份证掉落

4.1. 运动分析

身份证的质量较手机相比，非常的轻，所以在抛出阶段更容易收到风速的影响，对入水速度$v_1$ 有较大影响，从而对运动轨迹产生影响。由2.2.1已对增加风力影响后的入水速度${v^{\prime }}_1$ 重新计算，身份证入水后近似处于液体密封状态。且$F_{浮}\approx 0.046\text{N}$ ，$G\approx 0.066\text{N}$ ，因此身份证应该在水中不断移动，而非在运河表面漂浮。此外，若为轻质物品，则在水面流动会较大程度上受到风速影响，这主要源于风对其施加的

压力和摩擦力。依据空气动力学原理，风压力$F_{\text{wind-pressure}}$ 可通过公式$F_{\text{wind-pressure}}=\frac{1}{2}{C}_p{\rho }_{\omega }{v}_{\omega }^{2}A$ 计算，其中$C_p$

是压力系数，${\rho }_{\omega }$ 为空气密度，A是轻质物品迎风面积。风摩擦力$F_{\text{wind-friction}}={\mu }_{\omega }{v}_{\omega }A$ ，${\mu }_{\omega }$ 为空气与物品间的摩擦系数。我们需将这些力纳入物体在水面的受力分析体系，与原有的重力、浮力、水流作用力等综合考虑。由于分析知，身份证会在水中不断移动而不是在运河表面漂浮，因此我们暂不考虑这种情况。

${\text{d}}^2\left(x,y\right)=\frac{F+f}{M}\text{d}{t}^2=\frac{\rho Q{v}^{2}S+\rho \cdot \left(\frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)\cdot v\right)+\nabla p-\mu {\nabla }^{2}v}{M}\text{d}{t}^2$

代入代码后得到位置图为(图9)：

Figure 9. Location of ID card entering the water

图9. 身份证入水位置

4.2. 基于遗传算法的最优搜索路径分析

以身份证掉落一小时为背景，对每个概率高的点取点做遗传算法分析，找出最优搜索路径。本研究将种群规模设定为200，由于种群规模决定了每一代中包含的个体数量，较大的种群规模能提供更丰富的解空间多样性，使算法有更多机会搜索到全局最优解。本文参照文献[3]，将其应用于落水搜索路径中。

在TSP问题里，更多的个体意味着更多不同的路径组合被探索，降低陷入局部最优的可能性。但种群规模过大也会增加计算量和运行时间。此处设置200，是在计算效率和搜索全面性之间权衡的结果，既能保证一定的多样性，又不会使计算负担过重；将最大遗传代数设置为2000。由于遗传代数决定了算法的迭代次数，随着代数增加，算法不断进化种群，逐步优化解的质量。

较大的遗传代数能让算法有更充分的时间搜索最优解，但可能导致过度计算。设置为2000代，旨在给算法足够的迭代次数来收敛到较优解，同时避免不必要的长时间运算；在代码中，我们设置

$W=\text{ceil}\left(N/10\right)$

用其来确定交叉点个数，根据城市数目N动态确定交叉点个数，我们的这种设置使交叉操作的变化程度与问题规模相关。

例如，城市数量较多时，适当增加交叉点个数，能更充分地交换基因，促进种群的多样性和进化，城市数量较少时，较少的交叉点个数可避免过度打乱优良基因组合；我们再设

$p=\text{unidrnd}\left(N-W+1\right)$

随机选择交叉范围的起始点，从p到p + W的范围决定了父代个体中哪些基因片段会进行交换。随机选择交叉范围增加了算法的随机性和探索能力，使算法能在不同的基因组合间进行尝试，有助于发现更优的路径；最后在代码中，我们采用变异操作，通过随机选择两个位置

$p_1=\text{floor}\left(1+N\ast \text{rand}\left(\right)\right)$

$p_2=\text{floor}\left(1+N\ast \text{rand}\left(\right)\right)$

Figure 10. Genetic algorithm path

图10. 遗传算法路径

并将两个位置进行交换，来生成随机位置。由此，我们得到以下数据图(图10，图11)。

Figure 11. Fitness evolution curve

图11. 适应度进化曲线

5. 模型的评价与推广

5.1. 模型的优点

① 本文结合了纳维–斯托克斯方程来模拟水流对手机运动的影响，同时考虑流体的惯性力、压力梯度力、粘性力和外部作用力，提供了一个全面的流体动力学描述。

② 搜索模型结合流体力学模型与身份证的物体特性，使用了遗传算法寻找搜索方法，原理直观简洁富有灵活性，且编码容易实现。

5.2. 模型的缺点

① 流体力学模型基于理想化的假设，如不可压缩性、牛顿流体行为等，这些假设不一定适用于所有实际情况。

② N-S方程是非线性的，求解难度较大。特别是拱宸桥为活水，流动条件较为复杂，精确解较为难以计算。

5.3. 模型的推广

针对水中掉落物品，为提升成功打捞的效率，我们制作了一份搜索路径图(图12)。

Figure 12. Search path method

图12. 搜索路径方法