摘要
本文针对多波束测线问题，建立了二维几何模型、三维解析几何模型、单目标优化模型和多目标优化模型，分别解决了特定位置下的指标计算、覆盖宽度求解以及最优测线规划等问题。
一、问题重述
1.1 问题背景
由于陆地资源日益减少，人口规模迅速增长，经济发展的紧迫需求，人们迫切需要开采、开发和利用海洋资源。海洋已成为21世纪的现实和战略关注点，并承载着更重要的意义。而从事海岸工程建设、地质调查与开发等海洋活动，均需要海洋测绘勘探工作。海底地形测量作为这项工作的基石，需要综合运用多学科的理论、方法和技术。为了更加深入地了解海洋环境并精确掌握海底地形地貌，通常会采用多波束测深技术。由于其全覆盖、全水深、高精度和高分辨率等优势，多波束测深系统已成为当前主流的测深设备，在水下地形测量、海底地貌勘察和水下目标检测等领域得到广泛应用。
相较于单波束测深，多波束测深系统在垂直范围沿航迹上能够快速而准确地测量水下目标和地形起伏的变化。它将传统的点、线探测扩展至面的层次，实现了数字化记录和立体化自动测图。通过克服单波束测深的采样间距过大导致海底信息反映粗糙以及发射波束角度过大导致微小地形变化的深度误差等缺陷。因此，可以更准确、更广泛地获取航道的相关数据，从而全面提升航道测绘技术的精确度。
1.2 问题重述
多波束测深系统中，多波束测深条带的覆盖宽度W会随着换能器的开角θ和水深D的变化而发生变化。本文的数学模型主要解决以下问题：
- 在海底不平的情况下，当测量船测线方向垂直的平面和海底坡面的夹角为α时，建立多波束测深的覆盖宽度及相邻条带之间重叠率的数学模型。当开角θ=120°，坡度α=5°，海水深度为70米时，根据建立的数学模型求解相关参数，并将结果填入表1。
- 在海底不平的情况下，考虑一个矩形待测海域，当测量船测线方向与海底坡面的法向在水平面上投影的夹角为β时，建立多波束测深覆盖宽度的数学模型。当换能器的开角为θ=120°，坡度为α=1.5°，海域中心点处的海水深度为120 m时，根据覆盖宽度模型计算不同测量方向夹角的覆盖宽度，并将结果填入表2。
- 在海底不平的情况下，在一个南北长2海里、东西宽4海里的矩形海域内，海域中心点处的海水深度为110 m，呈西深东浅的坡度α=1.5°。多波束换能器的开角为θ=120°。现需要设计一组测量线，使其长度最短且可以完全覆盖整个海域，同时满足相邻条带之间重叠率在10% ~ 20%之间。
- 在海底不平的情况下，根据已知的海水深度数据，要求设计出满足如下要求的测线：沿测线扫描形成的条带尽可能地覆盖整个待测海域；相邻条带之间的重叠率尽量控制在20%以下；测线的总长度尽可能短。以该设计测线，要求计算出测线的总长度；漏测海区占总待测海域面积的百分比；在重叠区域中，重叠率超过20%部分的总长度。
二、问题分析
问题整体思路较为清晰，层层递进。问题一要求测量间距为200米是关于多波束测深覆盖宽度和相邻条带覆盖率参数的计算问题。问题二建立在问题一相关参数求解的基础上，增加一个矩形待测海域，测量间距由200米变为0.3海里，是关于在不同测线方向夹角下，覆盖宽度参数的计算问题。问题三是建立在问题二矩形待测海域的基础上，设置矩形海域宽度，是关于如何优化待测海域测线的问题。问题四是建立在问题三优化模型的基础上，是关于如何优化设计出满足要求的测线和相关指标的计算问题。
2.1 问题一的分析
问题一要求在海水起伏变化大，每隔200m测量，测量船测线方向垂直的平面和海底坡面的夹角为α的情况下，建立多波束测深覆盖宽度及相邻条带之间重叠率的数学模型，以求解特定位置下的指标。从题目内容可知，要求计算出测量船波束间的覆盖宽度W，即计算测量船波束边界点与斜坡交点之间的距离，需要得知波束与斜坡的交点坐标，再根据交点坐标实现重叠率的计算。因此，考虑使用二维几何模型来解决问题。假设海面近似为平面，以海域中心为原点，以航向垂直的方向为x轴，铅直方向为y轴建立直角坐标系，并列出波束边界点与测量点的坐标公式；其次，根据上述三点坐标，列出斜坡直线方程和给定位置处的波束边界方程，构建多束波覆盖宽度的数学模型。求解方程的一般计算方法有高斯消元法，二分法，迭代法等算法，本题宜于选用二分法求解。计算出波束与斜坡的交点坐标，以此实现测线重叠率的计算。最终求解出在不同测线距中心点处距离下的海水深度，覆盖宽度，重叠率。最后，分析了重叠率对斜坡角度、换能器开角和海水深度的灵敏度。利用黄金分割法对上述结果进行了重新求解，计算结果误差均较小，表明计算结果的可信度较高。
2.2 问题二的分析
问题二在问题一的基础上，测量间距变为0.3海里，要求考虑一个矩形待测海域，即测线方向与海底坡面的法向在水平面上投影夹角为β的平面的情况下，建立多波束测深覆盖宽度的数学模型以求解特定位置下覆盖的宽度。从题目内容可知，要求覆盖的宽度，求出波束直线与斜坡平面的交点坐标即可。因此，考虑使用三维几何模型来解决问题。首先，以海域中心为原点，以平行于海底交线方向为x轴，垂直交线方向为y轴，铅直方向为z轴建立直角坐标系，并列出波束边界点与测量点的坐标公式；其次，求出海底斜坡平面方程和给定位置波束边界的直线方程。采用牛顿迭代法求解波束直线与斜坡平面的交点坐标，再由交点坐标实现覆盖宽度的计算，最终求出不同测线方向夹角下，不同测量间距的覆盖宽度。最后，分析了覆盖宽度对测线间距、测线方向角和海水深度的灵敏度，验证了模型的合理性。
2.3 问题三的分析
问题三在问题二矩形待测海域的基础上，设置了具体的长宽矩形海域，要求以使测线累计距离最小为目标，以特定的海域中心处的海水深度，坡度，多波束换能器的开角，相邻条带重叠率满足10%到20%为约束条件，建立优化模型。从题目内容可知，该题属于测线长度的优化问题，考虑使用单目标优化模型来解决问题。首先，在问题二直角坐标系的基础上，给出平行测线的直线方程，并求出测线累计距离，以该距离最小为目标；其次，求出各测线与边界的交点，并求得这些交点处的重叠率，以重叠率满足10% ~ 20%为约束条件，建立单目标优化模型。解决此类模型的一般计算方法有梯度下降法，粒子群优化，差分进化算法等方法，本题采用差分进化算法进行求解，计算得出最优的测线总长度，此时最大最小的重叠率分别对应的海水深度。最后，利用追赶法对模型的求解结果进行检验，计算结果与差分进化算法的求解结果误差较小，验证了计算结果的可信度性。
2.4 问题四的分析
问题四在问题三矩形待测海域的基础上，给定另一不同矩形海域的海水深度数据，要求对该数据进行拟合处理，以此结果设计出距离最短，重叠率大于20%的比例最小的测线，给出最优的多波束测线方案。从题目内容可知，该题属于多波束测线方案的优化问题，考虑使用多目标优化模型进行求解。首先，利用二元三次幂函数对附件中的海深数据进行拟合，得到拟合的曲面方程；其次，利用航向角求出累计测线距离，并以该距离最小为目标1。由航线求出航线在拟合曲面上的投影方程，并利用波束边界直线方程确定覆盖宽度范围，进而计算重叠率，以重叠率大于20%的比例最小为目标2，建立多目标优化模型。在等权条件下，将多目标转为单目标优化模型。采用混合模拟退火禁忌搜索算法对优化模型进行求解，得出最优的测线总长度，漏测海区占总海域的面积比等参数，最后，采用灰狼优化算法对模型结果进行重新检验，误差较小，验证了计算结果的可靠性。
三、基本假设
- 假设海平面近似平面；
- 假设不考虑测量船的姿态；
- 假设海面平静；
- 假设海水折射率均匀分布。
四、符号说明
α：坡度
θ：换能器的开角
Pi：测线与海底坡面的左交点
Qi：测线与海底坡面的右交点
M0：沿海底坡面方向的直线
Mi：沿测线方向的第一条直线
Mi+1：沿测线方向的第二条直线
D：海水深度
a：测线距离中心点处的距离
五、模型建立与求解
5.1 问题一的模型建立与求解
本题为解决特定位置下的指标计算问题，建立二维几何模型。以海域中心为原点，航向方向为x轴，铅直方向为y轴，建立平面直角坐标系。其次，根据坐标系示意图，求出斜坡直线方程和特定位置处波束边界方程，利用二分法求解出测量船的波束边界和斜坡的交点坐标，以此计算出重叠率和条带的覆盖宽度。
5.1.1 模型建立
5.1.1.1 二维几何模型的确立
本题旨在建立二维几何模型，求解特定位置下的指标。因此，以海域中心为原点，以航向垂直的方向为x轴，铅直方向为y轴，建立如下图的直角坐标系：
图 1 测量船波束与斜坡示意图
由上图可知，在xoy平面上，测量船位于海平面上，设测量船i所在的坐标为(xi, 0)。坡地倾角为α，沿海底坡面方向的直线M0的斜率为tanα，故沿海底坡面方向的直线M0的表达式为：
M
0
:
y
=
tan
α
⋅
x
−
D
M_0: y = \tan\alpha \cdot x - D
M0:y=tanα⋅x−D
其中D=70m为海水深度。利用i点的横坐标代入该方程即可得到海水深度。
由于平面xOy垂直于测线方向，且垂线平分多波束换能器开角θ，故沿测线方向的第一条直线Mi与平面的夹角为π - θ，该直线的斜率为tan(π - θ)，故沿测线方向的第一条直线Mi的表达式为：
M
i
:
y
=
tan
(
π
−
θ
)
(
x
−
x
i
)
M_i: y = \tan(\pi - \theta)(x - x_i)
Mi:y=tan(π−θ)(x−xi)
同理，沿测线方向的第二条直线Mi+1与平面的夹角为π - θ + θ = π + θ，斜率为tan(π + θ)，故沿测线方向的第二条直线Mi+1的表达式为：
M
i
+
1
:
y
=
tan
(
π
+
θ
)
(
x
−
x
i
)
M_{i+1} : y = \tan(\pi + \theta)(x - x_i)
Mi+1:y=tan(π+θ)(x−xi)
由沿海底坡面方向的直线方程M0，与沿测线方向的第一条直线方程Mi和沿测线方向的第二条直线方程Mi+1联立，可分别求出两条测线与海底坡面的交点Pi (xpi, ypi)，Qi (xqi, yqi)的坐标，即：
Pi (xpi, ypi)：
{
y
=
tan
α
⋅
x
−
D
y
=
tan
(
π
−
θ
)
(
x
−
x
i
)
\begin{cases} y = \tan\alpha \cdot x - D \\ y = \tan(\pi - \theta)(x - x_i) \end{cases}
{y=tanα⋅x−Dy=tan(π−θ)(x−xi)
Qi (xqi, yqi)：
{
y
=
tan
α
⋅
x
−
D
y
=
tan
(
π
+
θ
)
(
x
−
x
i
)
\begin{cases} y = \tan\alpha \cdot x - D \\ y = \tan(\pi + \theta)(x - x_i) \end{cases}
{y=tanα⋅x−Dy=tan(π+θ)(x−xi)
利用交点Pi, Qi的坐标，由欧式距离计算出PiQi的长度，即为条带的覆盖宽度W，即
W
=
∣
P
i
Q
i
∣
=
(
x
p
i
−
x
q
i
)
2
+
(
y
p
i
−
y
q
i
)
2
W = |P_iQ_i| = \sqrt{(x_{p_i} - x_{q_i})^2 + (y_{p_i} - y_{q_i})^2}
W=∣PiQi∣=(xpi−xqi)2+(ypi−yqi)2
对于i + 1节点，可利用上述分析得到Pi+1, Qi+1节点的坐标，利用重叠关系可得重叠率表达式如下。根据题意可知，相邻条带之间的重叠率定义为
η
=
1
−
d
W
=
∣
P
i
Q
i
∣
−
∣
P
i
+
1
Q
i
∣
∣
P
i
Q
i
∣
\eta = 1 - \frac{d}{W} = \frac{|P_iQ_i| - |P_{i+1}Q_i|}{|P_iQ_i|}
η=1−Wd=∣PiQi∣∣PiQi∣−∣Pi+1Qi∣
5.1.1.2 二维几何模型表达式
综上所述，最终的二维模型表达式如下：
{ W = ∣ P i Q i ∣ = ( x p i − x q i ) 2 + ( y p i − y q i ) 2 η = 1 − d W = ∣ P i Q i ∣ − ∣ P i + 1 Q i ∣ ∣ P i Q i ∣ y = tan ( π − θ ) ( x − x i ) y = tan ( π + θ ) ( x − x i ) x i ∈ { − 800 , − 600 , ⋯ , 600 , 800 } \begin{cases} W = |P_iQ_i| = \sqrt{(x_{p_i} - x_{q_i})^2 + (y_{p_i} - y_{q_i})^2} \\ \eta = 1 - \frac{d}{W} = \frac{|P_iQ_i| - |P_{i+1}Q_i|}{|P_iQ_i|} \\ y = \tan\left(\pi - \theta\right) (x - x_i) \\ y = \tan\left(\pi + \theta\right) (x - x_i) \\ x_i \in \{-800, -600, \cdots, 600, 800\} \end{cases} ⎩ ⎨ ⎧W=∣PiQi∣=(xpi−xqi)2+(ypi−yqi)2η=1−Wd=∣PiQi∣∣PiQi∣−∣Pi+1Qi∣y=tan(π−θ)(x−xi)y=tan(π+θ)(x−xi)xi∈{−800,−600,⋯,600,800}
其中， D = 70 m D = 70\ \text{m} D=70 m， α = 1.5 ∘ \alpha = 1.5^{\circ} α=1.5∘， θ = 120 ∘ \theta = 120^{\circ} θ=120∘
5.1.2 模型求解
5.1.2.1 模型求解的原理
根据斜坡方程和测量船波束边界方程，求出两个的交点坐标，根据交点坐标计算出覆盖宽度和重叠率。本题采用二分法来求解交点坐标。
二分法又称二分搜索法、折半搜索法或对数搜索法，是在一组顺序排列的元素中查找某一特定元素的搜索方法。它通过不断对分数据区间，将不可能的解集排除，得到满足误差条件的近似解。在每次二分迭代时都会将范围缩小一半，迭代次数少，查找速度快，性能较高。
5.1.2.2 模型求解的过程
步骤一：先确认定义域的区间[a, b]，输入端点及精度值确定区间[-800, 800]，验证f(-800) · f(800) < 0，如果不成立则一般不存在根或是区间太过宽泛，再给定精度e=0.001；
步骤二：求区间(-800, 800)的中点X0 = (a + b)/2，即X0 = 0；
步骤三：计算f(X0)，验证f(X0)是否等于0，等于0则X0就是方程的解；计算端点f(a)、f(b)的值：若f(a) · f(X0) < 0，则令b = X0；若f(b) · f(X0) < 0，则令a = X0。
步骤四：判断是否达到精度e；即若|a − b| < e = 0.001，则得到零点近似值a (或b)，否则重复循环前述步骤计算中点值进行迭代。
用以上方法寻找每个方程的根，得到坐标，并计算面积，得到结果。
5.1.2.3 模型求解的结果
表 1 问题 1 的计算结果
| 测线距中心点处的距离/m | 海水深度/m | 覆盖宽度/m | 与前一条测线的重叠率/% |
|---|---|---|---|
| -800 | 90.96 | 315.79 | - |
| -600 | 85.73 | 287.59 | 33.64 |
| -400 | 80.43 | 279.43 | 29.63 |
| -200 | 75.18 | 261.24 | 25.03 |
| 0 | 70.00 | 242.96 | 19.75 |
| 200 | 64.79 | 224.86 | 13.71 |
| 400 | 59.56 | 206.70 | 6.84 |
| 600 | 54.25 | 188.53 | -1.39 |
| 800 | 49.07 | 170.36 | -11.16 |
由上表分析可知，随着测线距中心点处的距离逐渐增大，海水深度在逐渐降低，覆盖宽度也在逐渐减小，重叠率出现负数情况，即发生漏测现象。距测线中心 -600 m处为最大重叠率位置，重叠率为33.64%，此时覆盖范围为287.59m，距测线中心600 m和800 m处存在漏测现象，重叠率分别为-1.39%和-11.16%，此时的覆盖宽度分别为188.53 m和170.36 m。
5.1.3 问题结论
在测线距中心点处的距离为-800m时，对应的海水深度为90.96m，多波束测深条带的覆盖宽度 W = 315.79 m W=315.79\,\text{m} W=315.79m；在测线距中心点处的距离为-600m时，对应的海水深度为85.73m，多波束测深条带的覆盖宽度 W = 287.59 m W=287.59\,\text{m} W=287.59m，与前一条测线的重叠率为33.64%；在测线距中心点处的距离为-400m时，对应的海水深度为80.43m，多波束测深条带的覆盖宽度 W = 279.43 m W=279.43\,\text{m} W=279.43m，与前一条测线的重叠率为29.63%；在测线距中心点处的距离为-200m时，对应的海水深度为75.18m，多波束测深条带的覆盖宽度 W = 261.24 m W=261.24\,\text{m} W=261.24m，与前一条测线的重叠率为25.03%；在测线距中心点处的距离为0m时，对应的海水深度为70m，多波束测深条带的覆盖宽度 W = 242.96 m W=242.96\,\text{m} W=242.96m，与前一条测线的重叠率为13.75%；在测线距中心点处的距离为200m时，对应的海水深度!-- 本文件由 PDF 自动转换为 Markdown，可在 CSDN / Typora / Obsidian 等 Markdown 编辑器中继续编辑。 -->
多波束测线问题
摘要
论文PDF版下载
本文针对多波束测线问题，建立了二维几何模型，解决了特定位置下的指标计算问
题；建立了三维解析几何模型，完成了特定位置下覆盖宽度的求解；建立了单目标优化
模型，给出了最优的测线规划方案；建立了多目标优化模型，给出了最优的多波束测线
方案。
针对问题一，建立了二维几何模型，解决了特定位置下的指标计算问题。首先，以
海域中心为原点，以航向垂直的方向为 x 轴，铅直方向为 y 轴建立直角坐标系；其次，
求出斜坡直线方程和给定位置处的波束边界方程，进而利用二分法求解波束与斜坡的交
点坐标，再由交点坐标实现波束覆盖宽度及相邻条带重叠率的计算。计算结果表明：距
测线中心
−
600
-600
−600m 处为最大重叠率位置，重叠率为
33.64
%
33.64\%
33.64%，此时覆盖宽度为
287.59
287.59
287.59m，距
测线中心
600
600
600m 和
800
800
800m 处存在漏测现象，重叠率分别为
−
1.39
%
-1.39\%
−1.39% 和
−
11.17
%
-11.17\%
−11.17%，此时的
覆盖宽度为 188.53m 和 170.36m。重叠率对斜坡角度、换能开角和海水深度的灵敏度分
析表明：重叠率对海水深度最为灵敏，灵敏度为 3.722。最后，利用黄金分割算法对上
述结果进行了重新求解，计算结果误差均小于 0.001，表明计算结果的可信度。
针对问题二，建立了三维解析几何模型，完成了特定位置处覆盖宽度的求解。首先，
以海域中心为原点，以海底坡面法向量在海平面投影向量为 x 轴，海平面内与 x 轴垂
直向右方向为 y 轴，铅直方向为 z 轴建立直角坐标系；其次，求出海底斜坡平面方程和
给定位置波束边界的直线方程，进而利用牛顿迭代法求解波束直线与斜坡平面的交点坐
标，再由交点坐标实现覆盖宽度的计算。计算结果表明：测线方向的夹角为 270° 时，距
离海域中心 2.1 海里处覆盖宽度达到最大值 770.67m；同一点处测线方向夹角为 90° 时
覆盖宽度最小达到 63.35m。最后，分析了覆盖宽度对测线间距、测线方向角和海水深度
的灵敏度，结果表明：覆盖宽度对海水深度最为灵敏，灵敏度的值为 1.0085。
针对问题三，建立了单目标优化模型，给出了最优的测线规划方案。首先，在问题
二直角坐标系的基础上，给出平行测线的直线方程，并求出测线累计距离，以该距离最
小为目标；其次，求出各测线与边界的交点，并求得这些交点处的重叠率，以重叠率满
足 10% ∼ 20% 为约束条件，建立单目标优化模型，利用差分进化算法进行求解。计算结
果表明，最优的测线为非等间距平行线，航向为南北方向，测线总长度为 144.5km，整
个区域的平均重叠率为 15%。最后，利用追赶法对模型进行重新求解，对比差分进化算
法求得的测线总长度结果，长度误差小于 1.2%，验证了计算结果的可信度。
针对问题四，建立了多目标优化模型，给出了最优的多波束测线方案。首先，利用
二元三次幂函数对附件中的海深数据进行拟合，得到拟合的曲面方程；其次，利用航向
角求出累计测线距离，并以该距离最小为目标 1。由航线求出航线在拟合曲面上的投影
方程，并利用波束边界直线方程确定覆盖宽度范围，进而实现重叠率大小的计算，以重
叠率大于 20% 的比例最小为目标 2，建立多目标优化模型。在量纲归一化后，利用加权
法将多目标转为单目标优化模型，利用混合模拟退火禁忌搜索算法进行求解。计算结果
表明：测线总长度为 331.47km，漏测海区占总海域面积的 4.48%，重叠率超过 20% 的总
长度为 7.18km，占比约为 2.16%，主要位于海洋区域的西南角。最后，采用灰狼优化算
法对模型的求解结果进行检验，结果表明：测线总长度误差小于 5km，漏测面积占比和
重叠率总长度的误差率小于 3% 。
关键词：几何模型；黄金分割法；优化模型；差分进化算法；混合模拟退火禁忌搜索法
一、问题重述
一、问题重述
1.1 问题背景
由于陆地资源日益减少，人口规模迅速增长，经济发展的紧迫需求，人们迫切需要
开采、开发和利用海洋资源。海洋已成为 21 世纪的现实和战略关注点，并承载着更多的
重要意义。而从事海岸工程建设，地质调查与开发等海洋活动，均需要海洋测绘勘探工
作。海底地形测量作为这项工作的基石，需要综合运用多学科的理论、方法和技术。为
了更加深入地了解海洋环境并精确掌握海底地形地貌，通常会采用多波束测深技术。由
于其全覆盖、全水深、高精度和高分辨率等优势，多波束测深系统已成为当前主流的测
深设备，在水下地形测量、海底地貌勘察和水下目标检测等领域得到广泛应用。
相较于单波束测深，多波束测深系统在垂直范围沿航迹上能够快速而准确地测量水
下目标和地形起伏的变化。它将传统的点、线探测扩展至面的层次，实现了数字化记录
和立体化自动测图。通过克服单波束测深的采样间距过大导致海底信息反映粗糙以及发
射波束角度过大导致微小地形变化的深度误差等缺陷。因此，可以更准确、更广泛地获
取航道的相关数据，从而全面提升航道测绘技术的精确度。
1.2 问题重述
多波束测深系统中，多波束测深条带的覆盖宽度 W 会随着换能器的开角 θ 和水深
D 的变化而发生变化。本文的数学模型主要解决以下问题：
- 在海底不平的情况下，当测量船测线方向垂直的平面和海底坡面的夹角为 α 时，建立多波束测深的覆盖宽度及相邻条带之间重叠率的数学模型。当开角 θ=120°，坡度 α = 5 ∘ \alpha=5^{\circ} α=5∘，海水深度为 70 米时，根据建立的数学模型求解相关参数，并将结果填入表 1。
- 在海底不平的情况下，考虑一个矩形待测海域，当测量船测线方向与海底坡面的法向在水平面上投影的夹角为 β 时，建立多波束测深覆盖宽度的数学模型。当换能器的开角为 θ=120°，坡度为 α=1.5°，海域中心点处的海水深度为 120 m 时，根据覆盖宽度模型计算不同测量方向夹角的覆盖宽度，并将结果填入表 2。
- 在海底不平的情况下，在一个南北长 2 海里、东西宽 4 海里的矩形海域内，海域中心点处的海水深度为 110 m，呈西深东浅的坡度 α=1.5°。多波束换能器的开角为 θ=120°。现需要设计一组测量线，使其长度最短且可以完全覆盖整个海域，同时满足相邻条带之间重叠率在 10% 20% 之间。
- 在海底不平的情况下，根据已知的海水深度数据，要求设计出满足如下要求的测线：沿测线扫描形成的条带尽可能地覆盖整个待测海域；相邻条带之间的重叠率尽量控制在 20% 以下；测线的总长度尽可能短。以该设计测线，要求计算出测线的总长度；漏测海区占总待测海域面积的百分比；在重叠区域中，重叠率超过 20% 部分的总长度。
二、问题分析
问题整体思路较为清晰，层层递进。问题一要求测量间距为 200 米是关于多波束测深覆盖宽度和相邻条带覆盖率参数的计算问题。问题二建立在问题一相关参数求解的基础上，增加一个矩形待测海域，测量间距由 200 米变为 0.3 海里，是关于在不同测线方向夹角下，覆盖宽度参数的计算问题。问题三是建立在问题二矩形待测海域的基础上，
设置矩形海域宽度，是关于如何优化待测海域测线的问题。问题四是建立在问题三优化模型的基础上，是关于如何优化设计出满足要求的测线和相关指标的计算问题。
2.1 问题一的分析
问题一要求在海水起伏变化大，每隔 200m 测量，测量船测线方向垂直的平面和海
底坡面的夹角为
α
\alpha
α 的情况下，建立多波束测深覆盖宽度及相邻条带之间重叠率的数学模
型，以求解特定位置下的指标。从题目内容可知，要求计算出测量船波束间的覆盖宽度
W
W
W，即计算测量船波束边界点与斜坡交点之间的距离，需要得知波束与斜坡的交点坐标，
再根据交点坐标实现重叠率的计算。因此，考虑使用二维几何模型 [1]
来解决问题。假设
海面近似为平面，以海域中心为原点，以航向垂直的方向为
x
x
x 轴，铅直方向为
y
y
y 轴建立
直角坐标系，并列出波束边界点与测量点的坐标公式；其次，根据上述三点坐标，列出
斜坡直线方程和给定位置处的波束边界方程，构建多束波覆盖宽度的数学模型。求解方
程的一般计算方法有高斯消元法 [2]
，二分法 [3]
，迭代法 [4]
等算法，本题宜于选用二分
法求解。计算出波束与斜坡的交点坐标，以此实现测线重叠率的计算。最终求解出在不
同测线距中心点处距离下的海水深度，覆盖宽度，重叠率。最后，分析了重叠率对斜坡
角度、换能器开角和海水深度的灵敏度。利用黄金分割法 [5]
对上述结果进行了重新求
解，计算结果误差均较小，表明计算结果的可信度较高。
2.2 问题二的分析
问题二在问题一的基础上，测量间距变为
0.3
0.3
0.3 海里，要求考虑一个矩形待测海域，
即测线方向与海底坡面的法向在水平面上投影夹角为
β
\beta
β 的平面的情况下，建立多波束测深覆盖宽度的数学模型以求解特定位置下覆盖的宽度。从题目内容可知，要求覆盖的
宽度，求出波束直线与斜坡平面的交点坐标即可。因此，考虑使用三维几何模型来解决
问题。首先，以海域中心为原点，以平行于海底交线方向为
x
x
x 轴，垂直交线方向为
y
y
y 轴，
铅直方向为
z
z
z 轴建立直角坐标系，并列出波束边界点与测量点的坐标公式；其次，求出
海底斜坡平面方程和给定位置波束边界的直线方程。采用牛顿迭代法[6]求解波束直线
与斜坡平面的交点坐标，再由交点坐标实现覆盖宽度的计算，最终求出不同测线方向夹
角下，不同测量间距的覆盖宽度。最后，分析了覆盖宽度对测线间距、测线方向角和海
水深度的灵敏度 [7]
，验证了模型的合理性。
2.3 问题三的分析
问题三在问题二矩形待测海域的基础上，设置了具体的长宽矩形海域，要求以使测线累计距离最小为目标，以特定的海域中心处的海水深度，坡度，多波束换能器的开角，相邻条带重叠率满足
10
%
10\%
10% 到
20
%
20\%
20% 为约束条件，建立优化模型。从题目内容可
知，该题属于测线长度的优化问题，考虑使用单目标优化 [7]
模型来解决问题。首先，在
问题二直角坐标系的基础上，给出平行测线的直线方程，并求出测线累计距离，以该距离最小为目标；其次，求出各测线与边界的交点，并求得这些交点处的重叠率，以重叠率满足
10
%
∼
20
%
10\% \sim 20\%
10%∼20% 为约束条件，建立单目标优化模型。解决此类模型的一般计算方法
有梯度下降法，粒子群优化 [9]
，差分进化算法 [10]
等方法，本题采用差分进化算法进行
求解，计算得出最优的测线总长度，此时最大最小的重叠率分别对应的海水深度。最后，
利用追赶法 [11]
对模型的求解结果进行检验，计算结果与差分进化算法的求解结果误差
较小，验证了计算结果的可信度性。
2.4 问题四的分析
问题四在问题三矩形待测海域的基础上，给定另一不同矩形海域的海水深度数据，
要求对该数据进行拟合处理，以此结果设计出距离最短，重叠率大于 20% 的比例最小
的测线，给出最优的多波束测线方案。从题目内容可知，该题属于多波束测线方案的优化问题，考虑使用多目标优化模型进行求解。首先，利用二元三次幂函数对附件中的海水深度数据进行拟合，得到拟合的曲面方程；其次，利用航向角求出累计测线距离，并以该距离最小为目标1。由航线求出航线在拟合曲面上的投影方程，并利用波束边界直线方程确定覆盖宽度范围，进而计算重叠率，以重叠率大于20%的比例最小为目标2，建立多目标优化模型。在等权条件下，将多目标转为单目标优化模型。采用混合模拟退火禁忌搜索算法 [12] 对优化模型进行求解，得出最优的测线总长度，漏测海区占总海域的面积比等参数，最后，采用灰狼优化 [13] 算法对模型结果进行重新检验，误差较小，验证了计算结果的可靠性。
三、基本假设
- 假设海平面近似平面；
- 假设不考虑测量船的姿态；
- 假设海面平静；
- 假设海水折射率均匀分布。
四、符号说明
α
\alpha
α：坡度
θ
\theta
θ：换能器的开角
P
i
P_i
Pi：测线与海底坡面的左交点
Q
i
Q_i
Qi：测线与海底坡面的右交点
M
0
M_0
M0：沿海底坡面方向的直线
Mi: 沿测线方向的第一条直线
Mi+1: 沿测线方向的第二条直线
D: 海水深度
a: 测线距离中心点处的距离
五、模型建立与求解
5.1 问题一的模型建立与求解
本题为解决特定位置下的指标计算问题，建立二维几何模型。以海域中心为原点，航向方向为 x x x 轴，铅直方向为 y y y 轴，建立平面直角坐标系。其次，根据坐标系示意图，求出斜坡直线方程和特定位置处波束边界方程，利用二分法求解出测量船的波束边界和斜坡的交点坐标，以此计算出重叠率和条带的覆盖宽度。
5.1.1 模型建立
5.1.1.1 二维几何模型的确立
本题旨在建立二维几何模型，求解特定位置下的指标。因此，以海域中心为原点，
以航向垂直的方向为 x 轴，铅直方向为 y 轴，建立如下图的直角坐标系：
图 1 测量船波束与斜坡示意图
由上图可知，在 xoy 平面上，
测量船位于海平面上，
设测量船
i
i
i 所在的坐标为
(
x
i
,
0
)
(x_i, 0)
(xi,0)。
坡地倾角为
α
\alpha
α，沿海底坡面方向的直线
M
0
M_0
M0 的斜率为
tan
α
\tan\alpha
tanα，故沿海底坡面方向的直线
M
0
M_0
M0 的表达式为:
M
0
:
y
=
tan
α
⋅
x
−
D
M_0:\ y = \tan\alpha \cdot x - D
M0: y=tanα⋅x−D
其中
D
=
70
m
D=70\text{m}
D=70m 为海水深度。利用 i 点的横坐标代入该方程即可得到海水深度。
由于平面
x
O
y
xOy
xOy 垂直于测线方向，且垂线平分多波束换能器开角
θ
\theta
θ，故沿测线方向的第一条直线
M
i
M_i
Mi 与平面的夹角为
π
−
θ
\pi - \theta
π−θ，该直线的斜率为
tan
(
π
−
θ
)
\tan(\pi - \theta)
tan(π−θ)，故沿测线方向的第一条直线
M
i
M_i
Mi 的表达式为：
M
i
:
y
=
tan
(
π
−
θ
)
(
x
−
x
i
)
M_i:\ y = \tan(\pi - \theta)(x - x_i)
Mi: y=tan(π−θ)(x−xi)
同理, 沿测线方向的第二条直线
M
i
+
1
M_{i+1}
Mi+1 与平面的夹角为
π
−
θ
+
θ
=
π
+
θ
\pi - \theta + \theta = \pi + \theta
π−θ+θ=π+θ, 斜率为
tan
(
π
+
θ
)
\tan(\pi + \theta)
tan(π+θ)，故沿测线方向的第二条直线
M
i
+
1
M_{i+1}
Mi+1 的表达式为:
M
i
+
1
:
y
=
tan
(
π
+
θ
)
(
x
−
x
i
)
M_{i+1} : y = \tan(\pi + \theta)(x - x_i)
Mi+1:y=tan(π+θ)(x−xi)
由沿海底坡面方向的直线方程 M0，与沿测线方向的第一条直线方程 Mi 和沿测线方向的第二条直线方程 Mi+1 联立，可分别求出两条测线与海底坡面的交点 Pi (xpi, ypi)，Qi (xqi, yqi) 的坐标，即：
Pi (xpi, ypi) ：
y = tan α · x − D
y = tan
(
π
−
θ
)
(x − xi)
Qi (xqi, yqi) ：
y = tan α · x − D
y = tan
(
π
+
θ
)
(x − xi)
利用交点 Pi, Qi 的坐标，由欧式距离计算出 PiQi 的长度，即为为条带的覆盖宽度
W，即
W = |PiQi| =
√
(xpi − xqi)2 + (ypi − yqi)2
对于 i + 1 节点，可利用上述分析得到 Pi+1, Qi+1 节点的坐标，利用重叠关系可得重叠率表达式如下。根据题意可知，相邻条带之间的重叠率定义为
η
=
1
−
d
W
=
∣
P
i
Q
i
∣
−
∣
P
i
+
1
Q
i
∣
∣
P
i
Q
i
∣
\eta = 1 - \frac{d}{W} = \frac{|P_iQ_i| - |P_{i+1}Q_i|}{|P_iQ_i|}
η=1−Wd=∣PiQi∣∣PiQi∣−∣Pi+1Qi∣
5.1.1.2 二维几何模型表达式
综上所述，最终的二维模型表达式如下：
{ W = ∣ P i Q i ∣ = ( x p i − x q i ) 2 + ( y p i − y q i ) 2 η = 1 − d W = ∣ P i Q i ∣ − ∣ P i + 1 Q i ∣ ∣ P i Q i ∣ y = tan ( π − θ ) ( x − x i ) y = tan ( π + θ ) ( x − x i ) x i ∈ { − 800 , − 600 , ⋯ , 600 , 800 } \begin{cases} W = |P_iQ_i| = \sqrt{(x_{p_i} - x_{q_i})^2 + (y_{p_i} - y_{q_i})^2} \\ \eta = 1 - \frac{d}{W} = \frac{|P_iQ_i| - |P_{i+1}Q_i|}{|P_iQ_i|} \\ y = \tan\left(\pi - \theta\right) (x - x_i) \\ y = \tan\left(\pi + \theta\right) (x - x_i) \\ x_i \in \{-800, -600, \cdots, 600, 800\} \end{cases} ⎩ ⎨ ⎧W=∣PiQi∣=(xpi−xqi)2+(ypi−yqi)2η=1−Wd=∣PiQi∣∣PiQi∣−∣Pi+1Qi∣y=tan(π−θ)(x−xi)y=tan(π+θ)(x−xi)xi∈{−800,−600,⋯,600,800}
其中， D = 70 m D = 70\ \text{m} D=70 m， α = 1.5 ∘ \alpha = 1.5^{\circ} α=1.5∘， θ = 120 ∘ \theta = 120^{\circ} θ=120∘
5.1.2 模型求解
5.1.2.1 模型求解的原理
根据斜坡方程和测量船波束边界方程，求出两个的交点坐标，根据交点坐标计算出
覆盖宽度和重叠率。本题采用二分法来求解交点坐标。
二分法又称二分搜索法、折半搜索法或对数搜索法，是在一组顺序排列的元素中查找某一特定元素的搜索方法。它通过不断对分数据区间，将不可能的解集排除，得到满足误差条件的近似解。在每次二分迭代时都会将范围缩小一半，迭代次数少，查找速度快，性能较高。
5.1.2.2 模型求解的过程
步骤一：先确认定义域的区间
[
a
,
b
]
[a, b]
[a,b]，输入端点及精度值确定区间
[
−
800
,
800
]
[-800, 800]
[−800,800]，验证
f
(
−
800
)
⋅
f
(
800
)
<
0
f(-800) \cdot f(800) < 0
f(−800)⋅f(800)<0，如果不成立则一般不存在根或是区间太过宽泛，再给定精度
e
=
0.001
e=0.001
e=0.001；
步骤二：求区间
(
−
800
,
800
)
(-800, 800)
(−800,800) 的中点 X0 = (a + b)/2，即 X0 = 0；
步骤三：计算 f(X0)，验证 f(X0) 是否等于 0，等于 0 则 X0 就是方程的解；计算端点
f(a)、f(b) 的值：若 f(a) · f(X0) < 0，则令 b = X0；若 f(b) · f(X0) < 0，则令 a = X0。
步骤四：判断是否达到精度 e；即若 |a − b| < e = 0.001，则得到零点近似值 a (或 b)，否则
重复循环前述步骤计算中点值进行迭代。
用以上方法寻找每个方程的根，得到坐标，并计算面积，得到结果。
5.1.2.3 模型求解的结果
表 1 问题 1 的计算结果
| 测线距中心点处的距离/m | 海水深度/m | 覆盖宽度/m | 与前一条测线的重叠率/% |
|---|---|---|---|
| -800 | 90.96 | 315.79 | - |
| -600 | 85.73 | 287.59 | 33.64 |
| -400 | 80.43 | 279.43 | 29.63 |
| -200 | 75.18 | 261.24 | 25.03 |
| 0 | 70.00 | 242.96 | 19.75 |
| 200 | 64.79 | 224.86 | 13.71 |
| 400 | 59.56 | 206.70 | 6.84 |
| 600 | 54.25 | 188.53 | -1.39 |
| 800 | 49.07 | 170.36 | -11.16 |
由上表分析可知，随着测线距中心点处的距离逐渐增大，海水深度在逐渐降低，覆盖宽度也在逐渐减小，重叠率出现负数情况，即发生漏测现象。距测线中心 -600 m 处
为最大重叠率位置，重叠率为 33.64%，此时覆盖范围为 287.59m，距测线中心 600 m 和 800 m 处存在漏测现象，重叠率分别为
−
1.39
%
-1.39\%
−1.39% 和
−
11.16
%
-11.16\%
−11.16%，此时的覆盖宽度分别为 188.53 m 和 170.36 m。
5.1.3 问题结论
在测线距中心点处的距离为-800m 时，对应的海水深度为 90.96m，多波束测深条
带的覆盖宽度
W
=
315.79
m
W=315.79\,\text{m}
W=315.79m；在测线距中心点处的距离为-600m 时，对应的海水深度为
85.73m，多波束测深条带的覆盖宽度
W
=
287.59
m
W=287.59\,\text{m}
W=287.59m，与前一条测线的重叠率为 33.64%；在测
线距中心点处的距离为-400m 时，对应的海水深度为 80.43m，多波束测深条带的覆盖宽
度
W
=
279.43
m
W=279.43\,\text{m}
W=279.43m，与前一条测线的重叠率为 29.63%；在测线距中心点处的距离为-200m 时，
对应的海水深度为 75.18m，多波束测深条带的覆盖宽度
W
=
261.24
m
W=261.24\,\text{m}
W=261.24m，与前一条测线的重
叠率为 25.03%; 在测线距中心点处的距离为 0m 时，对应的海水深度为 70m, 多波束测深
条带的覆盖宽度
W
=
242.96
m
W=242.96\,\text{m}
W=242.96m, 与前一条测线的重叠率为 13.75%; 在测线距中心点处的距
离为 200m 时，对应的海水深度为 64.79m, 多波束测深条带的覆盖宽度 W=224.86m, 与
前一条测线的重叠率为 13.71%; 在测线距中心点处的距离为 400m 时，对应的海水深度
为 59.56m, 多波束测深条带的覆盖宽度 W=206.70m, 与前一条测线的重叠率为 6.84%; 在
测线距中心点处的距离为 600m 时，对应的海水深度为 54.25m, 多波束测深条带的覆盖
宽度 W=188.53m, 与前一条测线的重叠率为-1.3910%; 在测线距中心点处的距离为 800m
时，对应的海水深度为 49.07m, 多波束测深条带的覆盖宽度 W=170.36m, 与前一条测线
的重叠率为-11.16%。
5.1.4 检验分析
5.1.4.1 灵敏度分析
考虑在实际生活中，由于海底凹凸不平，海面起伏不定，导致换能器的开角与海水深度存在偏差，因此为验证斜坡角度
α
=
1.5
∘
\alpha = 1.5^{\circ}
α=1.5∘、换能器开角
θ
=
120
∘
\theta = 120^{\circ}
θ=120∘ 和海水深度 D 对重叠率
η
\eta
η 的影响，将重叠率的变化幅度作为指标，取测线距中心点处距离为 0，将斜坡角度、换能器的开角、海水深度依次增加 0.2°、0.2°、
0.2
m
0.2\ \text{m}
0.2 m，观察斜坡角度、换能器的开角、海水深度变化后重叠率的变化幅度。
在坡度为 1.5° 时，每次以 0.2° 为单位依次增加，经计算，坡度每次变换后，所得重叠率为 19.78、20.06、20.33、20.61、20.89、21.14。
在海水深度为
70
m
70\ \text{m}
70 m 时，每次以 0.2m 为单位依次增长，经计算得，海水深度每次变
换后，所得重叠率为 19.78、19.96、20.43、20.63、20.85。
在换能器开角为 120° 时，每次以 0.2° 为单位依次增加，计算结果如下表所示：
表 2 换能器的开角对重叠率的影响结果表
| 换能器开角/° | 120 | 120.2 | 120.4 | 120.6 | 120.8 | 121.0 |
|---|---|---|---|---|---|---|
| 重叠率/% | 19.78 | 20.09 | 20.40 | 20.71 | 21.01 | 21.32 |
根据上述换能器开角变化后的重叠率，将所得结果用如下的折线图表示：
图 2 换能器开角的变化对重叠率的示意图
由上图可知，以换能器开角 120° 为起始点，依次增加 0.2°，观察到相邻条带间的重
叠率在依次增长且幅度较大。灵敏性数据转换成百分比公式如下：
△y
y
△x
x
→
dy
dx
·
x
y
根据上述公式计算可知，在换能器、坡度和海底深度中，重叠率对换能器开角最为
灵敏，灵敏度的值为 3.772，即改变转换器开角时，重叠率变化的幅度最大。
5.1.4.2 方法比较
采用黄金分割法对该模型的结果进行重新求解，对比二分法的求解结果。其理论基础如下：
步骤一：选取初始数值，确定初始区间 [a0, b0]，即 [−400, 400]，给出精度
δ
>
0
\delta > 0
δ>0，即
δ
=
0.001
\delta=0.001
δ=0.001。
步骤二：计算初始的两个试点（规定在区间 [ak, bk] 上，黄金分割点用 xk+1 表示，黄金分割点的对称点用 x′
k+1 表示）。
计算
{
x
1
=
a
0
+
0.618
(
b
0
−
a
0
)
x
1
′
=
a
0
+
0.382
(
b
0
−
a
0
)
\begin{cases} x_1 = a_0 + 0.618 (b_0 - a_0) \\ x'_1 = a_0 + 0.382 (b_0 - a_0) \end{cases}
{x1=a0+0.618(b0−a0)x1′=a0+0.382(b0−a0)
并计算出相应试点的函数值
φ
(
x
1
)
\varphi(x_1)
φ(x1) 和
φ
(
x
1
′
)
\varphi(x'_1)
φ(x1′)。
步骤三：此时规定 k=0，比较目标函数值。当
φ
(
x
k
+
1
′
)
≤
φ
(
x
k
+
1
)
\varphi(x'_{k+1}) \leq \varphi(x_{k+1})
φ(xk+1′)≤φ(xk+1) 时，按照下列规则缩小搜索范围：
取
{
a
k
+
1
=
a
k
b
k
+
1
=
x
k
+
1
\begin{cases} a_{k+1} = a_k \\ b_{k+1} = x_{k+1} \end{cases}
{ak+1=akbk+1=xk+1
步骤四：算法在 k=k+1 下转至第三步。当
φ
(
x
k
+
1
′
)
>
φ
(
x
k
+
1
)
\varphi(x'_{k+1}) > \varphi(x_{k+1})
φ(xk+1′)>φ(xk+1) 时，按照下列规则缩小搜索范围：取
{
a
k
+
1
=
x
k
+
1
′
b
k
+
1
=
b
k
\begin{cases} a_{k+1} = x'_{k+1} \\ b_{k+1} = b_k \end{cases}
{ak+1=xk+1′bk+1=bk
若 ∆ = bk+1 − ak+1 / b0 − a0 < δ, 则算法终止。若 ∆ = bk+1 − ak+1 / b0 − a0 > δ, 则要计算新的试点, 即
{
x′
k+2 = xk+1
xk+2 = ak+1 + 0.618 (bk+1 − ak+1)
经过以上的对函数的重新计算，得到的结果误差为 0.000162，小于 0.001，证明了计算结果的准确性。
5.1.5 小结
本题建立了二维几何模型，解决了特定位置下的指标计算问题。首先，以海域中心为原点，以航向垂直的方向为 x 轴，铅直方向为 y 轴建立直角坐标系；其次，求出斜坡直线方程和给定位置处的波束边界方程，进而利用二分法求解波束与斜坡的交点坐标，再由交点坐标实现测线重叠率的计算。计算结果表明：距测线中心
−
600
m
-600\,\text{m}
−600m 处为最大重叠
率位置，重叠率为
33.64
%
33.64\%
33.64%，此时覆盖范围为
287.59
m
287.59\,\text{m}
287.59m，距测线中心 600m 和
800
m
800\,\text{m}
800m 处存在
漏测现象，重叠率分别为
−
1.39
%
-1.39\%
−1.39% 和 −11.17%，此时的覆盖范围为 188.53m 和 170.36m。
最后，分析了重叠率对斜坡角度、换能开角和海水深度的灵敏度，结果表明：换能开角
对海水深度最为灵敏，灵敏度达 3.722。利用黄金分割算法对上述结果进行了重新求解，
计算结果误差均小于 0.001，表明计算结果的可信度。
5.2 问题二的模型建立与求解
本题为求解特定位置下的覆盖宽度，建立了三维解析几何模型。首先，以海域中心
为原点，建立空间直角坐标系，其次，根据坐标系相关示意图，求出海底斜坡平面方程
和给定位置的波束边界直线方程，利用牛顿迭代法求解出波束边界直线和斜坡平面的交
点坐标，以此计算出在不同测向方向夹角
θ
\theta
θ下的不同测量间距的覆盖宽度。
5.2.1 模型建立与求解
5.2.1.1 三维解析几何模型的确立
本题旨在建立三维解析几何模型，求解出特定位置下的覆盖宽度。因此，以海域中
心为原点，坡面法线向量在水平上的投影正向量为 x 轴，垂直于海平面方向为 z 轴，垂
直于坡面法向量与投影所构成的平面为 y 轴，建立空间直角坐标系，如下图所示：
图 3 矩形待测海域示意图
由上图可知，测量船的测线方向为
v
⃗
\vec{v}
v，作两条平行于测线方向的直线
w
,
w
+
1
w,w+1
w,w+1。在海
平面上作一直线 w0，使其垂直于直线 w，w+1。z 轴垂直于海平面，w0 在海平面内，故
w0 垂直于 z 轴。
w
0
w_0
w0 垂直于
w
w
w，
w
w
w 平行于测线方向，故
w
0
w_0
w0 垂直于测线方向。根据图中的角度关系可得，直线
w
0
w_0
w0 的方向向量
s
⃗
=
(
cos
(
β
−
π
)
,
sin
(
β
−
π
)
,
0
)
\vec{s} = \left( \cos\left(\beta - \pi\right), \sin\left(\beta - \pi\right), 0 \right)
s=(cos(β−π),sin(β−π),0)，且该直线过坐标原点，故可计算出直线
w
0
w_0
w0 的方程为
x
cos
(
β
−
π
)
=
y
sin
(
β
−
π
)
=
z
\frac{x}{\cos\left(\beta - \pi\right)} = \frac{y}{\sin\left(\beta - \pi\right)} = z
cos(β−π)x=sin(β−π)y=z
为进一步确立测量船的位置，需要确立以海域原点为圆心，以
R
=
i
d
R=id
R=id 为半径，
d
d
d 为测量船距海域中心点处的距离，
i
d
id
id 表示第
i
i
i 条测量间距，
i
=
0
,
…
,
8
i=0,\dots,8
i=0,…,8，得出以下海平面方程：
{
x
2
+
y
2
=
R
2
Z
=
0
\begin{cases} x^2 + y^2 = R^2 \\ Z = 0 \end{cases}
{x2+y2=R2Z=0
将直线
w
0
w_0
w0 方程与海平面方程联立，求出测量船的位置坐标
R
(
x
i
,
y
i
,
z
i
)
R (x_i, y_i, z_i)
R(xi,yi,zi)，即
R
(
x
i
,
y
i
,
z
i
)
:
{
x
cos
(
β
−
π
)
=
y
sin
(
β
−
π
)
=
z
x
2
+
y
2
=
R
2
Z
=
0
R (x_i, y_i, z_i) : \begin{cases} \dfrac{x}{\cos\left(\beta - \pi\right)} = \dfrac{y}{\sin\left(\beta - \pi\right)} = z \\ x^2 + y^2 = R^2 \\ Z = 0 \end{cases}
R(xi,yi,zi):⎩
⎨
⎧cos(β−π)x=sin(β−π)y=zx2+y2=R2Z=0
由图可知，坡面法线向量与竖直方向的夹角为
α
\alpha
α，故坡面的法线向量为
n
⃗
=
(
sin
α
,
0
,
cos
α
)
\vec{n} = (\sin\alpha, 0, \cos\alpha)
n=(sinα,0,cosα)。
假设矩形待测海域上一点为 Y (0, 0, D),D 为海水深度，根据点法式方程，求得待测海域平面 E 的表达式为
E : sin α · x + cos α (z − D) = 0
根据图中的角度关系，求得沿测线方向的第一条直线 Mi 的方向向量
s
⃗
=
(
sin
θ
cos
(
3
π
−
β
)
,
sin
θ
sin
(
3
π
−
β
)
,
cos
θ
)
\vec{s} = (\sin\theta\cos(3\pi - \beta), \sin\theta\sin(3\pi - \beta), \cos\theta)
s=(sinθcos(3π−β),sinθsin(3π−β),cosθ)。根据测量船的位置坐标 R = (xi, yi, zi), 计算出沿测线方向的第一条直线 Mi 的表达式为
Mi :
x
−
x
i
sin
θ
cos
(
3
π
−
β
)
=
y
−
y
i
sin
θ
sin
(
3
π
−
β
)
=
z
−
z
i
cos
θ
\frac{x - xi}{\sin\theta\cos(3\pi - \beta)} = \frac{y - yi}{\sin\theta\sin(3\pi - \beta)} = \frac{z - zi}{\cos\theta}
sinθcos(3π−β)x−xi=sinθsin(3π−β)y−yi=cosθz−zi
同理可知，沿测线方向的第二条直线 Mi+1 的方向向量为
s
⃗
:
(
−
sin
θ
cos
(
3
π
−
β
)
,
cos
θ
)
\vec{s} : (-\sin\theta\cos(3\pi - \beta), \cos\theta)
s:(−sinθcos(3π−β),cosθ)。根据测量船的位置坐标 R = (xi, yi, zi)，可计算出沿测线方向的第二条直线 Mi+1 的直线方程为
Mi+1 :
x
−
x
i
−
sin
θ
cos
(
3
π
−
β
)
=
y
−
y
i
−
sin
θ
sin
(
3
π
−
β
)
=
z
−
z
i
cos
θ
\frac{x - xi}{-\sin\theta\cos(3\pi - \beta)} = \frac{y - yi}{-\sin\theta\sin(3\pi - \beta)} = \frac{z - zi}{\cos\theta}
−sinθcos(3π−β)x−xi=−sinθsin(3π−β)y−yi=cosθz−zi
将测量船的位置坐标 R，矩形待测海域平面 E，沿测线方向的第一条直线 Mi 联立方程组，求出测量船第一条波束直线 Mi 与斜坡平面 E 的交点坐标 Pi (xpi, ypi) ，联立方程如下所示：
Pi (xpi, ypi) ：
{
R
=
(
x
i
,
y
i
,
z
i
)
sin
α
⋅
x
+
cos
α
(
z
−
D
)
=
0
x
−
x
i
sin
θ
cos
(
3
π
−
β
)
=
y
−
y
i
sin
θ
sin
(
3
π
−
β
)
=
z
−
z
i
cos
θ
\begin{cases} R = (xi, yi, zi) \\ \sin α · x + \cos α (z − D) = 0 \\ \dfrac{x - xi}{\sin\theta\cos(3\pi - \beta)} = \dfrac{y - yi}{\sin\theta\sin(3\pi - \beta)} = \dfrac{z - zi}{\cos\theta} \end{cases}
⎩
⎨
⎧R=(xi,yi,zi)sinα⋅x+cosα(z−D)=0sinθcos(3π−β)x−xi=sinθsin(3π−β)y−yi=cosθz−zi
将测量船的位置坐标 R，矩形待测海域平面 E，沿测线方向的第二条直线 Mi+1 联
立方程组，求出测量船第二条波束直线 Mi+1 与斜坡平面 E 的交点坐标 Qi (xQi, yQi) ，联
立方程如下所示：
Qi (xQi, yQi) ：
R = (xi, yi, zi)
sin α · x + cos α (z − D) = 0
x − xi
− sin θ
cos
(3π
− β
) =
y − yi
− sin θ
sin
(3π
− β
) =
z − zi
cos θ
根据欧式距离，计算出 Pi, Qi 的距离，即条带的覆盖宽度。
W = |PiQi|
根据题目中相邻条带之间重叠率的定义，进一步推导出重叠率，即：
η =
|PiQi| − |PiQi+1|
|PiQi|
5.2.1.2 解析几何模型表达式
综上所述，最终的三维解析几何模型表达式如下：
{
x
i
cos
(
β
−
π
)
=
y
i
sin
(
β
−
π
)
=
z
i
x
i
2
+
y
i
2
=
d
2
,
其中
d
∈
{
0
,
0.3
,
⋯
,
2.1
}
\begin{cases} \frac{x_i}{\cos\left(\beta - \pi\right)} = \frac{y_i}{\sin\left(\beta - \pi\right)} = z_i \\ x_i^2 + y_i^2 = d^2,\ \text{其中}\ d \in \{0, 0.3, \cdots, 2.1\} \end{cases}
{cos(β−π)xi=sin(β−π)yi=zixi2+yi2=d2, 其中 d∈{0,0.3,⋯,2.1}
{ sin α ⋅ x + cos α ( z − D ) = 0 x − x i sin θ cos ( 3 π 2 − β ) = y − y i sin θ sin ( 3 π 2 − β ) = z − z i cos θ \begin{cases} \sin\alpha \cdot x + \cos\alpha \left(z - D\right) = 0 \\ \frac{x - x_i}{\sin\theta \cos\left(\frac{3\pi}{2} - \beta\right)} = \frac{y - y_i}{\sin\theta \sin\left(\frac{3\pi}{2} - \beta\right)} = \frac{z - z_i}{\cos\theta} \end{cases} {sinα⋅x+cosα(z−D)=0sinθcos(23π−β)x−xi=sinθsin(23π−β)y−yi=cosθz−zi
{ sin α ⋅ x + cos α ( z − D ) = 0 x − x i − sin θ cos ( 3 π 2 − β ) = y − y i − sin θ sin ( 3 π 2 − β ) = z − z i cos θ \begin{cases} \sin\alpha \cdot x + \cos\alpha \left(z - D\right) = 0 \\ \frac{x - x_i}{-\sin\theta \cos\left(\frac{3\pi}{2} - \beta\right)} = \frac{y - y_i}{-\sin\theta \sin\left(\frac{3\pi}{2} - \beta\right)} = \frac{z - z_i}{\cos\theta} \end{cases} {sinα⋅x+cosα(z−D)=0−sinθcos(23π−β)x−xi=−sinθsin(23π−β)y−yi=cosθz−zi
W = ∣ P i Q i ∣ W = |P_iQ_i| W=∣PiQi∣
η = ∣ P i Q i ∣ − ∣ P i Q i + 1 ∣ ∣ P i Q i ∣ \eta = \frac{|P_iQ_i| - |P_iQ_{i+1}|}{|P_iQ_i|} η=∣PiQi∣∣PiQi∣−∣PiQi+1∣
其中， D = 120 m D = 120\ \text{m} D=120 m， α = 1.5 ∘ \alpha = 1.5^{\circ} α=1.5∘， θ = 120 ∘ \theta = 120^{\circ} θ=120∘
5.2.2 模型求解
5.2.2.1 模型求解原理
牛顿迭代法，又被称为牛顿-拉弗森（拉夫逊）方法，也可以称之为牛顿切线法，应
用于数据分析和求解不同类型的方程，牛顿迭代法在计算机等领域具有非常重要的地
位。原始的牛顿法采用迭代的方法来求函数方程的根，从几何意义上简单地来说就是一
个不断求取切线的过程。传统的牛顿法需要计算并使用目标函数的一阶导数和二阶导
数，而简化牛顿法则只使用一阶导数，计算复杂度更低，避免了计算和存储二阶矩阵。
因此，本文将采用简化的牛顿法进行求解三维解析几何模型。
5.2.2.2 模型求解的过程
设
f
(
m
)
f(m)
f(m) 可微并连续的函数, 将 f(m) 在 m0 处 Taylor 展开, 即
f
(
m
)
=
f
(
m
0
)
+
f
′
(
m
0
)
(
m
−
m
0
)
+
f
′
′
(
m
0
)
2
!
(
m
−
m
0
)
2
+
⋯
+
f
(
k
)
(
m
0
)
k
!
(
m
−
m
0
)
k
+
⋯
f(m) = f(m_0) + f'(m_0)(m - m_0) + \frac{f''(m_0)}{2!}(m - m_0)^2 + \dots + \frac{f^{(k)}(m_0)}{k!}(m - m_0)^k + \cdots
f(m)=f(m0)+f′(m0)(m−m0)+2!f′′(m0)(m−m0)2+⋯+k!f(k)(m0)(m−m0)k+⋯
若
f
′
(
m
0
)
≠
0
f'(m_0) \neq 0
f′(m0)=0 ，取其线性部分近似替代 f (m) ，使得
f
(
m
)
=
0
f(m) = 0
f(m)=0 的近似方程
即
m
=
m
0
−
f
(
m
0
)
f
′
(
m
0
)
m
m = m_0 - \frac{f(m_0)}{f'(m_0)}m
m=m0−f′(m0)f(m0)m ，然后使用迭代
m1 = m0 −
f (m0)
f′ (m0)
, · · · , mk+1 = mk −
f (xm)
f′ (xm)
故相应的迭代函数为
φ(m) = m −
f(m)
f′(m)
设任意 2 个点 m 之间的误差为 ∆ (mn) ，然后确定所需要的精度 ε，当 |∆ (mn)| < ε
时，迭代停止，所得的 φ(m) 即为最终结果. 简化牛顿法是在牛顿迭代法的基础上进行
简化，迭代过程使用一个固定的 f′
(m0) ，而不是每次都计算 f′
(mk) ，简化牛顿法的公
式为
mk+1 = mk −
f (mk)
f′ (m0)
该方法躲避了复杂的计算, 同时也降低了收敛速度.
5.2.2.3 模型求解的结果
表 3 问题 2 的计算结果
| 测线方向夹角/° | 0海里 | 0.3海里 | 0.6海里 | 0.9海里 | 1.2海里 | 1.5海里 | 1.8海里 | 2.1海里 |
|---|---|---|---|---|---|---|---|---|
| 0 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 |
| 45 | 416.36 | 380.28 | 344.83 | 309.67 | 273.85 | 237.54 | 202.11 | 166.09 |
| 90 | 416.69 | 366.31 | 315.75 | 265.98 | 214.39 | 164.22 | 113.63 | 63.35 |
| 135 | 416.37 | 380.84 | 344.47 | 309.91 | 273.39 | 237.42 | 202.13 | 166.43 |
| 180 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 | 415.69 |
| 225 | 416.59 | 451.57 | 487.36 | 523.23 | 558.19 | 594.54 | 630.27 | 665.37 |
| 270 | 416.69 | 467.25 | 517.37 | 568.59 | 618.77 | 669.38 | 719.39 | 770.67 |
| 315 | 416.29 | 451.37 | 487.48 | 523.23 | 558.19 | 594.59 | 630.62 | 665.39 |
由上表分析可知，当测量夹角为 0° 和 180° 时，不论测量船距海域中心点处的距离如何变化，覆盖宽度不会发生改变，均为 415.69m。以测量方向夹角 180° 为分界，当测量船测线方向夹角小于 180°，随着测量船距离海域中心点距离的增加，其覆盖宽度在逐渐降低，最小覆盖范围为 63.35m, 测线方向夹角为 90 ∘ 90^{\circ} 90∘，测量船距离海域中心点处的距离为 2.1 海里；当测量船测线方向夹角大于 180°，随着测量船距离海域中心点处的距离逐渐增加，其覆盖宽度在逐渐增大，覆盖范围最大为 77.67m, 测线方向角为 270 ∘ 270^{\circ} 270∘，测量船距海域中心点处的距离为 2.1 海里。
5.2.3 问题结论
在考虑一个矩形待测海域，换能器开角为 120 ∘ 120^{\circ} 120∘，坡度为 1.5 ∘ 1.5^{\circ} 1.5∘，海水深度为 120m，测量间距为 0.3 海里，约为 500 米的情况下，计算出覆盖宽度，其具体结果如下：当测量方向夹角为 0°，不论测量船距离海域中心点处的距离如何改变，多束波的测深覆盖宽度均为 415.69，不会发生改变。当测量方向夹角为 45 ∘ 45^{\circ} 45∘ 时，测量船距离海域中心点处的距离在 0，0.3，…，1.8，2.1 海里处的覆盖宽度分别为 416.36 ，380.28，344.83 ，309.67 ，273.85，237.54 ，202.11 ，166.09。当测量方向夹角为 90 ∘ 90^{\circ} 90∘ 时，测量船距离海域中心点处的距离在 0，0.3，…，1.8，2.1 海里处的覆盖宽度分别为 416.69，366.31 ，315.75 ，265.98，214.39，164.22 ，113.63 ，63.35 。当测量方向夹角为 135 ∘ 135^{\circ} 135∘ 时，测量船距离海域中心点处的距离在 0，0.3，…，1.8，2.1 海里处的覆盖宽度分别为 416.37 ，380.84 ，344.47 ，
309.91 ，273.39 ，237.42 ，202.13 ，166.43。当测量方向夹角为 180^{\circ} 时，与方向角为 0° 的情况相同，不论测量船距离海域中心点处的距离如何改变，多束波的测深覆盖宽度均
为 415.69，不会发生改变。当测量方向夹角为 270^{\circ} 时，测量船距离海域中心点处的距离在 0，0.3，…，1.8，2.1 海里处的覆盖宽度分别为 416.69 ，467.25 ，517.37 ，568.59 ，618.77 ，669.38 ，719.39 ，770.67。当测量方向夹角为 315^{\circ} 时，测量船距离海域中心点处的距离在 0，0.3，…，1.8，2.1 海里处的覆盖宽度分别为 416.29 ，451.37 ，487.48 ，
523.23 523.23 523.23， 558.19 558.19 558.19， 594.59 594.59 594.59， 630.62 630.62 630.62， 665.39 665.39 665.39。
5.2.4 检验分析
5.2.4.1 灵敏度分析
对海水深度，测线方向角，测线间距进行灵敏度分析，探究覆盖宽度对哪一参数更
为灵敏。为此，当测线距中心点处距离为 0 时，分别取海水深度为 120m，测线方向角为
30°，测线间距为 0m 作为初始值，将覆盖宽度的变化幅度作为指标，来进行探究。在海
水深度为 120m，
β 为 0 时，
深度依次提高 1m，
2m，
3m，
4m，
5m，
经计算得，
海水深度每
次变换后，所得覆盖宽度为 416.19m、419.65m、423.12m、426.59m、430.06m、433.53m。
在测线方向角为 0° 时，角度依次增加 45°，90°，135°，180°，225°，经计算得，测线
方向角每次变换后，所得覆盖宽度为 415.69m、416.19m、416.70m、416.19m、415.69m、
416.19m。
在测线间距为 0m，β 为 45° 时，间距依次增加 100m，200m，300m，400m，500m，经计算得，测线间距每次变换后，所得覆盖宽度为 416.34m、409.79m、403.38m、396.02m、
390.58m、384.08m。
表 4 海水深度对覆盖宽度的影响结果表
| 海水深度/m | 120 | 121 | 122 | 123 | 124 | 125 |
|---|---|---|---|---|---|---|
| 覆盖宽度/m | 416.19 | 419.65 | 423.12 | 426.56 | 430.08 | 433.54 |
根据上述海水深度变化后的覆盖宽度，将所得结果用如下的折线图表示：
图 4 海水深度的变化对覆盖宽度示意图
由上图可知，以海水深度 120m 为起始点，依次增加 1m，观察到覆盖深度在依次增
长且幅度较大。灵敏性数据转换成百分比公式如下：
△y
y
△x
x
dy
dx
·
x
y
根据上述公式计算可知，覆盖宽度对海水深度最为灵敏，灵敏度的值为 1.0085 ，即
改变海水深度时，覆盖宽度变化的幅度最大。
5.2.5 小结
本题建立了三维解析几何模型，完成了特定位置下覆盖宽度的求解。首先，以海域
中心为原点，以平行于海底交线方向为 x 轴，垂直交线方向为 y 轴，铅直方向为 z 轴建
立直角坐标系；其次，求出海底斜坡平面方程和给定位置波束边界的直线方程，进而利
用简化的牛顿迭代法求解波束直线与斜坡平面的交点坐标，再由交点坐标实现覆盖宽度
的计算。计算结果表明：当测线方向夹角大于 180° 时，覆盖宽度逐渐增大，最大覆盖范
围为 770.67m；小于 180° 时，覆盖宽度逐渐减小，最小覆盖范围为 63.35m。其中，测交
方向的夹角为 270° 时，覆盖宽度达到最大值为 770.67m；测线方向夹角为 90° 时，覆盖
宽度最小达到 63.35m。最后，分析了覆盖宽度对换测线间距、测线方向角和海水深度的
灵敏度，结果表明：覆盖宽度对海水深度最为灵敏，灵敏度的值为 1.0085。
5.3 问题三的模型建立与求解
本题建立在问题二的直角坐标系的基础上，作出长宽分别为 2 海里，4 海里的海域
示意图，给出平行侧线的直线方程，并以求出测线累计距离，根据各测线与边界的交点
求出交点处的重叠率，故以测线的累计距离最小为目标，以重叠率满足百分之十到百分
之二十为约束条件，建立单目标优化模型，采用差分进化算法对模型进行求解，得出最
优测线长度。
5.3.1 模型建立与求解
5.3.1.1 模型建立的过程
本题为给出最优的测线规划方案，建立了单目标模型。基于问题二中的空间直角坐
标系，作出南北长 2 海里，东西宽 4 海里的矩形海域的相关分析示意图，如下所示：
图 5 长2海里，宽4海里矩形海域示意图
由上图可知，
以该海域中心为原点 O，
作出平行于测线方向的直线 w，
w+1，
w+2…w+i，
i=1，…，k。过海域平面原点作出一条可以完全覆盖这个待测海域的测线 L。为设计出符
合特定要求的测线，以测量长度 SiRi 之和最短为目标函数，重叠率满足 10% 到 20% 为
约束条件，建立单目标优化模型。
（1）目标函数
本题要求设计一组测量长度最短的测线，即多波束测深系统发射的所有波束在矩形
海域内测量长度之和最短。测量长度可由图 6 中的 SiRi 表示，故目标函数的表达式为：
min
k
∑
i=1
(SiRi)2
其中，SiRi 表示波束测量的长度。
（2）约束条件
本题是在问题二建立的直角坐标系的基础上，故控制 β 的范围为
0 ≤ β ≤ 2π
根据题意，测线需要覆盖海域平面，故南北长 2 海里，东西宽 4 海里的海域平面的
表达式为：
sin α · x + cos α · (z − D) = 0
待测海域测线 L 的直线方程满足以下条件：
x
cos
(
β − π
) =
y
sin
(
β − π
) =
z
沿测线方向的第一条波束的边界射线 Mi
x − xi
sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
沿测线方向的第二条波束的边界射线 Mi+1
x − xi
− sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
任意一条测线相邻条带都必须满足重叠率在 10% 与 20% 之间，即航线 Li 上的每一点都
满足下式子：
10% ≤ η(wi) ≤ 20%, i = 1，…，k
5.3.1.2 模型表达式
综上所述，最终的单目标优化模型如下所示：
min
k
∑
k=1
|SiRi|
s.t.
0 ≤ β ≤ 2π
x
cos
(
β − π
) =
y
sin
(
β − π
) =
z
x2
- y2
= R2
sin α · x + cos α · (z − D) = 0
x − xi
sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
x − xi
− sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
10% ≤ η(wi) ≤ 20%, i = 1…k
5.3.2 模型求解
5.3.2.1 模型求解原理
在自然界中，存在着遗传、变异和选择等自然机制，这些机制导致生物体之间进行
竞争和适应，
从而使得更适应环境的生物能够生存下来，
而不适应环境的生物会被淘汰。
这种优胜劣汰的过程推动了生物的进化，使得生物从低级向高级逐渐演化。差分进化算
法就是从这种模式中产生的一种智能优化算法。差分进化算法是基于群体的理论的优化
算法，与进化算法相比，保留了基于种群的全局搜索策略，采用实数编码，并结合差分
的简单变异操作和” 一对一” 竞争生存策略，差分进化算法降低了操作的复杂性。此外，
差分进化算法独特的记忆能力使其能够动态地追踪当前的搜索情况，并相应地调整搜索
策略。这使得差分进化算法具有较强的收敛能力和稳健性，并且无需依赖于问题的特定
信息。
5.3.2.2 模型求解的过程
步骤一：初始化数据。随机初始化数目为 NP (NP = 50) 的 D(D=10) 维参数向量 x，
x (i) 表示第 i 个解，
每个解参数可以表示为 x (i, j) ，
i = 1, 2, …, NP, j = 1, 2, …, D, (NP =
50; D = 10; )
步骤二：变异。对于每个解向量 x (i, j)，对应的变异向量 v 可以表示为 v (i) =
x (r0) + F ∗ (x (r1) − x (r2))，其中，r0,r1,r2 为属于 1，…，NP 的三个随机
数，并且 i,r0,r1,r2 都不相同。变异算子 F 取值范围为 [0,2]，F 过小可能陷入局部最优，F
过大则不容易收敛。若出现边界问题，即如果变异以后的值 v (i, j) 超出了边界，再随机
再选择一个数，或者直接去边界值。
步骤三：交叉。求取交叉向量 u，对于 u 的每个值，随机产生一个值，有：
rand () ≤ CR, u (i, j) = v (i, j) rand () > CR, u (i, j) = v (i, j)
其中，CR = 0.4 是交叉算子，rand () 是一个范围是 [0,1] 的随机数，用来控制选择变异
向量值还是原来的向量值。
步骤四：选择。把交叉向量和原向量对比，选择较优的那个，这里交叉向量之和对
应的原向量对比，也就是对比 u (i)（新解）和 x (i) （原来的解）哪个更优，就选择哪个
作为新的解向量，更新向量 x，进行下一步。
步骤五：终结条件。当最后的解满足条件，或者遍历次数达到最大，则结束，否则
重复步骤二到四，直到找到最优解。
5.3.2.3 模型求解的结果
由算法计算可知，方向角为 0° 时，测量长度为 1.34 × 105
m，为最短的测量长度，
符合题目的设计要求，因此，方向角为 0° 时是最佳方向角，沿测线方向；方向角为 167°
时，测量长度为 1.71 × 105
m, 此时，测量长度最大，不符合题目的设计要求。
由差分进化算法得到变化间距的示意图，如下：
图 6 变化间距的示意图
由图可知，测线的间距为变化间距，测线间距序号取 0 时，间隔距离最大为 519m；
测线间距序号取 38 时，间隔距离最小为 47m。根据上述分析，及差分进化算法的求解，
求得共有 39 条测线，测量长度最短为 144.457km。
5.3.3 问题结论
本题要求设计出最优的测线规划方案。以求出的测线累计距离最小为目标，以特定
的海域中心处的海水深度，坡度，多波束换能器的开角，相邻条带满足重叠率为 10% 到
20% 为约束条件，建立优化模型。由差分进化算法可求解出方向角为方向角为 0° 时，测
量长度为 1.34 × 105
m，为最短的测量长度，符合题目的设计要求，因此，方向角为 0°
时是最佳方向角，沿测线方向；方向角为 167° 时，测量长度为 1.71 × 105
m, 此时，测量
长度最大，不符合题目的设计要求。由于测线的间距为变化间距，故测线间距序号取 0
时，间隔距离最大为 519m，测线间距序号取 38 时，间隔距离最小为 47m。在总共 39 条
测线中，测量长度最短为 144.457km。
5.3.4 检验分析
5.3.4.1 方法比较
追赶法的数学原理如下：
追的过程
uk =
rk − Kk−1 − ak
bk − vk + ak
, k = 1, 2, · · · , n
vk =
bk − Vk−1ak
, k = 1, 2, · · · , n
赶的过程
{
xn = un
xk = uk − vkxk+1, k = n − 1, n − 2, · · · , 2, 1
步骤一：输入已知数据 ai, bi, ci, fi
步骤二：y1 =
f1
b1
; d = b1
步骤三：i 从 2 到 n 循环，得到 βi−1 =
ci−1
d
；βi−1 =
ci−1
d
;yi =
(fi − aiyi−1)
d
步骤四：进而求的 xn = yn，否则，返回第三步重新计算
步骤五：i 从（n-1）循环到 1 循环
步骤六：进而求的 xi = yi − βixi+1，否则，返回第五步重新计算
步骤七：输出 x1, x2, …, xn
经过以上的对函数的重新计算，得到共有 39 条测线，测量长度最短为 144.476km，
得到的结果误差为 0.000131，小于 0.001，证明了计算结果的准确性。
5.3.5 小结
本题建立了单目标优化模型，给出了最优的测线规划方案。首先，在问题二直角
坐标系的基础上，给出平行测线的直线方程，并求出测线累计距离，以该距离最小为
目标；其次，求出各测线与边界的交点，并求得这些交点处的重叠率，以重叠率满足
10% ∼ 20% 为约束条件，建立单目标优化模型，利用差分进化算法进行求解。计算结果
表明，最优的测线总长度为 14457m。最后，利用追赶法对模型的计算结果进行重新求
解，计算结果为 144476m，与差分进化算法求得的结果误差小于 1%，验证了计算的结
果可信性。
5.4 问题四的模型建立与求解
本题要求给出最优的多波束测线方案，建立了多目标优化模型。首先，对附件中的
深海数据进行拟合，得到拟合的曲面方程；其次，利用航向角求出累计测线距离，并以
该距离最小为目标 1。由航线求出航线在拟合曲面上的投影方程，并利用波束边界直线
方程确定覆盖宽度范围，进而实现计算重叠率大小的计算，以重叠率大于 20% 的比例
最小为目标 2，建立多目标优化模型。
5.4.1 模型建立与求解
5.4.1.1 模型建立的过程
本题为给出最优的测线方案，建立了多目标优化模型，基于问题二中的空间直角坐
标系，作出南北长 5 海里，东西宽 4 海里的矩形海域相关分析示意图，如下所示：
图 7 长5海里，宽4海里矩形海域示意图
由上图可知，
以该海域中心为原点 O，
作出平行于测线方向的曲线 v，
v+1，
v+2…v+i，
i=1，…，k。
（1）目标函数
为设计出符合特定要求的测线，需从测线总长度与重叠率这两个方面考虑, 既要满
足测线总长度尽可能最短，也要实现相邻条带之间的重叠率在 20% 以下。因此，本题以
测量总长度最短为第一个优化目标函数，相邻条带间的重叠率最小为第二个优化目标，
在该海域平面上重叠率满足 10% 到 20% 等为约束条件，建立多目标优化模型。
min
∑ ∫ Ri
Si
√
x′ (t)2
- y′ (t)2
dt
min
∑ ∫
n>20%
√
x′ (t)2 - y′ (t)2
dt
∑ ∫ Ri
Si
√
x′ (t)2 - y′ (t)2
dt
在等权重下将多目标转为单目标，表达式如下
min 0.5 ×
∑ ∫ Ri
Si
√
x′ (t)2 - y′ (t)2
dt + 0.5 ×
∑ ∫
n>20%
√
x′ (t)2 - y′ (t)2
dt
∑ ∫ Ri
Si
√
x′ (t)2 - y′ (t)2
dt
（2）约束条件
本题是在问题二建立的直角坐标系的基础上，故控制 β 的范围为
0 ≤ β ≤ 2π
根据题意，测线需要覆盖海域平面，故南北长 2 海里，东西宽 4 海里的海域平面的
表达式为：
sin α · x + cos α · (z − D) = 0
待测海域测线 L 的直线方程满足以下条件：
x
cos
(
β − π
) =
y
sin
(
β − π
) =
z
通过对数据的拟合，得到曲面的表达式：
24.4 − 14.4x − 4y + 14 · 4x2
- 3 − 2xy + 3.2y2
− 2 · 105 × 10−6
X3
−3 − 2x2
y2
− 9.708 × 10−6
× y2 - 7.417 × 10−0
y3
= 0
沿测线方向的第一条波束的边界射线 Mi
x − xi
sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
沿测线方向的第二条波束的边界射线 Mi+1
x − xi
− sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
任意一条测线相邻条带都必须满足重叠率在 10% 与 20% 之间，即：
10% ≤ η(Mi) ≤ 20%, i = 1，…，k
5.4.1.2 模型表达式
综上所述，最终的多目标优化模型的表达式为如下所示：
min 0.5 ×
∑ ∫ Ri
Si
√
x′ (t)2
- y′ (t)2
dt + 0.5 ×
∑ ∫
n>20%
√
x′ (t)2 - y′ (t)2
dt
∑ ∫ Ri
Si
√
x′ (t)2 - y′ (t)2
dt
s.t.
0 ≤ β ≤ 2π
x
cos
(
β − π
) =
y
sin
(
β − π
) =
z
x2 - y2
= R2
sin α · x + cos α · (z − D) = 0
x − xi
sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
x − xi
− sin θ
cos
(3π
− β
) =
y − yi
sin θ
sin
(3π
− β
) =
z − zi
cos θ
10% ≤ η(Mi) ≤ 20%
5.4.2 模型求解
5.4.2.1 模型求解原理
混合模拟退火禁忌搜索算法是一种组合了模拟退火算法和禁忌搜索算法的元启发
式优化算法。它被广泛应用于解决组合优化问题，该混合算法吸收了模拟退火算法简单
通用，具有逃离局部最优的陷阱的能力，效率不依赖初始解的优点和禁忌搜索算法灵活
跳出局部最优解，方法简单通用的优点，使得该混合算法既依赖于初始解，又有记忆功
能。其基本思想是：在进入模拟退火算法的之前，利用禁忌表构造一个候选邻域，从候
选邻域中选取最好的解，若这个解比当前最好的解还好，就接收这个解，否则，就进入
模拟退火算法的每次接收新的解都要更新禁忌表。
5.4.2.2 模型求解的过程
该算法步骤如下：
步骤一：设定初始温度 T0=100 ，T= T0；
步骤二：随机产生初始解 X0 ;x=X0；X*=X0；
（X* 为最优解）
步骤三：利用禁忌表构造候选邻域 N ∗ (x) ；
步骤四：在候选邻域中选择最好的解 x
′
，如果 f
(
x
′
)
< f （x∗
）, 则接收 x
′
，x = x
′
，更新禁忌表，并转到步骤 7；步骤五：重复下列步骤 L 次在邻域中随机产生解 x
′
；
df = f
(
x
′
)
− f (x) 如果 df < 0, 则 x = x
′
; 否则以概率 exp (−df/T) （接受解的概率）
接收更新禁忌表；如果 f (x) < f (x∗
), 则更新 x∗
= x；
步骤六：如果满足终止条件，则输出 x∗
；否则退火 T = aT 并返回步骤三。
步骤七：如果满足终止条件，则输出 x∗
；否则返回步骤三。
5.4.2.3 模型求解的结果
本题采用二元三次幂函数对附件中深海数据进行拟合，结果如下图所示：
图 8 深海数据拟合图
依据混合模拟退火禁忌搜索算法得到测量船航行轨迹图，如下：
图 9 测量船航行轨迹图
由上图可知测量船沿等深线行驶，测线间的间距不同, 由算法计算得出测线间的间
距最大值为 587.62m, 测线间的间距最小值为 32.15m。测量船的航行轨迹中间较为密集，
两端较为稀疏。
根据混合模拟退火禁忌搜索算法求解及上述分析得到测线的总长度为 331.4729km，
漏测区域占总待测海域面积的百分比为 4.4769%，在重叠区域中，重叠率超过 20% 部分
的总长度为 7.1772km。
5.4.3 问题结论
本题要求给出最优的多波束测线方案，建立了多目标优化模型。首先，对附件中的
深海数据进行拟合，得到拟合的曲面方程；其次，利用航向角求出累计测线距离，并以
该距离最小为目标 1。由航线求出航线在拟合曲面上的投影方程，并利用波束边界直线
方程确定覆盖宽度范围，进而实现计算重叠率大小的计算，以重叠率大于 20% 的比例
最小为目标 2，建立多目标优化模型, 根据混合模拟退火禁忌搜索算法求解及上述分析
得到测线的总长度为 331.4729km，漏测区域占总待测海域面积的百分比为 4.4769%，在
重叠区域中，重叠率超过 20% 部分的总长度为 7.1772km。
5.4.4 检验分析
5.4.4.1 仿真检验
灰狼优化算法是通过模拟自然界中灰狼群体的社会等级机制和捕食行为而提出的
一种新型群体搜索方法。灰狼群体的社会等级分为四层：首领狼 (头狼)α ，历史次最优
解记为副首领狼 β，
历史第三最优解记为普通狼 δ，
以及种群其他个体记为底层狼 ω 。
灰
狼群体在捕获猎物时，其他灰狼个体在头狼 α 的带领下有组织地对猎物进行围攻。第 i
只灰狼的位置记为 Xi = (xi1, , xi2, · · · , xid) ，猎物的位置对应于优化问题的全局最优解。
狼群猎食行动包括以下两个主要步骤:
(1) 包围猎物：描述灰狼逐渐接近并包围猎物的行为，满足如下公式：
D = |CXp (t) − Xt| ,
X (t + 1) = Xp (t) − A · D,
在上式中，t 表示当前迭代次数，X (t) 表示第 t 代灰狼个体的位置，Xp 为第 t 代猎物位
置，常数 C 为摆动因子，常数 A 为收敛因子，定义如下：
C = 2 · rand2,
A = 2a · rand1 − a。
其中 rand1,rand2 表示 [0, 1] 之间的随机变量；变量 a 称为随机因子，它随着迭代次数的
增大从 2 减小到 0，数学表达为：
a = 2 −
t
tmax
。
（2) 捕食狩猎：当灰狼判断出猎物所在位置时，通常由头狼 a 引领 6 和 o 狼发动捕食
行为。因此狼群可以根据 α,β,δ 三者的位置判断出猎物所在方向，进而更新的灰狼位置：
X1 = Xα − A1 · Dα,
X2 = Xβ − A2 · Dβ,
X3 = Xδ − A3 · Dδ,
Xp(t + 1) =
X1 + X2 + X3
。
X (t) 表示当前灰狼位置向量。
综上所述，灰狼算法流程可归结如下：
步骤 1: 初始化灰狼算法基本参数：NP 为种群个数，初始化参数 a,A 和 C，问题维
数 d 和最大迭代次数 Tmax 等参数；
步骤 2: 计算灰狼个体的目标函数值并进行排序，得到当前的历史最优解 Xα，次最
优解入 Xβ，第三最优解 Xδ 和 ω 普通狼。
步骤 3: 由式 (32) 计算群体中其他灰狼个体分别与 α,β, δ 狼的距离，并根据式上式
更新每个灰狼个体的位置。重新计算所有个体狼的适应度值并进行排序，更新最优值前
三的灰狼个体 Xα，Xβ，Xδ 的位置。
步骤 4: 对普通狼 ω 实施差分单纯形扰动策略，重新计算前三的交狼个体 Xα，Xβ，
Xδ 的位置。
步骤 5: 更新 a、A 和 C 等参数的值，算法结束。
经过以上的对函数的重新计算，测线的总长度为 335.8165km，漏测海区占总待测海
域面积的百分比为 4.59%，在重叠区域中，重叠率超过 20% 部分的总长度为 7.21km。得
到的结果误差为误差率小于 3%，证明了计算结果的准确性。
5.4.5 小结
本题建立了多目标优化模型，给出了最优的多波束测线方案。利用二元三次幂函数
对附件中的海深数据进行拟合，得到拟合的曲面方程，其次，以该累计测线距离距离最
小为目标 1, 以重叠率大于 20% 的比例最小为目标 2，建立多目标优化模型。量纲归一
化后，利用加权法将多目标转为单目标优化模型，利用混合模拟退火禁忌搜索算法进行
求解。计算结果表明：测线总长度为 331.47km，漏测海区占总海域面积百分比的 4.48%，
重叠率超过 20% 的总长度为 7.18km，占比约为 2.16%。最后，采用灰狼优化算法对模型
的求解结果进行检验，结果为：测线的总长度为 331.4769km，漏测区域占总待测海域面
积的百分比为 4.47%，在重叠区域中，重叠率超过 20% 部分的总长度为 7.16km。得到的
结果误差为 0.000163，小于 0.001，证明了计算结果的准确性。
六、模型评价与推广
6.1 模型评价
（1）本文一二问建立了几何模型，解决了特定位置下的指标计算和覆盖宽度问题。
几何模型具有精度高，可视化，方便修改和优化的特点，可以非常精确地描述测量船的
相关位置信息，
在设计工程和制造领域，
高精度的几何模型可以保证产品的质量和性能。
几何模型可以以图形形式展现出来，使得人们可以直观地了解和理解物体的形状、结构
和特点。这使得几何模型成为教育、科学研究和娱乐等领域的有用工具。在几何模型中，
人们可以非常方便地修改、删除或添加元素，以达到优化设计的目的。这使得几何模型
可以快速有效地进行设计和优化。
(2) 本文三四问建立了优化模型，给出了最优的测线规划方案。采用差分进化算法
和混合模拟算法对模型进行求解，既吸收了模拟退火算法简单通用，具有逃离局部最优
的陷阱的能力，效率不依赖初始解的优点和禁忌搜索算法灵活跳出局部最优解，方法简
单通用的优点。
6.2 模型推广
（1）几何模型可以通过计算机程序进行自动化处理和分析，例如计算物体的体积、
重心、惯性矩等参数，以及进行模拟分析和预测等操作。这使得几何模型成为计算机辅
助设计和工程领域的重要工具。
（2）通常会采用多波束测深技术来进行海洋测绘勘探工作, 多波束测深系统已成为
当前主流的测深设备，在水下地形测量、海底地貌勘察和水下目标检测等领域得到广泛
应用。
参考文献
[1] 彭斌, 刘慧鑫, 陶耀辉. 基于变径基圆渐开线涡旋压缩机的几何模型及优化研究 [J].
上海交通大学学报,2023,57(08):1046-1054.
[2] 万新儒, 刘单, 邵尉哲等. 高斯消元法计算技巧的研究及应用 [J]. 电力系统及其自动
化学报,2018,30(04):109-113.
[3] 尚勇敏, 宓泽锋, 周灿等. 中国城际低碳技术转移对碳排放的影响——基于知识学习
与技术学习“二分法”视角 [J]. 资源科学,2023,45(04):827-842.
[4] 赵化时, 宋智强, 黄耀辉等. 含多类型直流的交直流系统潮流统一迭代算法研究 [J].
华北电力大学学报 (自然科学版):2023,9:1-10
[5] 冯成, 龚晓峰, 雒瑞森. 基于 FFT 和黄金分割的快速 DOA 估计 [J]. 计算机仿
真,2021,38(09):159-163.
[6] 贾传果, 赵进级, 李兴等. 内嵌牛顿迭代的线性隐式算法及其在结构非线性动力分析
中的应用 [J]. 振动与冲击,2023,42(09):177-188.
[7] 丁永刚, 宋战炯, 陈科委等. 基于蒙特卡洛法的粮食侧压力下 SIW 墙板可靠性及灵敏
度分析 [J]. 河南工业大学学报 (自然科学版),2023,44(04):106-113.
[8] 毛西耶子, 孙若辰, 段青云. 单目标和多目标优化在 SWAT 模型率定中的对比 [J]. 南
水北调与水利科技 (中英文),2023,21(02):289-300.
[9] 韩宝辉, 赵起超, 常荣等. 叶绿素 a 浓度反演模型：堆栈自编码器粒子群优化 BP 神
经网络 [J]. 地球信息科学学报,2023,25(09):1882-1893.
[10] 王子实, 吴耀华. 基于差分进化算法的 FMS 中机器与多载 AGV 调度 [J]. 控制与决
策,2023,9(10):1-9.
[11] 张晓波, 闪丽洁, 张瑶兰. 追赶法在含闸分洪河道水动力计算中的改进研究 [J]. 中国
水利水电科学研究院学报 (中英文),2023,21(01):10-22.
[12] Wei G,Xiaoxu D,Chen Z, et al. HSAEP: a new platform to evaluate hybrid simulation
algorithms[J]. Frontiers of Computer Science,2021,15(5):1-15.
[13] 王海群, 邓金铭, 张怡等. 基于改进混合灰狼优化算法的无人机三维路径规划 [J]. 无
线电工程,2023,9:1-12
附 录
附录一 问题一的代码
下面是：问题一的代码：
main1.m 代码如下
% 清空工作区
clc, clear;
D = -70;% 设置参数
DI = 200;% 设置水深、间隔
% 设置坡度、换能器开角度、航向
jiaodu = [1.5, 90, 120]*pi/180;% 角度转弧度
%% 计算部分
RS = -800:DI:800;%每一个数据
for i=1:length(RS)%遍历矩阵
[rq1,rq2,HI] = solvepoint(RS(i),D,jiaodu);%调用函数来解决坐标
aiD(i,1)=double(rq1.x);%将算的的坐标存起来
aiD(i,2)=double(rq1.y);%y的坐标
aiD(i,3)=double(rq1.z);%z的坐标
biD(i,1)=double(rq2.x);%第二个x的坐标
biD(i,2)=double(rq2.y);%第二个y的坐标
biD(i,3)=double(rq2.z);%第三个z的坐标
deepi(i) = double(HI);%计算海水深度
end
[fugai1,chongdie1] = wideandetafun(aiD,biD);
%% 显示
disp('测线距中心点处的距离/m')
disp(RS);
disp('海水深度/m')
disp(deepi*(-1));
disp('覆盖宽度/m')
disp(fugai1');
disp('与前一条测线的重叠率/%')
disp(chongdie1);
以下是问题 1 和问题中附带的函数：
这个函数是求方程交点坐标：
solvepoint.m 代码如下
function [r1,r2,hi] = solvepoint(ri,H,abt)
syms x y z;%定义符号变量
if ri == 0%当船在远点的位置时的特殊情况
xi=0;
yi=0;
zi=0;
else
eqfun01 =
sin(abt(2)-0.5*pi)*(x-0)-cos(abt(2)-0.5*pi)*(y-0);%列出方程1
eqfun00 = x*x+y*y-ri*ri;%列出方程2
r0 = solve(eqfun01,eqfun00,x,y);%利用solve函数解方程
if ri > 0%判断当ri大于0
if cos(abt(2)-0.5*pi)*(r0.x(1)-0)+sin(abt(2)-0.5*pi)
*(r0.y(1)-0)<0
xi = r0.x(1);%x坐标
yi = r0.y(1);%y的坐标
zi = 0;
else
xi = r0.x(2);%同上来计算
yi = r0.y(2);%同上再次计算y的坐标
zi = 0;
end
else
if cos(abt(2)-0.5*pi)*(r0.x(1)-0)+sin(abt(2)-0.5*pi)
*(r0.y(1)-0)<0
xi = r0.x(2);%x(2)赋值到xi上
yi = r0.y(2);
zi = 0;
else
xi = r0.x(1);
yi = r0.y(1);
zi = 0;
end
end
end
eqfun1 = sin(abt(1))*(x-0)+0*(y-0)+cos(abt(1))*(z-H);
%坡度的方程
eqfun21 =
sin(abt(3)*0.5)*sin(1.5*pi-abt(2))*(x-xi)-sin(abt(3)*0.5)
*cos(1.5*pi-abt(2))*(y-yi);
eqfun22 = cos(abt(3)*0.5)*(x-xi)-sin(abt(3)*0.5)
*cos(1.5*pi-abt(2))*(z-zi);
%两个光束的方程
eqfun31 = -sin(abt(3)*0.5)*sin(1.5*pi-abt(2))*(x-i)+
sin(abt(3)*0.5)*cos(1.5*pi-abt(2))*(y-yi);
eqfun32 =
cos(abt(3)*0.5)*(x-xi)+sin(abt(3)*0.5)*cos(1.5*pi-abt(2))
*(z-zi);
%测线的方程
r1 = solve(eqfun1,eqfun21,eqfun22,x,y,z);
r2 = solve(eqfun1,eqfun31,eqfun32,x,y,z);
%利用solve来计算两个方程的坐标
eq5 = sin(abt(1))*(xi-0)+0*(yi-0)+cos(abt(1))*(z-H);
hi = solve(eq5,z);%最后再计算高度
end
这个函数是利用坐标来求面积
wideandetafun.m 代码如下
function [AABBi,chong] = wideandetafun(AXi,BYi)
%这个函数用来求解坐标来计算面积
temp = AXi-BYi;%表示差值
AABBi = sqrt(sum(temp.^2,2));%求根号的比
chong(1)=nan;%第一个为无
for i=2:length(AXi)%遍历矩阵
wi = AABBi(i-1);%最后再一次求wi
di = sqrt(sum((AXi(i,:)-AXi(i-1,:)).^2,2));%求一次di
chong(i) = (wi-di)*100/wi;%最后求重叠率
end
问题一的检验分析的图：
tu1.m 代码如下
x =[120 120.2 120.4 120.6 120.8 121.0];
y = [19.7822 20.0912 20.3966 20.7073 21.0145 21.321];
figure();
plot(x,y,'o');
hold on;
plot(x,y,'k-','LineWidth',1);
set(gcf,'Units','centimeters','Position',[5 5 20 12]);
%legend('角度','覆盖率','Location', 'southeast',FontSize=16);
set(gca,'FontName','Times New Rome','FontSize',16);
xlabel('换能器开角/°','FontName','宋体','FontSize',18);
ylabel('覆盖率/%','FontName','宋体','FontSize',18);
附录二 问题二的代码
main2.m 代码如下
% 清空工作区
clc, clear;
Hi = -120;
% 设置参数
DD = 0.3*1852;
% 设置水深、间隔
%计算部分
for betaijiaodu = 0:45:315%每次增加45度改变方向
% 设置坡度、换能器开角度、航向
jiaoduhe = [1.5, betaijiaodu,
120]*pi/180;%代入条件坡度、角度、转换器的开角
ri = 0:DD:2.1*1852;%初始化ri
for i=1:length(ri)%遍历数组
[R11,R22,HHIi] =
solvepoint(ri(i),Hi,jiaoduhe);%计算每个点的坐标
%计算每个点的坐标，然后在计算覆盖面积
aiI(i,1)=double(R11.x);
aiI(i,2)=double(R11.y);
aiI(i,3)=double(R11.z);
biI(i,1)=double(R22.x);
biI(i,2)=double(R22.y);
biI(i,3)=double(R22.z);
deepi(i) = double(HHIi);
end
[widedu,eta] = wideandetafun(aiI,biI);%根据坐标计算面积
% 显示结果
if betaijiaodu==0
disp('测量船距海域中心点处的距离/海里')
disp(0:0.3:2.1)% 每次增加0.3海里
end
% 开始输出结果
strtemp = sprintf('覆盖宽度/m,此时的 betaijiaodu = %d
',betaijiaodu);
disp(strtemp);
disp(widedu');
end
第二问的检验分析的图：
tu2.m 代码如下
x =[120 121 122 123 124 125];
y = [416.1915 419.6598 423.1280 426.5963 430.0646 433.5328];
figure();
plot(x,y,'o');
hold on;
plot(x,y,'k-','LineWidth',1);
set(gcf,'Units','centimeters','Position',[5 5 20 12]);
%legend('角度','覆盖率','Location', 'southeast',FontSize=16);
set(gca,'FontName','Times New Rome','FontSize',16);
xlabel('海水深度/m','FontName','宋体','FontSize',18);
ylabel('覆盖宽度/m','FontName','宋体','FontSize',18);
附录三 问题三的代码
此代码运行时间长
main31.m 代码如下
%% 清空工作区%%
%clc；
%clear;
%% 设置参数
% 设置水深、间隔
DD = -110;%海水的深度
WXX = 4*1852;%海域的面积计算
WYY = 2*1852;%海域的宽度计算
%% 计算数据的方法
disp("p=");
kaijiaodu = 0:60:360;%需要计算的角度。
for sss1 = 1:size(kaijiaodu,2)%遍历数组之后
betai = kaijiaodu(sss1);%计算气体中的值
dqun = 1;%将dqun赋值为1
dQ0 = 0.3;%这是他的间距
dQ1 = 0.05;%这是间隙
while dqun==1%下面进行循环
d01=(dQ0+dQ1)*0.5;
disp(d01);%来输出do1
di=d01*1852;%计算他的宽度
eee1 = [1.5, betai, 120]*pi/180;%初始化数据
[gx0,gy0,gx1,gy1] =
gridfun(di,WXX,WYY,eee1);%调用gridfun函数之后计算处结果。
dmin = (gx0+0.5*WYY).^2+gy0.^2;%求出测试线的最小值
[ixw,iyw] =
find(dmin==min(min(dmin)));%利用find的函数寻找范围。
dmin = (gx0-0.5*WYY).^2+gy0.^2;%找到了最小值最后求出来
[ixe,iye] = find(dmin==min(min(dmin)));%同理如上
xnode = [ixw(1),ixe(1)];%最后算出坐标
ynode = [iyw(1),iye(1)];
for li=1:2%遍历求出li的整个
i = xnode(li);%从xnode中取出i来
for j=1:size(gx0,2)%遍历gx0数列
xyz = [gx0(i,j),gy0(i,j),0];%jjj存储数据
[r1,r2,hi] =
solvegrid(xyz,DD,eee1);%调用函数solvegrid
AAXi(j,1)=double(r1.x);%赋值用浮点型写入新的矩阵AAIj1
AAXi(j,2)=double(r1.y);%赋值用浮点型写入新的矩阵AAIj2
AAXi(j,3)=double(r1.z);%赋值用浮点型写入新的矩阵AAIJ3
BBJi(j,1)=double(r2.x);%赋值用浮点型写入新的矩阵AAIJ1
BBJi(j,2)=double(r2.y);%赋值用浮点型写入新的矩阵AAIJ2
BBJi(j,3)=double(r2.z);%赋值用浮点型写入新的矩阵AAIJ3
deepi(j) =
double(hi);%赋值用浮点型写入新的矩阵,赋值带入海水深度
end
[wide,eta] =
wideandetafun(AAXi,BBJi);%调用函数来求解坐标的计算
geta(li) = eta(ynode(li));%用求差方法来测量长度
end
if min(geta)>10
dQ1=d01;
else
dQ0=d01;
end
if abs(dQ1-dQ0)/dQ1 <0.02
dqun=-1;
end
end
lineadd = sumli(di,WXX,WYY,eee1);
liarray(sss1) = lineadd(1) +sum(lineadd(2:end))*2;
end
plot(kaijiaodu,liarray);%换图表示出结果来。
main32.m 代码如下
% 清空工作区
clc, clear;
% 设置参数
% 设置水深、间隔
D = -110;%海水深度
wxX = 4*1852;%海域的宽度
wyY = 2*1852;%海域的长度
%% 计算
jiaodu = 90;
huduzhi = [1.5, jiaodu, 120]*pi/180;%角度转化弧度制。
syms x y;%定义符号变量
x1=-0.5*wxX;%计算x1
y1=tan(huduzhi(1))*x1+D;%计算y1
li = 1; %初始化li为 1
while li>0 %遍历数组；
eq01 = tan(huduzhi(1))*x-y+D;%方程1的求解
eqil = cos(0.5*huduzhi(3))*(x-x1)-sin(0.5*huduzhi(3))*(y-y1);
%方程2的求解
eqi0 = cos(0.5*huduzhi(3))*(x-x1)-sin(0.5*huduzhi(3))*(0-y1);
%方程3的求解
ril = solve(eqi0,x);%利用函数的求解方程
xi0 = double(ril);%存放x1的值
xi0all(li) = xi0;%去除他的值
eqir = -cos(0.5*huduzhi(3))*(x-xi0)-sin(0.5*huduzhi(3))*(y-0);
%函数值
rir = solve(eq01,eqir,x,y);
%函数求解得出
x2 = rir.x;%存出其中的函数x
y2 = rir.y;%存出其中的函数y
x1 = x1+(x2-x1)*0.8;%存出其中的函数坐标
y1 = y1+(y2-y1)*0.8;%存出其中的函数坐标值y
if x1>=0.5*wxX
disp("*****");
disp(li);
disp('总里程/m');
disp(li*wyY);
disp(li*wyY<1.33*10^5);
li = -1;
else
li=li+1;
end
end
plot([-0.5*wxX,0.5*wxX],[-0.5*wyY,-0.5*wyY]);
hold on
plot([-0.5*wxX,0.5*wxX],[0.5*wyY,0.5*wyY]);
hold on
plot([0.5*wxX,0.5*wxX],[-0.5*wyY,0.5*wyY]);
hold on
plot([-0.5*wxX,-0.5*wxX],[-0.5*wyY,0.5*wyY]);
for i=1:size(xi0all,2)
plot([xi0all(i),xi0all(i),],[-0.5*wyY,0.5*wyY]);
hold on;
end
diff(xi0all);
pp.m 代码如下
clc
clear
RSS=-800:200:800;
HHH=tan(1.5*pi/180)*(-RSS);
for QQ=60:65
D=HHH+70;
d=200;
l1=sin(QQ*pi/180)*D/sin(31.5*pi/180);
l2=sin(QQ*pi/180)*D/sin(28.5*pi/180);
w1=(l2+l1);
g=1-d./w1;
g(1)=NaN;
plot(RSS,g);
hold on
end
以下是第三小题附带的函数：
gridfun.m 代码如下
function [x0,y0,x1,y1] = gridfun(di,wx,wy,abt)%用来求解坐标的函数
rd = sqrt(wx*wx+wy*wy);%开根号求解长度
num = floor(0.5*rd/di)+1;%标记位置，数量
x = (-num:num)*di;%用数列来表示x
y = (num:-1:-num)*di;%计算出来y的值。
for i =1:size(x,2)%便利x矩阵
for j=1:size(y,2)%便利y的矩阵，然后如赋值
x1(j,i)= x(i);
y1(j,i) = y(j);
end
end
A = [cos(pi-abt(2)) -sin(pi-abt(2));%函数计算
sin(pi-abt(2)) cos(pi-abt(2))];%函数计算的方程
for i =1:size(x1,1)%遍历x1的函数
for j =1:size(x1,2)%计算函数的表
xij =x1(i,j);%赋值求解过程。
yij =y1(i,j);%赋值缘散
temp= A * [xij;yij];%矩阵个位相乘
x0(i,j)=temp(1);
y0(i,j)=temp(2);
end
end
end
slovegrid.m 代码如下
function [r1,r2,hi] = solvegrid(xyz,H,at)%计算方程的函数
syms x y z;%定义成字符变量
xiX= xyz(1);%赋值xI
yiY= xyz(2);%赋值yI
ziZ =xyz(3);%赋值zI
eqfun1 = sin(at(1))*(x-0)+0*(y-0)+cos(at(1))*(z-H);
%方程1
eqfun21 =
sin(at(3)*0.5)*sin(1.5*pi-at(2))*(x-xiX)-sin(at(3)*0.5)
*cos(1.5*pi-at(2))*(y-yiY);
eqfun22 =
cos(at(3)*0.5)*(x-xiX)-sin(at(3)*0.5)*cos(1.5*pi-at(2))
*(z-ziZ);
%方程2用代码表示
eqfun31 =
-sin(at(3)*0.5)*sin(1.5*pi-at(2))*(x-xiX)+sin(at(3)*0.5)
*cos(1.5*pi-at(2))*(y-yiY);
%方程3用代码表示
eqfun32 =
cos(at(3)*0.5)*(x-xiX)+sin(at(3)*0.5)*cos(1.5*pi-at(2))
*(z-ziZ);
%方程4用代码表示
r1 = solve(eqfun1,eqfun21,eqfun22,x,y,z);
r2 = solve(eqfun1,eqfun31,eqfun32,x,y,z);
%利用solve来求解函数算的结果
eq5 = sin(at(1))*(xiX-0)+0*(yiY-0)+cos(at(1))*(z-H);
hi = solve(eq5,z);
end
sumli.m 代码如下
function liadd = sumli(di,wWx,wYy,abFFt)%最后求得的结果
if abFFt(2) == 0%在特殊的位置，为0的时候求解
liadd = floor(wWx/di)*wYy;%最后算的结果
elseif abFFt(2) == 0.5*pi%一般情况
liadd = floor(wYy/di)*wWx;%得到的数据
else
li = 0;%赋值li为0
while li>=0%遍历数组
syms x y;%定义字符常量
eqfun1=
-tan(abFFt(2))*x-y+(li-1)*di/sin(abFFt(2)-0.5*pi);
%函数1的解
eqfun11= y+wYy*0.5;%函数二
eqfun12= y-wYy*0.5;%函数3
eqfun21= x+wWx*0.5;%函数4
eqfun22= x-wWx*0.5;%函数5
r1= solve(eqfun1,eqfun11,x,y);%solve来写出函数1的解值。
r2= solve(eqfun1,eqfun12,x,y);%solve来写出函数2的解值。
r3= solve(eqfun1,eqfun21,x,y);%solve来写出函数3的解值。
r4= solve(eqfun1,eqfun22,x,y);%solve来写出函数4的解值。
xZy(1,:) = [double(r1.x),double(r1.y)];
xZy(2,:) = [double(r2.x),double(r2.y)];
xZy(3,:) = [double(r3.x),double(r3.y)];
xZy(4,:) = [double(r4.x),double(r4.y)];
j=1;
for i=1:4%遍历数组之后
if abs(xZy(i,1))<=0.5*wWx & abs(xZy(i,2))<=0.5*wYy
temp(j,:) = xZy(i,:);
j=j+1;%求得他的面积
end
end
if size(temp,1) == 2%求出size特殊情况
liadd(li+1) = sqrt(sum((temp(1,:)-temp(2,:)).^2));
%算出结果来
li = li+1;%得到数据
temp = [];%赋空
elseif size(temp,1) == 0%判断分支
li= -1;
else
disp("li警告");
size(temp,1)
disp(temp);
end
end
end
end
附录四 问题四的代码
main4.m 代码如下
% 清空工作区
clc, clear;
% 设置参数
% 设置水深、间隔
wx = 5*1852;%计算矩阵的宽
wy = 4*1852;%计算矩阵的长
% 计算结果的过程
jiaodu = 90;%jiaodu
kaijiao = [1.5, jiaodu, 120]*pi/180;%初始化参数
% 绘曲面图
figure(1);
x=(0:0.02:4.0)*1852;%计算长度
y=(0:0.02:5.0)*1852;%计算宽度
[xx,yy]=meshgrid(x,y);%计算大小
seah = -1*xlsread("附件.xlsx",'C3:GU253');%读取文件
sdaaas2 = 6;
sc =
surf(xx(1:sdaaas2:size(xx,1),1:sdaaas2:size(xx,2)),yy(1:sdaaas2:
size(xx,1),1:sdaaas2:size(xx,2)),seah(1:sdaaas2:size(xx,1),1:sdaaas2
:size(xx,2)));
ax = gca;
hold on
xlabel('东西方向');%x标
ylabel('南北方向');%y标
zlabel('水深')%z坐标
%% 拟合曲面
[fitr, gof] = createFit3(xx, yy, seah);%拟合的函数
fitr(1,2);%画出面
pars = coeffvalues(fitr);%计算大小
syms x y z%定义符号常量
eq1 = [1,x,y,x^2,x*y,y^2,x^3,x^2*y,x*y,y^3]*pars'-z;%造成函数，具体
nx = diff(eq1,x);%写入x坐标
ny = diff(eq1,y);%写入y坐标
nz = diff(eq1,z);%写入z坐标
[ix,iy] = find(seah==min(min(seah)));%寻找大小
nxfun111 = matlabFunction(nx);%利用模型来求解
nyfun222 = matlabFunction(ny);%同理用模型来求解
nxyASDmax = [nxfun111(xx(ix(1),iy(1)),yy(ix(1),iy(1))),
nyfun222(xx(ix(1),iy(1)) ,yy(ix(1),iy(1))),-1];
asin(nxyASDmax(3)/sqrt(sum(nxyASDmax.^2)))
%最后写出了算出他的值
% 计算等高线
zhi = seah(ix(1),iy(1));%计算他的值
itQQQemp = 1;%为itemp赋值1
zhiall(itQQQemp) = zhi;%将值写入矩阵中
itQQQemp = 2;%为itemp赋值2
while itQQQemp>0%循环
zhiall(itQQQemp) = zhi+39.8*0.896;%写出他的大小
if zhiall(itQQQemp) < 0%判断他是否小于零
if itQQQemp==2%若果是就计算
zhiall(itQQQemp) =
zhiall(itQQQemp-1)+13.3179;%赋值计算加上前面
else
zhiall(itQQQemp) =
zhiall(itQQQemp-1)-(zhiall(itQQQemp-2)-
zhiall(itQQQemp-1))*0.93;
end
itQQQemp = itQQQemp +1;
end
if zhiall(itQQQemp-1) > -20
itQQQemp = -1;
end
end
%初始化circle为0
circle = zeros(size(seah));
%利用三层循环来计算他的坐标与位置
for k =1:size(zhiall,2)-1
for i=1:size(seah,1)
for j=1:size(seah,2)
if abs( -seah(i,j)+0.37*zhiall(k)+ 0.63*zhiall(k+1)
)<0.3
circle(i,j) =circle(i,j)+1;
else
circle(i,j) =circle(i,j)+0;
end
end
end
end
%初始化circle为全面的几为0
circAHFle2 = zeros(size(circle));
itmp = 0;%为itmp赋值为0
for j=1:size(circle,2)%遍历数组
if circle(1,j)==1%写出来他的值
circAHFle2(1,j)=1;
end
if circle(i,end)==1%判断他的大小是否为1
circAHFle2(i,end)=1;%是1则，赋值
end
cri= diff(circle(:,j)');
for i=1:size(cri,2)
if cri(i)==1
tempq1 = i;
itQQQemp=1;
elseif cri(i)==-1
if itQQQemp==1
indexi = floor(0.37*i+0.63*tempq1);
circAHFle2(indexi,j)=1;
itQQQemp=0;
end
end
end
end
ASDE =0;
djx = 0.02*sqrt(2)*1852;
figure('Name','测量路线图')
for i=1:size(circAHFle2,1)
for j=1:size(circAHFle2,2)
if circAHFle2(i,j)==1%判断他们是否为1
ASDE =ASDE + djx*2.1;%计算差值，他们的和
scatter(xx(i,j),yy(i,j),10,'MarkerEdgeColor','k',
'MarkerFaceColor','k');
hold on
end
end
end
xlim([0,wx]);%x轴的范围
ylim([0,wy]);%y轴的范围
tempq1 = sum(circle,1);%zuihou 写出他的值
loss0= sum(tempq1<=20)*100/size(circle,2);
up20 = sum(30<=tempq1)*djx;
%% 显示
disp('总里程/m')
disp(ASDE/1000);
disp('漏测比例/%')
disp(loss0);
disp('重叠率超过20%的总里程/m')
disp(up20/1000);
createFit3.m 代码如下
function [fitresult, gof] = createFit3(xx, yy, data)
%CREATEFIT(XX,YY,DATA)
%“3次拟合”的数据如下所示：
%fitresult：代表拟合结果的拟合对象。
%gof：包含拟合质量信息的结构体。
%% Fit: '3次拟合'
[xData, yData, zData]= prepareSurfaceData( xx, yy, data );
% 建立例如迭代次数、收敛容限
ft= fittype( 'poly33' );
%将模型与数据进行拟合
[fitresult, gof]= fit( [xData, yData], zData, ft );
% 创建一个图表用于绘制图形。
figure( 'Name', '3次拟合' );
% 绘制拟合曲线与数据的图形。
subplot( 2, 1, 1 );
h = plot( fitresult, [xData, yData], zData );
legend( h, '3次拟合', 'data vs. xx, yy', 'Location', 'NorthEast',
'Interpreter', 'none' );
% Label axes
xlabel( 'x', 'Interpreter', 'none' );
ylabel( 'y', 'Interpreter', 'none' );
zlabel( '水深', 'Interpreter', 'none' );
grid on
view( 5.0, 16.6 );
% 绘制残差图。
subplot( 2, 1, 2 );
%自定义的函数类型
h = plot( fitresult, [xData, yData], zData, 'Style', 'Residual' );
legend( h, '3次拟合 - residuals', 'Location', 'NorthEast',
'Interpreter', 'none' );
%标注坐标轴。
xlabel( 'xx', 'Interpreter', 'none' );
ylabel( 'yy', 'Interpreter', 'none' );
zlabel( 'data', 'Interpreter', 'none' );
zlim([-0.1,0.1]);
grid on
view( 5.0, 16.6 );
扫一扫
684
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
