CN103591942B - 基于声速最大偏移的声速剖面快速精简与自动优选方法 - Google Patents

Abstract

本发明公开了基于声速最大偏移的声速剖面快速精简与自动优选方法，给出了详细的综合技术流程，以解决因原始声速剖面数据量大而严重影响多波束探测与数据处理的工作效率问题。提出声速最大偏移量方法(MOV，MaximumOffsetofsoundVelocity)，用于自动、快速删除原始声速剖面的冗余点，并通过射线追踪和误差分析评估精简后的声速剖面对多波束测深的精度影响。经实际测试，处理后的声速剖面精简率达90%，经评估声速剖面精度能满足多波束脚印计算要求。该方法可大幅减少多波束勘测与数据处理时间，可有效提升多波束工作效率达3倍以上。该方法在海洋测绘、多波束勘测、海洋地理信息系统、计算机图形学和海底科学研究等方面具有重要的实际应用价值，可推广使用。

Description

基于声速最大偏移的声速剖面快速精简与自动优选方法

技术领域

本发明涉及海底地形地貌制图、海洋测绘、海洋地理信息系统、计算机图形学和海底科学等技术领域。

背景技术

声速剖面（Sound Velocity Profile,SVP）是实现多波束准确勘测必须的基本参数。声速剖面一般通过直接法（仪器测定）和间接法（声速经验模型）确定。目前，声速剖面测量设备的采样率普遍可达到20Hz，如果按1m/s的仪器下降速度算，100m水深所采集的声速剖面数据点的数量将达到2000组。在深水区进行多波束勘测，SVP数据量将非常庞大。高采样密度的声速剖面数据，导致射线追踪与波束脚印归算的运算时间大幅增加，从而降低了多波束勘测工作效率，有时甚至导致多波束系统不能正常运行。很多多波束测深系统对所使用SVP的数量提出了限制，如：最早引进我国的深水多波束系统SeaBeam 2112就限制声速剖面数据文件中声速数据不超过30行；又如：Kongsberg EM系列多波束系统对输入到处理单元用于实时射线追踪计算的SVP文件进行限制，要求声速剖面文件不能大于30kB，对于EM710、EM302和EM 122，最大采用数据点数量不超过1000组，对于更老型号的系统，则不超过570组。为了提升多波束勘测与数据处理的时效性，必须筛除原始声速剖面中的冗余点，同时需要评估和控制精简后的声速剖面所带来的误差。

目前对于声速剖面简化的方法多采取人工方法挑选声速剖面特征点，但人工方法挑选的SVP精度难以评估，也易遗漏SVP的特征点，同时效率也非常低。也有方法通过选定一定深度值大小的窗口，对声速数据进行滑动平均后，将所求的平均声速值作为窗口中心位置深度所对应的声速值，但此方法可能会遗漏原始声速剖面的特征点。

在申请专利“一种适用海底探测与假地形处理的实测SVP重构方法（申请号：201310152512）”中提到“使用D-P算法，对原始SVP剖面进行特征提取与拟合，调整曲线偏离因子D来拟合并简化原始SVP，保留原始SVP曲线的拐点”，该专利主要目的是对SVP曲线的特征提取而非删除曲线冗余点，同时该专利未给出SVP曲线简化的技术细节，最重要的一点是因为声速剖面的坐标维度不同，无法直接采用基于D-P算法的距离判断的方法来进行剖面精简，因为该算法所依赖的物理模型与数学模型已经发生改变，但该专利未给出如何解决坐标维度不同问题的具体方案。

发明内容

本发明公开了一种基于声速最大偏移的声速剖面快速精简与自动优选方法，以解决因原始测量的声速剖面数据量大而严重影响多波束测深的工作效率问题。基于声速维最大偏移量的计算来快速精简声速剖面，并通过射线追踪和均方差分析评估简化后的声速剖面精度，同时以阈值自动调节为纽带，实现了声速剖面的自动精简、评估与优选。

本发明是通过下述技术方案得以实现的：

基于声速最大偏移的声速剖面快速精简与自动优选方法，包括下列步骤：

1）形成原始声速剖面数据集

1.1)若有声速剖面，直接形成声速剖面数据集合SVP_in＝{in_svp_i}_i＝1,n，i为声速剖面序号，n为采集的声速剖面个数，i和n均自然数；

1.2)若无声速剖面，可采用声速剖面采集仪器，在海洋中获取原始声速剖面，并形成原始声速剖面数据集合SVP_in＝{in_svp_i}_i＝1,n；

1.3)每个声速剖面in_svp_i＝{P_j＝(d_j,v_j)}_j＝1,m，其中，P_j为声速剖面点，d_j和v_j为每个声速剖面点P_j所对应的深度值和声速值，m为声速剖面中的有效点数，j和m均为自然数；

1.4)输出一个声速剖面in_svp_i，转入步骤2）；

2)确定优选阈值区间

2.1）输入声速剖面in_svp_i；

2.2）遍历声速剖面in_svp_i，获取声速项的最小值v_s和最大值v_e；

T_step＝0.001×(v_e-v_s)，T_step为阈值自动计算步长；

初始化T_k＝0，T_k为当前声速精简阈值；

2.3）自动设置阈值：T_k＝T_k+T_step；

2.4)初始化当前声速阈值T_cur＝T_k，V_cur∈in_svp_i，V_cur是当前处理声速剖面段，V_cur＝{P_j＝(d_j,v_j)}_j＝a,b，a和b为自然数，是当前声速剖面段的首点和尾点，初始化V_cur＝in_svp_i＝{P_j＝(d_j,v_j)}_j＝1,m；

2.5)删除声速剖面冗余点：

2.5.1)提取当前剖面点数据集V_cur的首点P_a＝(d_a,v_a)和尾点P_b＝(d_b,v_b)；

2.5.2)遍历当前声速剖面数据集V_cur，依次取出每个声速剖面点P_j，使用公式(1)计算该点在声速维方向的偏移值D_j：

$D_j=\mathrm{abs}[\frac{(v_j-v_a)*(d_j-d_a)}{v_b-v_a}+d_a-v_j]---\left(1\right)$

将最大偏移值D_j存入D_max，并将相应的声速剖面点P_j存入P_k，P_k为临时声速剖面变量；

本发明的特征之一是基于声速维度最大偏移量的计算来保留声速剖面之特征点，以解决二维声速剖面中纵横坐标维度不同问题，因此本发明也称之为MOV方法(Maximum Offset ofsound Velocity)；

2.5.3)若D_max＞T_cur，添加P_k至V_tmp，V_tmp为过程数据集合，将剖面点数据集从P_k处分为两段，即V_cut1＝{P_j}_j＝a,k和V_cut2＝{P_j}_j＝k,b，将声速剖面段V_cut1和V_cut2分别赋值给V_cur，并分别返回步骤2.5.1)重新运算；

2.5.4)若D_max≤T_cur，添加P₁和P_m至V_tmp；

2.6)输出精简声速剖面：out_svp_i＝V_tmp＝{P_j＝(d_j,v_j)}_j＝1,mo；其中，mo为自然数，是精简后的声速剖面点数；out_svp_i是与in_svp_i对应在阈值为T_cur下减冗处理后形成的新声速剖面；

2.7)输出简化率： ${\mathrm{Per}}_i=(1-\frac{\mathrm{mo}}{m})\×100\%;$

2.8)获得简化率参数Per_k，添加至数据集合

2.9)当T_cur＜v_e-v_s时，返回步骤2.3）；

2.10)以T_k为横轴，以Per_k为纵轴，获取精简率曲线，并计算其二阶导数，获得二阶导数曲线

2.11）遍历二阶导数曲线获取其绝对值区间[f_min,f_max]，设置曲线截断值f_cut＝0.1×|f_max-f_min|；

2.12）根据二阶导数的曲线形态和震荡特征，仅保留二阶导数值小于f_cut的曲线段，并获取该段的优选阈值区间TT[T_min,T_max]；

2.13）输出优选阈值区间T＝[T_min,T_max]，转入步骤3）；

3）精简声速剖面

3.1）输入声速剖面in_svp_i和阈值区间T＝[T_min,T_max]；

3.2）设置T_step＝0.01×(T_max-T_min)，T_k＝T_min；

3.3）初始化T_cur＝T_k；初始化V_cur＝in_svp_i＝{P_j＝(d_j,v_j)}_j＝1,m

3.4)删除声速剖面冗余点：

3.4.1)提取V_cur的P_a＝(d_a,v_a)和P_b＝(d_b,v_b)；

3.4.2)遍历V_cur，依次取出P_j，使用公式（1）计算D_j，将最大偏移值D_j存入D_max，并将相应的声速剖面点P_j存入P_k；

3.4.3)若D_max＞T_cur，添加P_k至V_tmp；将剖面点数据集从P_k处分为两段，即V_cut1＝{P_j}_j＝a,k和V_cut2＝{P_j}_j＝k,b，将V_cut1和V_cut2分别赋值给V_cur，并分别返回步骤3.4.1)重新运算；

3.4.4)若D_max≤T_cur，添加P₁和P_m至V_tmp；

3.5)输出in_svp_i和out_svp_i，转入步骤4）；

4）评估声速剖面精度

4.1)输入原始声速剖面V_orig和精简声速剖面V_simp；

4.2)输入波束角集合B＝{θ_i}_i＝1,nb，nb为波束数，是自然数；

4.3)采用公式(2)，分别计算原始声速剖面V_orig和精简声速剖面V_simp的坐标值:(Orig_F_x_i ,Orig_F_d_i)和(Simp_F_x_i ,Simp_F_d_i)；

$\left\{\begin{array}{c}{F}_{\_x_i}=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{v}_j\×\sin \left({\mathrm{\α}}_i\right)\\ F\_d_i=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{v}_j\×\cos \left({\mathrm{\α}}_i\right)\end{array}\right.---\left(2\right)$

其中，α_i为波束角，其初始值为θ_i；v_j为声速值；

4.4)使用公式(3)，计算横向误差百分比ε_x_i和垂向误差百分比：

$\left\{\begin{array}{c}\mathrm{\ϵ}\_x_i=\frac{(\mathrm{Orig}\_F\_x_i-\mathrm{Simp}\_F\_x_i)}{\mathrm{Orig}\_F\_x_i}\×100\%\\ \mathrm{\ϵ}\_d_i=\frac{(\mathrm{Orig}\_F\_d_i-\mathrm{Simp}\_F\_d_i)}{\mathrm{Orig}\_F\_d_i}\×100\%\end{array}\right.---\left(3\right)$

4.5)对每个入射角{θ_i＝B_i}_i＝1,nb依次使用步骤4.3)至步骤4.5)，获得横向偏移误差数据集合{ε_x_i}_i＝1,nb和垂向误差数据集合{ε_d_i}_i＝1,nb;

4.6)使用公式(4)，计算横向误差百分比平均值μ_x及横向均方差百分比σ_x：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=\frac{{\mathrm{\Σ}}_{i=1}^{\mathrm{nb}}\mathrm{\ϵ}\_x_i}{\mathrm{nb}}\\ {\mathrm{\σ}}_x=\sqrt{\frac{1}{\mathrm{nb}}\underset{i=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}{(\mathrm{\ϵ}\_x_i-{\mathrm{\μ}}_x)}^2}\end{array}\right.---\left(4\right)$

4.7)使用公式(5)，计算垂向误差平均百分比μ_d及垂向均方差百分比σ_d：

$\left\{\begin{array}{c}{\mathrm{\μ}}_d=\frac{{\mathrm{\Σ}}_{i=1}^{\mathrm{nb}}\mathrm{\ϵ}\_d_i}{\mathrm{nb}}\\ {\mathrm{\σ}}_d=\sqrt{\frac{1}{\mathrm{nb}}\underset{i=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}{(\mathrm{\ϵ}\_d_i-{\mathrm{\μ}}_d)}^2}\end{array}\right.---\left(5\right)$

4.8)精度评定

若σ_d>0.1%，则T_k＝T_k-T_step，返回步骤3.4）；

若σ_d＜0.1%，则T_k＝T_k+T_step，返回步骤3.4）；

若σ_d=0.1%，输出V_simp，转入步骤5）；

5）依次处理声速剖面

5.1)将计算的最优声速剖面V_simp存入输出声速剖面数据集合SVP_out＝{out_svp_i}_i＝1,n，out_svp_i＝V_simp;

5.2)从原始声速剖面数据集合SVP_in＝{in_svp_i}_i＝1,n中顺序调入一个声速剖面，并转入步骤2），直至所有声速剖面处理完成，转入步骤6）；

6）使用声速剖面

将精简处理后的声速剖面数据集合SVP_out导入多波束勘测系统和数据处理系统，进行多波束测深勘测及数据处理。

本发明的有益效果

本发明提出并实现了基于声速最大偏移的声速剖面快速精简与自动优选方法。该方法的优势在于整体性、快速性和智能性，处理后的声速剖面的具有保真与高效特征。该方法能最大限度地保留原始声速剖面的最大特征点，能快速地从多种声速剖面组合中评估最佳的声速剖面，具有人工智能之特征。通过自动优选阈值，可删除原始声速剖面90%的冗余数据点，并通过射线追踪和均方差精度评估，证明精简后的声速剖面符合多波束探测精度要求。因此，该方法最大的优势是在保证声速剖面精度的前提下，删除90%的冗余点，从而达到大幅减少计算时间、有效提升工作效率之目的，在多波束勘测与数据处理中具有重要的工程应用价值。

附图说明

图1本发明实施例中的技术流程图；

图2是图1中的确定优选阈值区间技术流程图；

图3是图1中的精简声速剖面技术流程图；

图4是图1中的评估声速剖面精度流程图；

图5本发明实施例中不同阈值下声速剖面精简情况；

图6本发明实施例中射线追踪误差分析法示意图；

图7本发明实施例中使用的原始声速剖面数据集；

图8本发明实施例中阈值和声速剖面精简率以及均值二阶导数的关系；

图9本发明实施例中阈值vs水深标准误差百分比示意图；

图10本发明实施例中阈值与精简率和运算时间关系图。

具体实施方式

下面结合附图和实施例对本发明做进一步的说明。

实施例1

基于声速最大偏移的声速剖面快速精简与自动优选方法，包括下列步骤：

本发明完整的技术流程图见图1；

步骤一:形成原始声速剖面数据集

1.1)若有声速剖面，直接形成声速剖面数据集合SVP_in＝{in_svp_i}_i＝1,n，i为声速剖面序号，n为采集的声速剖面个数，i和n均自然数；

1.2)若无声速剖面，可采用声速剖面采集仪器，在海洋中获取原始声速剖面，并形成原始声速剖面数据集合SVP_in＝{in_svp_i}_i＝1,n；

1.3)每个声速剖面in_svp_i＝{P_j＝(d_j,v_j)}_j＝1,m，其中，P_j为声速剖面点，d_j和v_j为每个声速剖面点P_j所对应的深度值和声速值，m为声速剖面中的有效点数，j和m均为自然数；

1.4)输出一个声速剖面in_svp_i，转入步骤2）；

步骤二:确定优选阈值区间（详细技术流程如图2所示）

2.1）输入声速剖面in_svp_i；

2.2）遍历声速剖面in_svp_i，获取声速项的最小值v_s和最大值v_e；

T_step＝0.001×(v_e-v_s)，T_step为阈值自动计算步长；

初始化T_k＝0，T_k为当前声速精简阈值；

2.3）自动设置阈值：T_k＝T_k+T_step；

2.4)初始化当前声速阈值T_cur＝T_k，V_cur∈in_svp_i，V_cur是当前处理声速剖面段，V_cur＝{P_j＝(d_j,v_j)}_j＝a,b，a和b为自然数，是当前声速剖面段的首点和尾点，初始化V_cur＝in_svp_i＝{P_j＝(d_j,v_j)}_j＝1,m；

2.5)删除声速剖面冗余点：

2.5.1)提取当前剖面点数据集V_cur的首点P_a＝(d_a,v_a)和尾点P_b＝(d_b,v_b)；

2.5.2)遍历当前声速剖面数据集V_cur，依次取出每个声速剖面点P_j，使用公式(1)计算该点在声速维方向的偏移值D_j：

$D_j=\mathrm{abs}[\frac{(v_j-v_a)*(d_j-d_a)}{v_b-v_a}+d_a-v_j]---\left(1\right)$

将最大偏移值D_j存入D_max，并将相应的声速剖面点P_j存入P_k，P_k为临时声速剖面变量；

2.5.3)若D_max＞T_cur，添加P_k至V_tmp，V_tmp为过程数据集合，将剖面点数据集从P_k处分为两段，即V_cut1＝{P_j}_j＝a,k和V_cut2＝{P_j}_j＝k,b，将声速剖面段V_cut1和V_cut2分别赋值给V_cur，并分别返回步骤2.5.1)重新运算；

2.5.4)若D_max≤T_cur，添加P₁和P_m至V_tmp；

2.6)输出精简声速剖面：out_svp_i＝V_tmp＝{P_j＝(d_j,v_j)}_j＝1,mo；其中，mo为自然数，是精简后的声速剖面点数；out_svp_i是与in_svp_i对应在阈值为T_cur下减冗处理后形成的新声速剖面，如图5所示是在不同阈值大小下相应的减冗处理后形成的新声速剖面；

2.7)输出简化率： $\mathrm{Pe}{r}_i=(1-\frac{\mathrm{mo}}{m})\×100\%;$

2.8)获得简化率参数Per_k，添加至数据集合

2.9)当T_cur＜v_e-v_s时，返回步骤2.3）；

2.10)以T_k为横轴，以Per_k为纵轴，获取精简率曲线，并计算其二阶导数，获得二阶导数曲线

2.11）遍历二阶导数曲线获取其绝对值区间[f_min,f_max]，设置曲线截断值f_cut＝0.1×|f_max-f_min|；

2.12）根据二阶导数的曲线形态和震荡特征，仅保留二阶导数值小于f_cut的曲线段，并获取该段的优选阈值区间T＝[T_min,T_max]；

2.13）输出优选阈值区间T＝[T_min,T_max]，转入步骤3）；

步骤三:精简声速剖面（详细技术流程如图3所示）

3.1）输入声速剖面in_svp_i和阈值区间T＝[T_min,T_max]；

3.2）设置T_step＝0.01×(T_max-T_min)，T_k＝T_min；

3.3）初始化T_cur＝T_k；

初始化V_cur＝in_svp_i＝{P_j＝(d_j,v_j)}_j＝1,m

3.4)删除声速剖面冗余点：

3.4.1)提取V_cur的P_a＝(d_a,v_a)和P_b＝(d_b,v_b)；

3.4.2)遍历V_cur，依次取出P_j，使用公式（1）计算D_j，将最大偏移值D_j存入D_max，并将相应的声速剖面点P_j存入P_k；

3.4.3)若D_max＞T_cur，添加P_k至V_tmp；将剖面点数据集从P_k处分为两段，即V_cut1＝{P_j}_j＝a,k和V_cut2＝{P_j}_j＝k,b，将V_cut1和V_cut2分别赋值给V_cur，并分别返回步骤3.4.1)重新运算；

3.4.4)若D_max≤T_cur，添加P₁和P_m至V_tmp；

3.5)输出in_svp_i和out_svp_i，转入步骤4）；

步骤四:评估声速剖面精度（详细技术流程如图4所示）

4.1)输入原始声速剖面V_orig和精简声速剖面V_simp；

4.2)输入波束角集合B＝{θ_i}_i＝1,nb，nb为波束数，是自然数；

4.3)采用公式(2)，分别计算原始声速剖面V_orig和精简声速剖面V_simp的坐标值:(Orig_F_x_i ,Orig_F_d_i)和(Simp_F_x_i ,Simp_F_d_i)；

$\left\{\begin{array}{c}{F}_{\_x_i}=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{v}_j\×\sin \left({\mathrm{\α}}_i\right)\\ F\_d_i=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{v}_j\×\cos \left({\mathrm{\α}}_i\right)\end{array}\right.---\left(2\right)$

其中，α_i为波束角，其初始值为θ_i；v_j为声速值；

4.4)使用公式(3)，如图6所示，计算横向误差百分比ε_x_i和垂向误差百分比：

$\left\{\begin{array}{c}\mathrm{\ϵ}\_x_i=\frac{(\mathrm{Orig}\_F\_x_i-\mathrm{Simp}\_F\_x_i)}{\mathrm{Orig}\_F\_x_i}\×100\%\\ \mathrm{\ϵ}\_d_i=\frac{(\mathrm{Orig}\_F\_d_i-\mathrm{Simp}\_F\_d_i)}{\mathrm{Orig}\_F\_d_i}\×100\%\end{array}\right.---\left(3\right)$

4.5)对每个入射角{θ_i＝B_i}_i＝1,nb依次使用步骤4.3)至步骤4.5)，获得横向偏移误差数据集合{ε_x_i}_i＝1,nb和垂向误差数据集合{ε_d_i}_i＝1,nb;

4.6)使用公式(4)，计算横向误差百分比平均值μ_x及横向均方差百分比σ_x：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=\frac{{\mathrm{\Σ}}_{i=1}^{\mathrm{nb}}\mathrm{\ϵ}\_x_i}{\mathrm{nb}}\\ {\mathrm{\σ}}_x=\sqrt{\frac{1}{\mathrm{nb}}\underset{i=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}{(\mathrm{\ϵ}\_x_i-{\mathrm{\μ}}_x)}^2}\end{array}\right.---\left(4\right)$

4.7)使用公式(5)，计算垂向误差平均百分比μ_d及垂向均方差百分比σ_d：

$\left\{\begin{array}{c}{\mathrm{\μ}}_d=\frac{{\mathrm{\Σ}}_{i=1}^{\mathrm{nb}}\mathrm{\ϵ}\_d_i}{\mathrm{nb}}\\ {\mathrm{\σ}}_d=\sqrt{\frac{1}{\mathrm{nb}}\underset{i=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}{(\mathrm{\ϵ}\_d_i-{\mathrm{\μ}}_d)}^2}\end{array}\right.---\left(5\right)$

4.8)精度评定

若σ_d>0.1%，则T_k＝T_k-T_step，返回步骤3.4）；

若σ_d＜0.1%，则T_k＝T_k+T_step，返回步骤3.4）；

若σ_d=0.1%，输出V_simp，转入步骤5）；

步骤五:依次处理声速剖面

5.1)将计算的最优声速剖面V_simp存入输出声速剖面数据集合SVP_out＝{out_svp_i}_i＝1,n，out_svp_i＝V_simp;

5.2)从原始声速剖面数据集合SVP_in＝{in_svp_i}_i＝1,n中顺序调入一个声速剖面，并转入步骤2），直至所有声速剖面处理完成，转入步骤6）；

步骤六:使用声速剖面

将精简处理后的声速剖面数据集合SVP_out导入多波束勘测系统和数据处理系统，进行多波束测深勘测及数据处理。

实施例2

为了评估“基于声速最大偏移的声速剖面快速精简与自动优选方法”精简声速剖面对于多波束勘测与数据处理的影响，使用实施例1中之技术流程，采用实测的声速剖面进行精简与评估处理，并使用实测的多波束水深数据进行处理，评估不同阈值下的声速剖面对于多波束处理效率的影响，具体流程：

（1）形成原始声速剖面数据集：在实例中使用的实测的原始声速剖面点数据集进行方法检验。采集仪器为XR-420型温盐深仪，实测声速剖面共11个；对获取的11个声速剖面数据进行分析与分类，可分为3种类型（图7）：（1）正梯度型：声速随深度的增加而增大；（2）负梯度型：声速随深度的增加而减少；（3）小梯度型：声速随深度的增加而几乎不变；对所有数据形成方法可处理的声速剖面数据集。

（2）确定优选阈值区间：按如图2所示步骤，将实测的11个声速剖面按不同阈值进行处理，得到阈值和精简率以及二阶导数的关系图8。图8中右下方为对11个声速剖面精简率均值（图8黑色实线）的梯度再求导数，即二阶导数。到当阈值为0.04m/s，二阶导数的值突然减小到0.1左右；当阈值为0.18m/s时，此后的二阶导数的值逐渐趋近于零。由此可见，针对一个特定的SVP进行精简，其阈值存在一个区间[T₁,T₂]，T₁和T₂就是SVP精简率梯度发生突变的位置。分析这条曲线，我们就可以得出优选阈值区间。

（3）精简声速剖面：按如图3所示步骤，根据上一步所得的优选阈值区间，就可以得到11个声速剖面相应的精简声速剖面。

（4）评估声速剖面精度：图9为对所有实测声速剖面在选取不同阈值精简后，产生的误差百分比均值统计的结果。当阈值T∈[0,1]，随着阈值的增大，阈值和水深数据的标准差百分比呈线性关系，其斜率为0.2；当阈值T∈[1,7]，随着阈值的增大，阈值和水深数据的标准差百分比呈非线性关系。我们以水深标准误差百分比等于0.1%为条件进行精简，获得一系列声速剖面。

（5）使用声速剖面：为了检验精简前后的声速剖面对数据处理效率的影响，选取若干用于对比的原始声速剖面和精简声速剖面，多波束数据的采集设备为ElacBottomChart1180/1050双频中浅水多波束系统，水深范围为40-50m，选取其中40个测线文件，文件大小共390Mb，测线长度为498Km，共583.16万个有效波束点；采用多波束处理软件Caris HIPS 7.1中的声速改正模块进行数据处理，并统计运算时间。如图10所示，当阈值为0，即使用原始声速剖面时，声线追踪所需的时间为58s；随着阈值的不断增大，剖面精简率的不断提高，相应声线追踪所需的时间不断减少，且呈现非线性下降的趋势；当精简率为90%时，声线追踪所需的时间为17s。使用本发明中阈值与精简率和运算时间关系图见图10。

由此可见,使用本方法所精简的声速剖面在保证数据精度的前提下，能大幅减少数据处理的计算时间，数据处理的工作效率提升了3.41倍，工作效率提升显著；工作效率的提升对于多波束探测与处理的工程应用意义重要。

Claims (1)

1.基于声速最大偏移的声速剖面快速精简与自动优选方法，其特征在于，包括下列步骤：

1)形成原始声速剖面数据集

1.1)若有声速剖面，直接形成声速剖面数据集合SVP_in={in_svp_i}_i=1,n，i为声速剖面序号，n为采集的声速剖面个数，i和n均自然数；

1.2)若无声速剖面，可采用声速剖面采集仪器，在海洋中获取原始声速剖面，并形成原始声速剖面数据集合SVP_in={in_svp_i}_i=1,n；

1.3)每个声速剖面in_svp_i={P_j=(d_j,v_j)}_i=1,n，其中，P_j为声速剖面点，d_j和v_j为每个声速剖面点P_j所对应的深度值和声速值，m为声速剖面中的有效点数，j和m均为自然数；

1.4)输出一个声速剖面in_svp_i，转入步骤2)；

2)确定优选阈值区间

2.1)输入声速剖面in_svp_i；

2.2)遍历声速剖面in_svp_i，获取声速项的最小值v_s和最大值v_e；

T_step=0.001×(v_e-v_s)，T_step为阈值自动计算步长；

初始化T_k=0，T_k为当前声速精简阈值；

2.3)自动设置阈值：T_k=T_k+T_step；

2.4)初始化当前声速阈值T_cur=T_k，V_cur∈in_svp_i，V_cur是当前处理声速剖面段，V_cur{P_j=(d_j,v_j)}_j=a,b，a和b为自然数，是当前声速剖面段的首点和尾点；初始化V_cur=in_svp_i={P_j=(d_j,v_j)}_j=1,m；

2.5)删除声速剖面冗余点：

2.5.1)提取当前剖面点数据集V_cur的首点P_a=(d_a,v_a)和尾点P_v=(d_b,v_b)；

2.5.2)遍历当前声速剖面数据集V_cur，依次取出每个声速剖面点P_j，使用公式(1)计算该点在声速维方向的偏移值D_j：

$D_j=\mathrm{abs}[\frac{(v_j-v_a)*(d_j-d_a)}{v_b-v_a}+d_a-v_j]---\left(1\right)$

将D_j的最大偏移值存入D_max，并将相应的声速剖面点P_j存入P_k，P_k为临时声速剖面变量；

2.5.3)若D_max>T_cur，添加P_k至V_cmp，V_tmp为过程数据集合，将剖面点数据集从Pk处分为两段，即V_cut1={P_j|_j=,,k和V_cut2={P_j}_j=k,b，将声速剖面段V_cut1和V_cut2分别赋值给V_cur，并分别返回步骤2.5.1)重新运算；

2.5.4)若D_max≤T_cur，添加P₁和P_m至V_tmp；

2.6)输出精简声速剖面：out_svp_i==_tmp={P_j=(d_j,v_j)}_j=1,mo；其中，mo为自然数，是精简后的声速剖面点数；out_svp_i是与in_svp_i对应，是使用阈值为T_cur下减冗处理out_svp_i后形成的新声速剖面；

2.7)输出简化率： ${\mathrm{Per}}_i=(1-\frac{\mathrm{mo}}{m})\×100\%;$

2.8)获得简化率参数Per_k，添加至数据集合

2.9)当T_cur<v_e-v_s时，返回步骤2.3)；

2.10)以T_k为横轴，以Per_k为纵轴，获取精简率曲线，并计算其二阶导数，获得二阶导数曲线

2.11)遍历二阶导数曲线获取其绝对值区间[f_min,f_max]，设置曲线截断值f_cut=0.1×|f_max-f_min|；

2.12)根据二阶导数的曲线形态和震荡特征，仅保留二阶导数值小于f_cut的曲线段，并获取该段的优选阈值区间T=[T_min,T_max]；

2.13)输出优选阈值区间T=[T_min,T_max]，转入步骤3)；

3)精简声速剖面

3.1)输入声速剖面in_svp_i和优选阈值区间T=[T_min,T_max]；

3.2)设置T_step=0.01×(T_max-T_min)，T_k=T_min；

3.3)初始化当前阈值变量T_cur=T_k；

初始化当前声速剖面V_cur=in_svp_i={P_j=(d_j,v_j)}_j=1,m

3.4)删除声速剖面冗余点：

3.4.1)提取V_cur的P_a=(d_a,v_a)和P_b=(d_b,v_b)；

3.4.2)遍历V_cur，依次取出P_j，使用公式(1)计算D_j，将D_j的最大偏移值存入D_max，并将相应的声速剖面点P_j存入P_k；

3.4.3)若D_max>T_cur，添加P_k至V_cmp；将剖面点数据集从P_k处分为两段，即V_cut1={P_j}_j=2,k和V_cut2={P_j}_j=k,b，将V_cut1和V_cut2分别赋值给V_cur，并分别返回步骤3.4.1)重新运算；

3.4.4)若D_max≤T_cur，添加P₁和P_m至V_tmp；

3.5)输出in_svp_i和out_svp_i，转入步骤4)；

4)评估声速剖面精度

4.1)输入原始声速剖面V_orig和精简声速剖面V_simp；

4.2)输入波束角集合B=Pθ_i}_i=1,nb，nb为波束数，是自然数；

4.3)采用公式(2)，分别计算原始声速剖面V_orig和精简声速剖面V_simp的坐标值:(Orig_F_x_i,Orig_F_d_i)和(Simp_F_x_i,Simp_F_d_i)；

$\left[\begin{array}{c}{F}_{\_x_i}=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{v}_j\&times;\sin \left({\mathrm{\&alpha;}}_i\right)\\ F_{\_}{d}_i=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{v}_j\&times;\cos \left({\mathrm{\&alpha;}}_i\right)\end{array}\right.---\left(2\right)$

其中，α_i为波束角，其初始值为θ_i；v_j为声速值；

4.4)使用公式(3)，计算横向误差百分比ε_x_i和垂向误差百分比：

$\left\{\begin{array}{c}{\mathrm{\&epsiv;}}_{\_}{x}_i=\frac{(\mathrm{Orin}\_F\_x_i-\mathrm{Simp}\_F\_x_i)}{\mathrm{Orig}\_F{x}_i}\&times;100\%\\ {\mathrm{\&epsiv;}}_{\_}{d}_i=\frac{(\mathrm{Orih}\_F\_d_i-\mathrm{Simp}\_F\_d_i)}{\mathrm{Orig}\_F\_d_i}\&times;100\%\end{array}\right.---\left(3\right)$

4.5)对每个入射角{θ_i=B_i}_i=1,nb依次使用步骤4.3)至步骤4.5)，获得横向偏移误差数据集合{ε_x_i}_i=1,nb和垂向误差数据集合{ε_d_i}_i=1,nb；

4.6)使用公式(4)，计算横向误差百分比平均值μ_x及横向均方差百分比σ_x：

$\left\{\begin{array}{c}{\mathrm{\&mu;}}_x=\frac{{\mathrm{\&Sigma;}}_{i=1}^{\mathrm{nb}}\mathrm{\&epsiv;}\_x_i}{\mathrm{nb}}\\ {\mathrm{\&sigma;}}_x=\sqrt{\frac{1}{\mathrm{nb}}\underset{i=1}{\overset{\mathrm{nb}}{\mathrm{\&Sigma;}}}(\mathrm{\&epsiv;}\_x_i-{\mathrm{\&mu;}}_x{)}^2}\end{array}\right.---\left(4\right)$

4.7)使用公式(5)，计算垂向误差平均百分比μ_d及垂向均方差百分比σ_d：

$\left\{\begin{array}{c}{\mathrm{\&mu;}}_d=\frac{{\mathrm{\&Sigma;}}_{i=1}^{\mathrm{nb}}\mathrm{\&epsiv;}\_d_i}{\mathrm{nb}}\\ {\mathrm{\&sigma;}}_d=\sqrt{\frac{1}{\mathrm{nb}}\underset{i=1}{\overset{\mathrm{nb}}{\mathrm{\&Sigma;}}}(\mathrm{\&epsiv;}\_d_i-{\mathrm{\&mu;}}_d{)}^2}\end{array}\right.---\left(5\right)$

4.8)精度评定

若σ_d>0.1％，则T_k=T_k-T_step，返回步骤3.4)；

若σ_d＜0.1％，则T_k=T_k+T_step，返回步骤3.4)；

若σ_d=0.1％，输出V_simp，转入步骤5)；

5)依次处理声速剖面

5.1)将计算的最优声速剖面V_simp存入输出到声速剖面数据集合SVP_out=={{t_svp_i}_i=1,n，out_svp_i=V_simp；

5.2)从原始声速剖面数据集合SVP_in={in_svp_i}_i=1,n中顺序调入一个声速剖面，并转入步骤2)，直至所有声速剖面处理完成，转入步骤6)；

6)使用声速剖面

将精简处理后的声速剖面数据集合SVP_out导入多波束勘测系统和数据处理系统，进行多波束测深勘测及数据处理。