利用流形学习揭示声场数据的固有物理自由度

2021-07-13张新耀宋文华王宁胡涛

哈尔滨工程大学学报 2021年7期

张新耀，宋文华，王宁，胡涛

(1.中国海洋大学信息科学与工程学院，山东青岛 266100；2.大连测控技术研究所，辽宁大连116013; 3.中国科学院水声环境特性重点实验室，北京 100190)

浅海水声环境复杂，存在多样的海洋动力学过程及未知的海底地形和分层结构，浅海声场呈现出显著的不确定性，这使得模基(model-based)方法，如匹配场(matched field processing)算法[1-3]，性能大大降低。为了抑制不确定性对声学方法的影响，国内外学者提出了许多改进算法，如聚焦方法[4-5]和贝叶斯方法[6-7]，其思想是将环境参数纳入到搜索空间中，即对环境参数和目标参数同时搜索。然而对于复杂的浅海水声环境，可能影响声场的参数众多，这使得上述2种方法的计算复杂度呈几何式增长，难以满足实际应用中的时效性要求。因此需要在满足精度要求的前提下尽可能缩小搜索空间，即减少搜索参数的数目。

从物理角度而言，声场对于不同参数的敏感程度不同，因此国内外学者开展了声场对不同参数灵敏度的研究[8-10]，回答了主导浅海声场变化的物理参数有哪些问题。本文从微分流形的角度出发，利用流形学习方法揭示不确定浅海声场数据的固有物理自由度，来回答控制声场的参数有几个的问题。对于海试数据，在判定其固有物理自由度之后，本文进一步结合实验中得到的水文数据初步分析了控制声场的真实物理参数。

1 基于流形学习的声场数据固有物理自由度判定方法

不同于传统的函数分析视角，本文从微分流形视角来看待不确定浅海中的声传播问题。数学上，微分流形刻画了一大类几何结构：其局部类似于欧几里德空间，但是在全局上，一般具有更为复杂的几何结构。对于流形的描述有2种方式，即外在描述和内蕴描述。外在描述是将流形嵌入到更高维的空间中，该高维空间通常称为外围空间或观测空间；相反，内蕴描述不依赖于外围空间，仅依赖于流形的内蕴坐标(参数)。

1.1 目标-环境参数流形

在不考虑噪声影响的条件下，图1给出了微分流形视角下声传播问题示意图。从声场观测角度而言，当利用N阵元阵列对声场进行观测时，可以得到一组N维声场矢量，记为：

图1 微分流形视角下的声传播问题

p(γi)=[p1(γi)p2(γi)…pN(γi)]T,i=1,2，…,m

(1)

式中：γi表示第i个声场矢量相应的物理参数矢量，m表示观测快拍总数。这些N维声场矢量可以看作是N维观测空间中的点，其相应的N维坐标可以记作(p1,p2,…,pN)，这些点构成了嵌入在观测空间中的声场流形。声场是由声源与接收器之间的格林函数唯一决定，而格林函数是目标-接收器位置和环境参数的函数，因此有：

(p1,p2,…,pN)=(G1(γi),G2(γi),…,GN(γi))

(2)

式中：Gj(j=1,2，…,N)表示第j个接收器与声源之间的格林函数；γi表示与第i个声场矢量相对应的物理参数矢量。式(2)表明目标-接收器位置和环境参数即为声场流形的内蕴参数。因此，本文将声场流形称为目标-环境参数流形。

从微分流形角度而言，格林函数是一种微分映射，它将目标-环境参数空间映射到刻画声场的观测空间。尽管目标-环境参数空间本身一般视为欧几里德空间。然而从声场观测角度，声场对不同目标-环境参数的灵敏度不同，换言之，由浅海声场格林函数这种微分映射诱导的目标-环境参数空间并非欧几里德空间，而是一种微分流形或一种“超曲面”结构。曲面上的一点对应一组目标-环境参数下的声场，在该点的局地邻域内，曲面的弯曲程度刻画了声场在该目标-环境参数条件下，对不同参数的灵敏度：弯曲越大的方向表示声场对该参数越敏感，反之，越平坦的方向表示声场对该参数越不敏感；该点的局地切空间的维数定义了微分流形的局地维度。曲面的联络系数和曲率张量刻画不同参数之间的参数耦合。

1.2 利用流形学习方法判定固有维度

在不确定浅海中，相应的目标-环境参数流形也是未知的，因此声场数据所对应的格林函数一般包含不确定成份。但是这些声场数据是从流形上采样获得的，如果采样的声场数据足够多，则可以利用流形学习[11]重构出目标-环境参数流形。流形学习最早应用于图像处理领域[12-15]，在水声学领域鲜有应用，本文作者曾将其应用于声场数据的非线性降维[16]。判定数据固有维度是流形学习方法中的关键一步。等距嵌入方法[12]利用残差随维度增长而降低的收敛趋势来判定数据的固有维度。由于残差的计算依赖于数据点对之间的测地距离矩阵，其计算复杂度较高。Brand在文献[17]中给出了另外一种判定数据固有维度的方法，其指出数据流形的固有维度可以通过数据点对之间距离的增长过程来估计。该方法具有直观的几何图像，且易于算法实现，本文将这种方法应用于声场数据来判定其固有维度。

考虑N阵元水听器阵列，其接收到的一组声场矢量可以表示为p1,p2,…,pm∈RN，其中：

pi=p(γi)=[p1(γi),p2(γi),…,pN(γi)]T,

i=1,2，…,m

(3)

式中：γi表示第i个声场矢量相应的物理参数矢量，其对应的维度一般是未知的，记为l，即γi∈Rl。一般情况下，l≪N。

假设声场矢量所处的目标-环境参数流形M的维度为d，由于控制声场的不同物理参数之间可能存在参数耦合，因此：

d≤l≪N

(4)

假设该流形M是光滑流形，这意味着在某个空间尺度上，流形M上的一个邻域到Rd的映射是充分线性的，Brand[15]将该空间尺度称为局地线性尺度。在高维观测空间(此处观测空间维度为N)中，考虑以某个数据点为球心，以r为半径的球，将其包围的数据点数记为n(r)。在局地线性尺度下，n(r)将以rd的规律增大，即n(r)∝rd；当空间尺度小于局地线性尺度时，由于噪声的存在，数据点在各个方向的分布几乎是等概率的，此时则有n(r)∝rN，将该空间尺度称为噪声尺度；当空间尺度大于局地线性尺度时，流形曲率的影响会变得显著，因为流形将不再垂直于球的表面，此时n(r)的增长速度要快于rd规律，将该空间尺度称为曲率尺度。

为了定量地刻画球面包含数据点数目n(r)随空间尺度r增长的过程，Brand引入函数：

(5)

根据前文的分析论述，在噪声尺度下，n(r)∝rL，从而c(r)≈1/N<1/d；在局地线性尺度下，n(r)∝rd，从而c(r)=1/d；在曲率尺度下，由于n(r)的增长速度要快于rd规律，从而c(r)<1/d。因此，可以通过函数c(r)的最大值来估计流形的局地固有维度。需要注意的是，当球的半径r增大到包含整个流形时，会产生边界效应，此时c(r)会有1个小幅增长，所以在实际处理时，可依次以每个数据点为球心进行一次球扩展过程，寻找c(r)的第1个峰值，然后近邻的数据点取平均作为c(r)峰值的估计值，从而判定数据的局地固有维度。

2 数值仿真

本节所用到的仿真数据均是利用声场计算程序KRAKEN[18]生成的。假设考虑的不确定浅海波导为三维轴对称的水平均匀波导，此外，仿真生成的声场数据不包含噪声。声源频率固定为500 Hz，相应的波长λ≈3 m。

2.1 目标参数变化情形

图2给出了该仿真算例采用的浅海波导几何示意图。该浅海波导由一层水体和液态半空间海底组成，水深为100 m。垂直阵距离坐标原点的水平距离为5.6 km，由101个阵元组成，均匀地分布在整个水体深度上，阵元间隔为1 m。假设声源位置在图中矩形阴影区域I内变化，其相应的深度变化范围为50～100 m，水平变化范围为0～600 m。

图2 目标参数变化情形浅海波导几何示意

为了刻画垂直阵接收声场随声源位置的变化，将矩形阴影区域I以水平间隔Δr=3 m≈1λ和垂直间隔Δz=1 m≈1/3λ进行均匀采样，对于每一组声源位置γi=(rsi,zsi),rsi∈[3,600],zsi∈[50,100]，计算垂直阵处的仿真声场，记为：

pi=[p(z1;γi)p(z2;γi)…p(z101;γi)]T

i=1,2，…,10 200

(6)

共计得到10 200组复声场矢量。对其进行能量归一化处理并取实部，记为：

(7)

将这10 200组预处理之后的声场数据作为算法的输入，来判定其固有维度。

图3给出了数据集中某个数据点的球扩展点数增长曲线。函数c(r)即为球扩展点数增长曲线的斜率。由于仿真中未引入噪声，所以图3中的曲线未现出噪声尺度的特征。当半径r较小时，球包围的数据点数较少，所以呈现出阶梯式增长；随着半径r增大，进入局地线性尺度，c(r)达到峰值；当半径r继续增大时，进入曲率尺度，此时由于流形曲率的影响，c(r)减小；当r逐渐增大到到包含整个数据集时，c(r)≈1，随后曲线会有一个陡升。可以看出在该仿真中，声场数据的球扩展点数增长曲线与1.2节中描述的基本一致。

根据1.2节中的方法，函数c(r)峰值的倒数即为声场目标-环境参数流形的局地固有维度。图4给出了不同数据点处c(r)的峰值，图中给出的是相邻50个数据点平滑的结果。从图4中可以看出，大部分数据点对应的c(r)的峰值处于0.45左右。对整个数据集求平均，可得c(r)峰值的平均值为0.43，因此该仿真情形下，声场数据的固有维度为d≈1/0.43≈2.3，取d=2，这与真实控制声场的物理参数矢量γi=(rsi,zsi)的维度是一致的。

图4 目标参数变化情形下不同数据点处c(r)的峰值

2.2 环境参数变化情形

本节将利用声速剖面的变化来模拟环境的变化。图5给出了该情形下浅海波导的示意图，其由一层水体和液态半空间海底组成，水深为35 m。垂直阵距离坐标原点的水平距离为5 km，由36个阵元组成，均匀地分布在整个水体深度上，阵元间隔为1 m。声源深度固定在35 m。为了定量地刻画声速剖面的变化，利用2007年夏季黄海海洋环境与声传播实验中实测的一组声速剖面，通过经验正交函数(empirical orthonormal function, EOF)构建多个声速剖面来计算仿真声场。

图5 环境参数变化情形下浅海波导示意

图6和图7分别给出了该海上实验测得的声速剖面的背景声速剖面和前2阶经验正交函数，根据公式来构建仿真中采用的多组声速剖面：

图6 2007年黄海实验背景声速剖面

图7 声速剖面前2阶经验正交函数

(8)

式中：C0(z)表示背景声速剖面；Ψj(z)为第j阶经验正交函数；βij为相应的系数且βi1∈[-10,10]，βi2∈[-5,5]。

通过控制βi1和βi2的取值来实现声速剖面的变化，将其对应的2个取值区间分别以步长Δβ1=0.2和Δβ2=0.1进行离散化。对于每一组参数γi=(βi1,βi2)，计算垂直阵处的仿真声场：

pi=[p(z1;γi)p(z2;γi)…p(z36;γi)]T

i=1,2,…,10 201

(9)

共计得到10 201组复声场矢量。对其进行能量归一化处理并取实部，记为：

(10)

将这10 201组预处理之后的声场数据作为算法的输入。

图8给出了该仿真情形下数据集中某个数据点的球扩展点数增长曲线。由于该仿真中同样未引入噪声，所以图8中的曲线未呈现出噪声尺度的特征。当半径r较小时，球包围的数据点数较少，所以呈现出阶梯式增长；随着半径r增大，进入局地线性尺度，球扩展点数增长曲线的斜率c(r)达到峰值；当半径r继续增大时，进入曲率尺度，此时由于流形曲率的影响，c(r)减小；当r逐渐增大到包含整个数据集时，c(r)≈1，随后曲线会有一个陡升。可以看出在该仿真中，声场数据的球扩展点数增长曲线同样与1.2节中描述的基本一致。

图8 声速剖面不确定情形下球扩展数据点数增长曲线

图9给出了该仿真情形下不同数据点处c(r)的峰值，图中给出的是相邻50个数据点平滑的结果，整个数据集上c(r)峰值的平均值为0.51，因此该仿真情形下，声场数据的固有维度为d≈1/0.51≈1.96，取d=2。这与真实控制声场的物理参数矢量γi=(βi1,βi2)的维度是一致的。

图9 声速剖面不确定情形下不同数据点处c(r)的峰值

将上述2次仿真的结果在表1中进行汇总。

表 1 数值仿真结果汇总

3 海上实验数据处理

3.1 2007年黄海海洋环境和声传播实验

2007年黄海海洋环境和声传播实验是由中国科学院声学研究所组织开展的。实验中对北纬35°东经121°附近海域的海洋环境进行了连续观测，同时还进行了8.7 km的声传播实验，以观测内波等海洋物理过程对声场的影响。本次实验采用定点观测，本节中采用的数据主要是16阵元垂直阵记录的声传播数据和1号温度链记录的温度数据。实验中未对海底底质和海深进行原位测量，但通过TD和CTD记录的数据可以估计出实验海域海深约为36 m。16个阵元大致分布在水深5～28 m，阵元间距为1.5 m。发射换能器中心频率为300 Hz，采用坐底式发射，换能器距海底约65 cm。发射信号以60 s为周期进行循环发射，包括：5 s的线性调频信号，带宽为260～340 Hz，继之以6.8 s的间歇；1 s的300 Hz单频信号，继之以6.9 s的间歇；5 s的300 Hz单频信号，继之以10 s的间歇；8.5 s的300 Hz调制的8阶M序列伪随机信号，继之以7.9 s的间歇；0.1 s的短脉冲信号，继之以9.8 s的间歇。共计发射223组信号，期间，部分时间监视信号波形失真，对此部分信号进行剔除，最终有效发射周期为212个。

由于此次声传播时间实验持续时间较短，为了充分利用发射周期内的不同信号，对接收数据作如下处理：对于每个周期的信号，截取前5 s的线性调频信号，经匹配滤波后，进行傅里叶变换，取300 Hz频率的能量归一化信号的实部作为第1组数据；截取11.8～12.8 s的300 Hz单频信号，进行傅里叶变换，对其能量归一化之后取实部作为第2组数据；对于19.7～24.7 s的300 Hz单频信号，从20.2～24.2 s，每1 s截取一组信号，进行傅里叶变换，对其能量归一化之后取实部作为第3～6组数据。经过这样的处理之后，共计得到6组300 Hz的单频声场数据集，每组包含212个16维的声场数据，总计1 272个数据样本。利用1.2节中的方法判定该声场数据集的固有维度。

图10给出了以数据集中某个样本点为中心的球扩展数据点数增长曲线。由于数据样本较少，在r较小时，曲线并未呈现出噪声尺度的特征；当logn(r)的取值在1～2时，曲线的斜率c(r)达到峰值，约为0.3；当r继续增大时，进入曲率尺度，此时曲线的斜率c(r)减小。

图10 2007年黄海声传播实验球扩展数据点数增长曲线

图11给出了以不同的数据点为中心得到的c(r)的峰值的变化，在整个数据集上，c(r)峰值的平均值为0.31，因此该声场数据集的固有维度为d≈1/0.31≈3.2，取d=3。

图11 2007年黄海声传播实验不同数据点处c(r)的峰值

在确定了声场数据集的固有维度d=3后，将结合实验中记录的声源位置数据及声速剖面数据，对潜在的控制声场的物理参数进行分析。

首先考虑声源位置的变化。在声传播实验过程中，发射船“金星二号”上的GPS记录了发射船与垂直阵之间水平距离的变化，如图12所示。从图中可以看出，在信号发射期间，发射船与垂直阵之间的距离变化最大约为50 m，占总的声传播距离(8.7 km)的0.57%。

图12 声传播实验期间发射船与垂直阵之间距离的变化

其次，在声传播实验期间，实验海域约有3 m的潮差，图13为垂直阵上TD记录的深度变化，在声传播实验中深度变化达到了4 m，超过了水体深度(约为36 m)的10%，这必然会对声场数据产生显著影响。

图13 声传播实验期间垂直阵上TD记录的深度变化

第三,考虑水体声速剖面的变化。图14给出了声传播实验期间利用垂直阵附近的1号温度链采集的温度数据结合声速经验公式得出的声速剖面变化。从图中可以看出，在声传播实验期间存在强烈的内波活动。通过对1号温度链记录的4 600多组数据进行主成分分析，其前2阶主成分的能量占比达到91.75%，相应的前2阶经验正交函数如图7所示。声传播实验期间，垂直阵附近声速剖面的前2阶主成分β1和β2的变化如图15所示，二者具有可比的变化幅度。

图14 声传播实验期间的声速剖面变化

图15 声传播实验期间声速剖面前2阶EOF系数β1和β2的变化

通过以上分析，可以看出控制此次声场实验数据的物理参数有4个：声源与垂直阵之间的距离、水深、垂直阵附近声速剖面的前2阶主成分β1和β2。然而，前文中判定的声场数据的固有维度d=3。这说明这4个物理参数不是互相独立的。Del Balzo[19]曾分析过水深失配对匹配场定位结果的影响，其指出：过度估计的水深会导致声源定位结果相比于真实距离较远，反之，低估的水深会导致声源定位结果较近。这从侧面说明，声源距离与水深之间存在着耦合关系。此外，考虑到此次声传播实验中声源距离的变化(50 m)相比于声源距离(8.7 km)仅为0.57%。综合以上几点分析，控制2007年黄海声传播实验中声场数据的潜在的3个物理自由度应为水深、声速剖面的前2阶主成分β1和β2。

3.2 2019年9月南海北部三维声场获取实验

2019年9月南海北部三维声场获取实验是由中科院声学研究所、中国海洋大学等多个单位联合开展的海上实验，旨在获取南海北部的三维声场信息，观测南海内波对三维声场的影响。实验中各站位的分布如图16所示。

图16 2019年南海北部三维声场获取实验站位分布

O站位为发射站位，水深84.6 m，布放有300 Hz低频换能器和温度链。B1、A2和A3站位均为接收站位，本节采用的实验数据为9月9日0∶20至9月13日23∶21 O站位发射B1站位接收的定点声传播实验数据。O站位与B1站位之间水平距离约为24.5 km，B1站位水深87.6 m，布放的垂直阵为不等间距的32阵元充油阵(实验期间有一个水听器发生故障)，分布在水深16～63.5 m。发射换能器发射信号以2 min为一个周期，包括：4 s线性调频信号，带宽为280～350 Hz，继之以76 s的间歇；2 s的310 Hz的单频信号，继之以38 s的间歇。

对数据作如下处理：对线性调频信号进行匹配滤波，然后进行傅里叶变换，为了与单频信号频率保持一致，取310 Hz频率的数据，能量归一化后取其实部作为第1组数据；对310 Hz单频信号作傅里叶变换，能量归一化后取其实部作为第2组数据。共计得到7 200组31维矢量，利用1.2节中的方法判定该数据集的固有维度。

图17给出了以数据集中某个样本点为中心的球扩展数据点数增长曲线。在r较小时，曲线并未呈现出噪声尺度的特征；当logn(r)的取值在1～1.5时，曲线的斜率c(r)达到峰值，约为0.2；当r继续增大时，进入曲率尺度，此时曲线的斜率c(r)减小。

图17 2019年南海实验球扩展数据点数增长曲线

图18给出了以不同的数据点为中心得到的c(r)的峰值的变化，在整个数据集上，c(r)峰值的平均值为0.19，因此该声场数据集的固有维度为d≈1/0.19≈5.2，取d=5。

图18 2019年南海实验不同数据点处c(r)的峰值

在确定了声场数据集的固有维度d=5后，类似于对2007年黄海实验数据的处理，将结合实验中记录的声源位置数据及声速剖面数据，对潜在的控制声场的物理参数进行分析。

首先考虑声源位置的变化。图19给出了定点声传播实验期间发射船与垂直阵之间的距离变化，距离最近点与最远点之间的差值为276.7 m，占总传播距离(24.5 km)的1.1%。

图19 2019年南海实验发射船与垂直阵之间的距离变化

其次考虑水深的变化。图20给出了B1站位垂直阵上端压力计记录的深度变化，定点声传播实验期间水深变化的最大值达到4.3 m，占总水深(87.6 m)的4.95%, 这会对声场产生显著影响。

图20 2019年南海实验垂直阵上TD记录的深度变化

图21 B1站位由温度链记录的声速剖面变化

图22 B1站位声速剖面的前4阶经验正交函数

通过以上分析，可以看出控制此次南海声场实验数据的物理参数有6个：声源与垂直阵之间的距离、水深、垂直阵附近声速剖面的前4阶主成分β1、β2、β3和β4。然而，前文中判定的此次声场实验数据的固有维度d=5。类似于对2007年黄海实验数据的分析，声源与垂直阵之间的距离与其他5个参数是耦合的，而且其变化范围相对较小，因此将控制2019年9月南海定点声传播实验中声场数据的潜在的5个物理自由度判定为水深、声速剖面的前4阶主成分β1、β2、β3和β4是合理的。将2次海试数据处理结果在表2中进行汇总。