普氏蹄蝠声呐脉冲序列与结构形变的模式匹配

刘彦茗，张舒心，杨璐慧

(1. 山东大学机械工程学院，山东济南250100；2. 山东大学物理学院，山东济南250100；3. 山东建筑大学机电工程学院，山东济南250101)

0 引 言

声呐系统是蝙蝠在漫长进化过程中取得成功的突出优势之一，目前利用声呐系统进行回声定位的蝙蝠约有1 200种[1]，且分布在地球的绝大部分地区[2]。而目前基于生物声呐系统研制的工程产品，如声呐、雷达，在信息编码效率、探测性能和结构尺寸等方面较生物声呐系统仍有巨大差异[3]。因此，蝙蝠声呐系统是改善现有工程传感系统的绝佳生物模型，可以为研发更优异的仿生声呐系统和工程声呐系统提供灵感[4]。

蹄蝠科和菊头蝠科蝙蝠能够在非常复杂的自然环境中自由穿梭并高效地完成目标识别、猎物捕捉等传感任务，它们都是由鼻孔发射超声声呐信号，利用耳朵接收回声信号，通过分析回声信号获取关于周围环境的传感信息。其声呐系统都具有独特的动态特性，主要表现为：(1) 蝙蝠能够根据周围环境和传感任务的不同，动态调整超声声呐信号的强度、方向、持续时间和发声频率等参数[5-6]，如在植被较密集的区域或猎物的捕捉阶段，蝙蝠发射的超声脉冲信号更多，发声频率更高[7]；(2) 蝙蝠的声呐系统结构，即鼻孔周围挡板形状的褶皱结构(即“鼻叶”)和耳廓，如图1所示，在超声信号发射和回声信号接收过程中受肌肉的控制可以发生不同形态的动态形变，如鼻叶可发生闭合-开启运动(即合开运动)和随机运动[8]；耳廓可发生刚性运动和柔性运动[9-10]。数值模拟和仿生研究表明鼻叶和耳廓运动能够显著地改变声场分布，并为蝙蝠提供额外有用的传感信息[9,11]，生物研究发现鼻叶运动和耳廓运动存在一定的耦合关系，即超声发射系统和回声接收系统的结构动态特性之间存在耦合关系[12]，而蝙蝠声呐系统的两种动态特性(即动态声呐信号和动态结构形变)之间是否存在耦合关系仍然是未知的。

本文以普氏蹄蝠为研究对象，利用高速摄像机和超声麦克风同步采集蝙蝠的鼻叶/耳廓运动数据及其对应的声呐脉冲序列数据；针对出现频率较高的两种鼻叶-耳廓耦合运动模式，即鼻叶合开-耳廓柔性运动、鼻叶随机-耳廓刚性运动，利用一种专门用于神经脉冲序列模式分类(Spike train SIMilarity Space, SSIMS)的方法[13]，将声呐脉冲序列中每个声呐脉冲看作神经元的一个脉冲，对两种耦合运动模式的声呐脉冲序列进行分类，以研究动态结构形变与动态声呐信号之间是否存在耦合关系。

1 实验方法

1.1 实验材料

本文以14只雄性普氏蹄蝠为实验对象，它们来自江西省永丰县梅坑村宗祠(N27°6', E115°37')和福建省沙县大佑山脚洞(N26°25', E117°39')，被分别安置在两个长6.9 m、宽1.3 m、高3.2 m的饲养室中。为给蝙蝠提供健康稳定的生活环境，饲养室常年保持相对恒温(22℃～25℃)、恒湿(60%～70%)，且蝙蝠每天摄入定量的大麦虫和添加维生素的矿泉水，以保证充足营养和微量元素。因蝙蝠属于昼伏夜出动物，为了保证白天实验时蝙蝠处于兴奋状态，颠倒饲养室的照明时间，即当地时间6:00-18:00关灯；18:00-6:00(次日)照明。普氏蹄蝠鼻叶/耳廓结构及其运动度量如图1所示。

图1 普氏蹄蝠鼻叶/耳廓结构及其运动度量Fig.1 Noseleaf/Pinna structure in Hipposideros pratti and their motion metrics

1.2 实验装置

为了研究蝙蝠发射的声呐脉冲序列与鼻叶/耳廓耦合运动模式之间的关系，本文使用图2所示的实验装置，同步采集蝙蝠发射的声呐脉冲序列、鼻叶运动和耳廓运动。该实验装置主要由3个高速摄像机(GigaView, Southern Vision Systems, Inc, 美国)、1个超声麦克风(Type 40DP，G.R.A.S Sound &Vibration, 丹麦)、NI数据采集系统(PXIe-6358)和辅助设备等组成。其中，高速摄像机的帧频均设置为400帧·s-1，并保证至少2个高速摄像机能清晰地捕获鼻叶或耳廓运动，以便后续对鼻叶和耳廓运动轨迹的三维重建；超声麦克风用于同步采集蝙蝠发射的声呐脉冲序列；普氏蹄蝠使用的超声声呐信号为恒频-调频(Constant Frequency-Frequency Modulation, CF- FM)信号，能量主要集中在恒频(Constant Frequency, CF)段，调频(Frequency Modulation, FM)段的能量很少，且第二谐波(CF约为 61 kHz，FM约49～61 kHz)的能量最强。根据采样定理，对声音数据的采样频率设置为250 kHz。

图2 采集蝙蝠声呐脉冲序列和鼻叶/耳廓运动的实验装置Fig.2 Experimental setup for collecting sonar pulses and noseleaf/pinna motions

1.3 实验步骤

采集蝙蝠声呐脉冲序列和鼻叶/耳廓运动的实验步骤如下：

(1) 将蝙蝠从饲养室带入实验室，并在蝙蝠左、右耳廓、鼻叶、头部分别做白色标记点，其中，左、右耳廓各8个，鼻叶9个，如图1(a)所示，用以三维重建鼻叶和耳廓的运动轨迹，头部标记点作为判断蝙蝠头部是否运动的参考点；

(2) 打开高速摄像机和麦克风等实验设备，检查设备、Matlab程序是否运行正常；

(3) 将蝙蝠放置在蝙蝠固定支架上，调整支架高度使蝙蝠鼻孔与超声麦克风处于同一高度；

(4) 在蝙蝠正前方 60 cm左右，用扬声器发出线性调频信号(持续时间0.3 s，频率范围8～18 kHz)或人工飞蛾等诱导蝙蝠发射超声信号并产生鼻叶和耳廓运动；

(5) 当目测蝙蝠鼻叶和耳廓均有较大形变时，启动Matlab程序同时触发高速摄像机和麦克风，分别记录并保存3 s的鼻叶/耳廓运动图像数据和超声数据；

(6) 实验结束后，清除蝙蝠身上的所有标记点，并喂入2～3只大麦虫和充足矿泉水，以补充能量和水分，送回饲养室。

2 实验数据处理

2.1 图像数据处理

为了精确地获取蝙蝠的鼻叶运动和耳廓运动，只对蝙蝠头部静止不动且鼻叶或耳廓运动的图像数据进行分析，即首先根据蝙蝠头部标记点是否运动，筛选出蝙蝠头部静止不动的图像数据；然后根据鼻叶运动和耳廓运动的差异，对筛选出的图像数据进行分类。其中，鼻叶运动主要有合开运动和随机运动两种，合开运动是指鼻叶各部分结构(上鼻叶、鼻孔两侧叶瓣、下鼻叶等)同时朝向鼻孔方向运动或背离鼻孔方向运动，如图 3(a)所示，即鼻叶宽度和高度同时减小或增大，如图4(a)所示。鼻叶随机运动是指鼻叶各部分结构运动方向不一致，区别于合开运动和鼻叶静止(最大鼻叶宽度和高度变化均小于0.02 mm[10])，其部分运动模式如图 3(b)所示，鼻叶宽度和高度也会出现同时增大或此变彼不变等现象，如图4(b)所示。耳廓运动主要有柔性运动和刚性运动，柔性运动是指耳廓结构发生非刚性形变，如图3(c)所示，即耳廓宽度发生明显的增大或减小，如图4(a)所示。刚性运动是指耳廓前后摆动，同时结构形状几乎不变，如图3(d)所示，即耳廓宽度几乎不变，如图4(b)所示，耳廓宽度归一化距离变化量小于0.1。

图3 鼻叶运动和耳廓运动分类Fig.3 Classification of noseleaf and pinna motions

图4 鼻叶运动和耳廓运动定量描述Fig.4 Quantitative description of noseleaf and pinna motions

研究发现[12]，鼻叶合开运动常伴随耳廓柔性运动，而鼻叶随机运动则大多伴随耳廓刚性运动，因此，本文选用鼻叶合开-耳廓柔性运动、鼻叶随机-耳廓刚性运动两种鼻叶/耳廓耦合运动模式为研究对象，分析蝙蝠声呐系统结构动态形变与其对应声呐脉冲序列是否存在耦合关系。

2.2 声呐数据处理

蝙蝠声呐信号的参数主要有声呐脉冲强度、时长、频率结构、发声频率等，由于本文旨在研究蝙蝠声呐脉冲序列与其动态结构形变之间的耦合关系，侧重声呐脉冲序列的时域结构，即只研究鼻叶/耳廓耦合运动期间蝙蝠发射声呐脉冲数目和激发时刻，忽略其他参数。因此，声呐数据处理的具体步骤如下：

(1) 根据高速摄像机的帧频(400帧·s-1)，将2.1小节筛选出的两种鼻叶/耳廓耦合运动模式的图片集，转换成对应声呐信号的时间段；

(2) 从采集的原始声音数据(时长为 3 s)中截取上述时间段的声呐数据，剔除不存在声呐脉冲的数据，并对存在脉冲的声呐数据进行滤波处理(基于汉明窗设计的FIR带通滤波器，512阶，中心频率为61 kHz，通带宽度4 kHz)，以去除环境噪声并保留第二谐波的恒频成分(61 kHz)，结果如图5(a)所示；

(3) 对滤波后的声呐数据进行希尔伯特变换后取绝对值，得到该声呐数据的包络线，结果如图5(b)所示；

(4) 以10倍于噪声有效值的阈值，截取该声呐数据中每个脉冲的起止时刻，结果如图 5(c)所示，以起止时刻的平均值作为每个脉冲的激发时刻。

图5 声呐数据处理示意图Fig.5 Schematic diagrams of sonar data processing

通过上述图像数据和声呐数据处理，建立两种鼻叶/耳廓耦合运动模式及其声呐脉冲序列的数据库，该数据库共有来自14只雄性普氏蹄蝠的66次实验的102组数据，其中，53组鼻叶合开-耳廓柔性运动，49组鼻叶随机-耳廓刚性运动，如图 6所示。不同鼻叶/耳廓耦合运动模式中声呐脉冲序列的差异主要表现为脉冲数目和脉冲激发时刻的差异，其中，脉冲数目的差异最为显著，所有鼻叶合开-耳廓柔性运动期间都只存在1个脉冲，而所有鼻叶随机-耳廓刚性运动期间都存在2个脉冲；脉冲激发时刻差异相对较小，每种耦合运动模式中脉冲激发时刻都比较分散，如图6(b)所示，图中灰色部分为两种运动重叠所致，黑色在下，白色在上。

图6 两种鼻叶/耳廓耦合运动模式的声呐脉冲序列Fig.6 Sonar pulse sequences for two noseleaf/pinna motion patterns

2.3 基于SSIMS的数据分析

2.3.1 SSIMS简介

神经脉冲模式分类(Spike train SIMilarity Space,SSIMS)是神经科学中对神经活动模式分类的方法之一，该方法由两部分组成：一是对神经脉冲序列间的距离度量，得到脉冲序列两两之间的相似度矩阵；二是利用非线性降维算法 t-SNE (t-distributed Stochastic Neighbor Embedding)，将高维的相似度矩阵投影至低维空间，以实现脉冲序列的可视化和分类。该方法几乎适用于所有类型的神经脉冲(时间序列)数据。

其中，脉冲序列间的距离度量是利用Victor和Purpura提出的距离度量(记为“V-P距离”[14])，该度量使用一个代价函数(cost function)量化通过脉冲添加、删除或时移等操作将一个脉冲序列Sa转变成另一个脉冲序列Sb所需要付出的代价，示意图如图7所示。同时赋予每种操作相应的代价值，其中，脉冲添加和删除的代价值均为1。脉冲时移的代价值为q×Δt，Δt为时移距离，q为时移参数，用来衡量脉冲时移的敏感度，当q非常大时，即使脉冲时移很小距离，时移操作的代价值会很大，即对脉冲时移非常敏感，而当q=0时，脉冲时移的代价值始终为0，即脉冲时移不产生任何代价，此时两个脉冲序列转变的代价仅来源于脉冲添加或删除。而V-P距离则是所有可能操作组合代价值的最小值，例如图7所示的将具有3个脉冲的序列Sa转变成具有3个脉冲的序列Sb，可以通过分别时移3个脉冲(即“时移-时移-时移”操作组合)，或删除其中一个或多个脉冲同时在相应位置添加脉冲，并时移剩余脉冲(如图7所示的“删除-添加-时移-时移”操作组合)的方法实现。但每种方法总的代价值不同，取决于移动参数q的选择，当q=2/0.3时，意味着当Δt ＜ 0 .3时，脉冲时移产生的代价值小于脉冲“删除-添加”组合的代价值(1+1=2)，此时应优先选用时移操作(图7序列Sa的第2和第3个脉冲)，而当Δt ＞ 0 .3时，优先选用“删除-添加”组合(图7序列Sa的第1个脉冲)，因此，图7所示两个脉冲序列间的V-P距离(相似度)为4。图7用来表征脉冲序列间的相似度，即从一个脉冲序列转变成另一个脉冲序列所需要的代价，代价通过脉冲删除、添加和移动所赋予的权重表征。图7中，α1表示一脉冲序列，有三个脉冲，激发时刻分别为0，0.8，1。α5该脉冲序列为另一个脉冲序列。α2～α4是中间的转换过程。

图7 脉冲序列间的距离度量示意图Fig.7 Schematic diagram of spike sequences’ distance metrics

2.3.2 声呐脉冲序列分析

利用 SSIMS方法对两种蝙蝠鼻叶/耳廓耦合运动模式对应的声呐脉冲序列进行分析，分析框图如图8所示，具体步骤如下：

图8 基于SSIMS方法的声呐脉冲序列分析框图Fig.8 Block diagram of SSIMS based sonar pulse analysis

(1) 随机选取34组鼻叶合开-耳廓柔性运动、34组鼻叶随机-耳廓刚性运动的声呐脉冲序列作为训练数据X1，剩余19组鼻叶合开-耳廓柔性运动、15组鼻叶随机-耳廓刚性运动的声呐脉冲序列作为测试数据X2；

(2) 计算训练数据X1中两两脉冲序列之间的V-P距离，得到相似度矩阵S1(68×68 矩阵)；

(3) 利用 t-SNE算法将相似度矩阵S1进行降维，得到低维矩阵L1(68×2矩阵)和转换矩阵 T1(68×2 矩阵)；

(4) 计算测试数据X2中每组声呐脉冲序列分别与训练数据中的声呐脉冲序列之间的V-P距离，得到相似度矩阵S2(34×68 矩阵)；

(5) 将相似度矩阵S2与转换矩阵 T1相乘，得到低维矩阵L2(34×2矩阵)，实现测试数据X2到低维空间的映射；

(6) 利用k近邻(k-th Nearest Neighbor, kNN)算法(k=3)，对低维矩阵L2分类，并与测试数据X2中声呐脉冲序列的真实耦合运动模式对比，以正确个数占总数的比例作为模式识别的正确率。

2.3.3 分析结果

利用 SSIMS方法对声呐脉冲序列模式识别的正确率取决于V-P距离的脉冲时移参数q，如图9所示。当q=0时，脉冲时移为无代价操作，此时脉冲序列间的相似度仅取决于序列中脉冲数目差异，模式识别正确率高达100%。这表明：在102组数据中，蝙蝠在所有鼻叶合开-耳廓柔性运动期间都只发射一个脉冲，在所有鼻叶随机-耳廓刚性运动期间都发射两个脉冲，这与从数据直接观察的结果(图 6)相符。当 q＞0时，脉冲时移操作开始要付出代价，即序列中脉冲激发时刻差异开始逐渐影响脉冲序列间的相似度，而此时模式识别正确率开始逐渐降低(图9中阶段 I)，这是由于两种鼻叶/耳廓耦合运动模式在脉冲激发时刻差异上并不十分显著，如图6所示。当4＜q＜36时，此时模式识别正确率整体开始上升且局部波动(图9中阶段Ⅱ)，这是由于脉冲激发时刻差异对相似度的影响开始小于脉冲数目差异对相似度的影响。当q足够大(约36)时，对大多数脉冲序列，时移操作的代价高于“删除-添加”组合操作的代价，序列间的相似度主要取决于脉冲数目差异，因此，此时模式识别正确率非常高(约100%)，如图9中的阶段Ⅲ。

图9 脉冲时移参数q对模式识别正确率的影响Fig.9 Effect of time shift q on pattern recognition accuracy

对测试数据低维投影L2的分类结果如图10所示。当q=0时，每种耦合运动模式的声呐脉冲序列都重合成一个点，如图10(a)所示，模式识别正确率为100%；当q=8时，虽然测试数据的数据点比较分散，模式识别正确率也能达到 100%，分类结果如图10(b)所示。

图10 生物声呐脉冲序列的分类结果Fig.10 Classification result of biosonar pulse sequences

3 结 论

本文首先通过生物实验采集普氏蹄蝠的鼻叶运动、耳廓运动及其同步的声呐脉冲序列；然后针对其中出现频率较高的两种鼻叶/耳廓耦合运动模式，即鼻叶合开-耳廓柔性运动和鼻叶随机-耳廓刚性运动，利用脉冲激发时刻定量表征声呐脉冲序列；最后利用 SSIMS方法对两种鼻叶/耳廓耦合运动模式的声呐脉冲序列分类。对两种鼻叶/耳廓耦合运动模式的102组声呐数据的分析结果表明：普氏蹄蝠声呐脉冲序列的变化与其声呐系统结构的动态形变相匹配，表现为不同鼻叶/耳廓耦合运动模式对应不同脉冲数目和脉冲激发时刻的声呐脉冲序列。这种耦合关系将蝙蝠声呐系统的动态特性提升至二维，因此，可利用两种动态特性的匹配关系改进动态仿生声呐系统的设计，推进蝙蝠动态传感信息编译机制的研究。