声波在水中的良好传播特性使其成为获取和传递水下信息最为有效的手段^[1-3]。声速是影响水下多波束系统作业精度的重要外部影响因素, 通过影响声线跟踪的精度, 最终影响到测深精度。多波束测深通过换能器实时接收其发射出的各波束经海底反射和散射后返回的到达角和旅行时^[3], 利用声速剖面数据, 由公式计算得到不同波束点对应的水深值。声速剖面是声速的垂直结构分布, 海水的介质特性导致声波的传播轨迹发生弯曲, 要获取波束脚印的确切位置, 需要沿着波束的传播路径追踪声线, 计算波束脚印的水平位移和深度, 即声线跟踪^[3]。因此, 声速剖面的正确与否直接影响多波束测深结果的精度和可靠性^[4]。

海洋中的声速剖面通常由以下两种方法获得: 直接测量法和间接测量法^[4]。直接测量法通过测量声波在海水中已知的固定距离内往返所需的时间或相位来计算声速值。间接测量法利用温盐深测量系统(conductivity, temperature, depth, CTD) 测量海水的温度、盐度及深度等物理量, 再通过声速经验公式计算作业水域内各采样深度的声速值, 进而构成声速剖面。然而在深水大洋进行实际作业时, 受到水深和洋流的影响, 停船投放CTD采集声速剖面数据时间成本较高。采用走航式的抛弃式声速剖面计虽然可以快速有效的获得声速剖面, 但其作业成本随之大大增加。相较于以上方法, 声速剖面重构只需少量的实测声速数据即可实现全海深声速剖面的重构, 提高效率的同时还有效的降低了成本。

由于海洋环境的复杂性, 声速剖面难以用一个简单的函数表示。Davis等^[5]证明EOF是在最小均方差意义下反映声速剖面最有效的基函数; 孙文川等^[6]论证了EOF重构的声速剖面具有较高的内符合精度, 能较好地描述实际声速剖面; 张志伟等^[7]研究表明取前6阶及以上的EOF表示声速剖面可满足多波束测深的精度要求; 张镇迈等^[8]采用部分深度的实测剖面, 并基于测区海域内已有的数据, 对全海深声速剖面实现了重构; 张孝首等^[9]利用实测CTD数据通过有理拟合延伸和EOF方法实现了深海5 000 m声速剖面场的重构; 李洪超等^[10]对不同深度位置的部分声速剖面的重构效果进行了探讨; 张维等^[11]在对残缺样本声速合理外延的基础上, 选取声速剖面变化较剧烈深度的3个声速值实现了剖面重构。以上研究表明, 基于EOF方法利用部分实测声速数据重构全海深剖面是可行有效的, 但未对参与重构的实测数据的采样深度给出明确依据。

本文将测区Argo浮标的温盐深数据通过声速经验公式计算得到的声速剖面作为测区历史声速剖面资料, 基于EOF空间函数的方差贡献率选取K(EOF阶次)个最大梯度对应深度处的声速值进行声速剖面重构, 在保证多波束测深精度的前提下实现了对全海深声速剖面的重构。

1 声速剖面EOF表示原理 1.1 声速剖面EOF分解

经验正交函数分析方法也称特征向量分析或主成分分析, 是分析矩阵数据的特性、提取主要数据特征量的一种方法^[12]。

在某海域采集N个声速剖面: C = [c₍₁₎, c₍₂₎, …, c_(N)], 内插成为M个垂直标准层, 得到声速矩阵C_M×N:

$\begin{array}{*{20}{c}} {{\boldsymbol{C}_{M \times N}} = \left[ {\begin{array}{*{20}{c}} {{c_{1(1)}}}&{{c_{2(1)}}}& \ldots &{{c_{N(1)}}} \\ {{c_{1(2)}}}&{{c_{2(2)}}}& \ldots &{{c_{N(2)}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{c_{1(M)}}}&{{c_{2(M)}}}& \ldots &{{c_{N(M)}}} \end{array}} \right]} , \end{array}$

(1)

式中, 每一列为一个声速剖面的标准层插值, 每一行为所有剖面在同一深度的声速。将N条声速剖面每一行进行平均, 得到平均声速剖面${\bar C_{M \times 1}}$:

$\begin{array}{*{20}{c}} {{{\bar C}_{M \times 1}} = \frac{1}{M}\mathop \sum \limits_{i = 1}^M {c_{\left( i \right)}}} , \end{array}$

(2)

扩展可得平均声速矩阵${\bar C_{M \times N}}$, 将声速矩阵与其相减, 得到每个剖面相对于平均剖面的扰动矩阵ΔC_M×N.

$\begin{array}{*{20}{c}} {\Delta {\boldsymbol{C}_{M \times N}} = {\boldsymbol{C}_{M \times N}} - {{\bar {\boldsymbol{C}}}_{M \times N}}}, \end{array}$

(3)

扰动矩阵的协方差矩阵为:

$\begin{array}{*{20}{c}} {{\boldsymbol{R}_{M \times M}} = \frac{{\Delta {\boldsymbol{C}_{M \times N}}\Delta \boldsymbol{C}_{M \times N}^T}}{N}} , \end{array}$

(4)

其中每个元素r_ij为:

$\begin{array}{*{20}{c}} {{r_{ij}} = \frac{1}{N}\mathop \sum \limits_{n = 1}^N \left[ {{c_n}\left( {{z_i}} \right) - \bar c\left( {{z_i}} \right)} \right]\left[ {{c_n}\left( {{z_j}} \right) - \bar c\left( {{z_j}} \right)} \right]} , \end{array}$

(5)

将R_M×M特征分解, 有:

$ \boldsymbol{R}_{M \times M} \boldsymbol{F}_{M \times N}=\boldsymbol{D}_{M \times M} \boldsymbol{F}_{M \times N}, $

(6)

式中, D_M×M为特征值矩阵, 令特征值λ_i按从大到小顺序排列:

$ \boldsymbol{D}_{M \times M}=\operatorname{diag}\left[\lambda_{1}, \lambda_{2}, \cdots, \lambda_{M}\right], $

(7)

F_M×N为特征值对应的特征向量矩阵, 也就是EOF空间函数, 可表示为:

$ \boldsymbol{F}_{M \times N}=\left[\boldsymbol{f}_{(1)}, \boldsymbol{f}_{(2)}, \cdots, \boldsymbol{f}_{(k)}\right], $

(8)

式中, f_(k)(k = 1, 2, …, N)为M×1的列向量, 即为确定的第k阶EOF。特征值λ_i对应的分量f_(i)即为第i阶EOF。特征值的大小代表了其对应EOF空间函数对线性空间的影响权重^[7], 即越大的特征值对应的特征向量中包含的声速重构信息越丰富。

1.2 声速剖面EOF重构

测区内的声速剖面完成EOF分解后, 利用前几阶EOF即可完成测区内任一声速剖面的重构:

$\begin{array}{*{20}{c}} {c\left( z \right) = \bar c\left( z \right) + \mathop \sum \limits_{i = 1}^K {\alpha _i}{f_i}\left( z \right), i = 1, 2, \cdots , M} , \end{array}$

(9)

式中, c(z)为重构的声速剖面; z为各层海水深度; $\bar c\left( z \right)$为平均声速; K为EOF阶次; α_i为重构系数; f_i(z)为EOF空间函数。移项得:

$\begin{array}{*{20}{c}} {\Delta {\boldsymbol{C}_{M \times N}} = {\boldsymbol{C}_{M \times N}} - {{\bar {\boldsymbol{C}}}_{M \times N}} = {\boldsymbol{F}_{M \times K}}{\boldsymbol{A}_{K \times N}}} , \end{array}$

(10)

则重构系数矩阵A_K×N为:

$\begin{array}{*{20}{c}} {{\boldsymbol{A}_{K \times N}}{\rm{ = }}\boldsymbol{F}_{M \times K}^T\Delta {\boldsymbol{C}_{M \times N}}} \end{array}{\rm{ = }}{\left[ {{a_1}, {a_2}, \cdots , {a_N}} \right]_{K \times N}}, $

(11)

式中, a_i = [a_i(1), a_i(2), …, a_i(K)]^T为每个声速剖面对应的重构系数。得到重构系数后, 即可按式(9)对测区内任一声速剖面进行EOF重构。

2 声速剖面EOF重构的数据采样深度选取方法

本文给出了用于声速剖面EOF重构的数据采样深度选取的方法步骤。算法流程图见图 1, 详细步骤和说明如下:

(1) 从中国Argo实时资料中心(http://www.argo.org.cn/)获取某海域温盐压数据, 预处理后利用声速经验公式计算得到采样位置处的声速剖面。目前国内外较为认可的声速经验公式有Chen-Millero, Wilson, Del Grosso及Leroy等, 各模型的温度、盐度、压力的适用范围有所差异。根据文献[13]的研究, Del Grosso, Chen-Millero, C.C.Leroy、Coppens四个经验声速模型在全球海域具有较高的精度及适用性, 文献[14]指出, 在350~1 000 m的深度上, Chen-Millero公式的误差最小, 因此本文选用Chen-Millero声速经验公式计算的声速剖面作为样本声速剖面, 利用Akima插值^[15]将声速剖面数据内插到垂直标准层。

(2) 按照1.1节中的方法进行EOF分解, 提取EOF函数, 计算各阶EOF的方差贡献率和累计方差贡献率。方差贡献率^[16]指的是每一阶EOF包含声速场信息的百分比。其计算公式为:

$\begin{array}{*{20}{c}} {{\sigma _i} = \frac{{{\lambda _i}}}{{\sum\limits_{i = 1}^M {{\lambda _i}} }} \times 100\% } , \end{array}$

(12)

式中, σ_i表示第i阶方差贡献率, λ_i表示协方差矩阵R经降序排列后的第i个特征值。

前K阶EOF对区域的累计贡献率可表示为:

$\begin{array}{*{20}{c}} {w(K) = {{\mathop \sum \limits_{i = 1}^K {\lambda _i}} \mathord{\left/ {\vphantom {{\mathop \sum \limits_{i = 1}^K {\lambda _i}} {\mathop \sum \limits_{j = 1}^M {\lambda _j}}}} \right. } {\mathop \sum \limits_{j = 1}^M {\lambda _j}}}} , \end{array}$

(13)

一般地, 如果w(K)的值大于等于95%, 则认为利用前K阶EOF能较好的表示该区域的主要信息^[5]。

(3) 计算测区内每一声速剖面每一深度层的声速梯度, 然后计算平均分层声速梯度:

${G_{i, n}} = \frac{{{c_{i + 1, n}} - {c_{i, n}}}}{{{z_{i + 1, n}} - {z_{i, n}}}}, $

(14)

$\overline {{G_i}} = \frac{1}{N}\sum\limits_{n = 1}^N {{G_{i, n}}} , $

(15)

式中, c_{i, n}和c_{i+1, n}分别表示第n条声速剖面第$i$层和第i+1层的声速值, z_i+1和z_{i+1, n}分别表示第n条声速剖面第i 层和第i+1层的深度值。按照平均分层声速梯度由大到小, 对声速层进行降序排列, 形成新的声速层矩阵。

(4) 从新构建的声速层矩阵中选取K个深度处的声速值(即选取K个声速变化最剧烈处的声速值), 结合K阶EOF计算得到重构系数矩阵A_K×N, 重构全海深声速剖面。

(5) 基于实测声速剖面, 预设水深值H和初始入射角, 根据分层常梯度声线跟踪法^[2], 以H作为约束迭代计算声波的传播时间T:

$ T={\displaystyle \sum _{i=1}^{Y-1}\frac{{\theta }_{i+1}-{\theta }_{i}}{p{g}_{i}^{2}\Delta {z}_{i}}\mathrm{ln}\left[\frac{{c}_{i+1}}{{c}_{i}}\right]}, $

(16)

式中, Y表示实测声速剖面的层数, θ_i和θ_i+1分别表示实测剖面第i层和第i+1层的入射角, g_i表示实测剖面第i层的梯度, c_i和c_i+1分别表示实测剖面第i层和第i+1层的声速, Δz_i表示第i层的深度差。p为Snell系数:

$ p=\frac{\sin{\theta }_{i}}{{C}_{i}}=\frac{\sin{\theta }_{i+1}}{{C}_{i+1}}.$

(17)

(6) 用重构的声速剖面和T反算波束水深值^[6]H′。声波在重构的声速剖面第i层的传播时间t_i为:

$ {t}_{i}=\frac{{{\theta }^{\prime }}_{i+1}-{{\theta }^{\prime }}_{i}}{p{{g}^{\prime }}_{i}^{2}\Delta {z}_{i}}\mathrm{ln}\left[\frac{{{c}^{\prime }}_{i+1}}{{{c}^{\prime }}_{i}}\right], $

(18)

式中, ${\theta '_i}$和${\theta '_{i + 1}}$分别表示重构声速剖面第i层和第i+1层的入射角, ${g'_i}$表示重构剖面第i层的梯度, ${c'_i}$和${c'_{i + 1}}$分别表示重构剖面第i层和第i+1层的声速。

$ T-{\displaystyle \sum _{i=1}^{J}{t}_{i}} =0,$

(19)

式中, J表示使等式成立的t_i的个数。所以:

$H' = \sum\limits_{i = 1}^J {\Delta {z_i}} .$

(20)

(7) 计算测深误差$\sigma $:

$\sigma = \frac{{\left| {H - H'} \right|}}{H} \times 100\% .$

(21)

定义有效波束比为满足深度限差的波束占总波束的比率^[17], 首先统计满足0.25%水深限差要求的波束数数, 然后判断多波束测深的有效波束比是否达到100%, 如果是, 则输出测深偏差, 实验结束, 反之, 则将选取阶数K加1, 重新完成重构系数矩阵的计算以及之后的步骤, 直至符合满足测深精度的要求。

3 实验分析

实验数据源于中国Argo资料实时中心发布的《全球Argo浮标剖面观测资料质量再控制数据集》^[18], 选取2006—2009年每年1月份的Argo浮标数据。将153.5°W—155W°, 27.5°N—29°N的海区作为区域1; 153.5°W—155.5W°, 22.5°N—24.5°N的海区作为区域2。去除不合格的剖面后, 区域1选取了18个剖面, 区域2选取了15个剖面。

典型深海声速剖面结构如图 2所示^[2], 主要划分为三部分: 混合层(由表面层和季节跃变层构成)、主跃层和深海等温层。其中, 声速在混合层和主跃层的变化情况复杂, 而在深海等温层中, 声速随深度的增加呈趋势稳定的正梯度增加, 接近线性变化。因此本文重点对混合层和主跃层这两部分声速变化情况更为复杂的水层的声速剖面进行研究, 统一截取海面至海深1 000 m的深度之间作为本次实验的剖面范围。

随机选取区域内1条剖面用作检核重构精度, 其余剖面作为已知剖面数据。实验区域内的Argo剖面浮标分布如图 3(a)和(b)所示, 图中“×”点代表已知的样本剖面, “▲”代表用于检核的剖面。

3.1 数据预处理

(1) Argo剖面浮标测量的是海水的温度、盐度和压力数据, 采用Chen-Millero声速经验公式^[19]计算得到相应深度处的声速值:

$\begin{gathered} c = {c_W}\left( {T, P} \right) + {c_C}\left( {T, P} \right) + A\left( {T, P} \right)S \\ + B\left( {T, P} \right){S^{3/2}} + D\left( P \right){S^2}, \\ \end{gathered} $

(22)

式中,

其中, T为温度, ℃; S为盐度; P为压力, bar。

(2) 将计算得到的声速数据利用Akima插值函数内插成1 m的等间距标准层。样本声速剖面结构如图 4。

3.2 确定参与重构的数据深度及声速值

分别对两个区域的样本声速剖面进行EOF分解, 得到测区EOF空间函数的方差贡献率和累计方差贡献率, 如表 1所示。可以看出, 在区域1中, 前5阶EOF的累计方差贡献率可达97.985%, 在区域2中, 前5阶EOF的累计方差贡献率达95.813%, 均已大于95%, 因此利用前5阶EOF能够较好的表示当前区域的主要信息。

图 5给出的是由两个实验海域内的样本声速剖面提取的1~5阶EOF, 可以明显看出, 第1阶EOF空间函数的变化较小, 整体呈趋势稳定, 这表明前几阶EOF中包含了实验海域内的声速场的主要信息, 而2~5阶EOF空间函数的变化明显, 变化集中在浅层区域, 尤其是在越靠后的阶次中, 这反映了声速场相对于平均SSP的细节变化^[7]。

计算每个样本声速剖面每层的梯度变化情况, 得到该区域的平均分层声速梯度及对应的深度位置, 按照平均分层声速梯度由大到小对声速层降序排列, 自上而下选取5个深度处的声速值, 选取结果如表 2所示。

3.3 声速剖面重构

选取表 2提取的声速值计算重构系数矩阵, 进而重构全海深声速剖面。两个实验海域内声速剖面的重构结果如图 6所示, 重构误差如图 7所示, 区域1的均方根误差为0.514 m/s, 区域2的均方根误差为0.609 m/s。可以看出: 重构声速剖面与实测声速剖面整体符合程度较好; 在两个区域的浅层部分, 尤其是区域1的100 m深度附近和区域2的350 m深度附近, 重构声速剖面与实测声速剖面有明显出入, 这可能是由于只选取了前几阶的EOF空间函数进行重构, 部分阶次的EOF空间函数并不能包含测区声速场的全部信息, 造成了部分信息的损失, 因而出现了一定的偏差。此外, 以500 m水深作为分界, 深层部分的声速剖面重构结果较为稳定, 整体的重构误差较小; 浅层部分的声速剖面重构结果相比较差, 造成该现象的原因可能是浅层部分的海水包含混合层和主跃层, 海洋要素在其中的变化情况复杂, 而在深层海水中已趋于稳定。

3.4 多波束测深检验

在多波束测深中, 利用EOF重构的声速剖面最终都要应用于声线改正。因此对于利用重构的声速剖面进行声线改正能否满足实际应用的要求, 还需要作进一步的检验。美国国家海洋与大气局(national oceanic and atmospheric administration, NOAA)对声速剖面引起的水深误差规定^[20]: “若使用的声速剖面与实际声速剖面对测深改正造成的互差超过0.25%水深则视为超限”。本文实验以NOAA的测深限差作为检验声速剖面重构精度的标准, 并结合有效波束比来判断利用EOF重构声速剖面的有效性。设置水深和波束的初始入射角, 基于实测声速剖面, 分层常梯度声线跟踪波束的往返时间, 利用该时间与重构的声速剖面再通过分层常梯度声线跟踪法反算水深, 并与真实水深值进行对比, 比较两者差值是否满足0.25%水深限差, 最后统计有效波束比。

表 3为使用两则区域重构的声速剖面进行多波束测深的结果, 可以看出, 使用重构的声速剖面进行声速改正得到的各波束点水深均符合限差要求, 有效波束比均达到了100%。

4 结论

本文对重构全海深声速剖面所需实测数据的采样深度进行了研究, 实验证明在对历史声速剖面资料进行分析后, 按照K阶EOF函数方差累计贡献率大于95%, 选取K个最大平均分层声速梯度处的声速值作为已知值, 重构的声速剖面能够满足NOAA对多波束测量声速剖面造成的误差要求, 为实际测量作业中声剖数据的采样深度提供了参考, 提高了工作效率。由于海洋声速易受各种环境复杂性的影响, 因此其他海区的选取阶次及深度会有所不同, 但本文介绍的方法和所得规律性结论具有一定的普适性。