技术领域
[0001] 本
发明涉及目标
定位技术领域，特别涉及一种基于二维正交非均匀线阵的水下波达方向估计方法及装置。
背景技术
[0002] 空间
信号的波达方向(简称：DOA)估计是阵列
信号处理研究中的重点。阵列信号处理是一种在空间的不同
位置放置多个
传感器组成的阵列，利用阵列对空间信号进行接收和处理，抑制无用的干扰和噪声，提取有用的信号特征或信息的技术；DOA估计便是一种利用这种技术来估计信号入射
角的方法。而在水面放置传感器阵列利用阵列信号处理技术来对水下目标物进行方位估计的方法就是水下DOA估计。
[0003] 现有的水下DOA估计方法主要有MUSIC
算法和ESPRIT算法。它们同属于特征子空间类算法，此类算法通过对阵列接收数据的数学分解(如奇异值分解、特征分解及QR分解等)，将接收数据划分为两个相互正交的子空间：一个与信号源的阵列流型空间一致的信号子空间和一个与信号子空间正交的噪声子空间，它利用两个子空间的正交性可以大大提高算法的分辨
力。特征子空间类算法从处理方式上分为两类：一类是以MUSIC算法为代表的噪声子空间类算法；另一类是便以ESPRIT算法为代表的信号子空间类算法。由于非均匀线阵相对于均匀线阵具有角度分辨力好，侧向
精度高，抗模糊性强的优势，并且非均匀线阵在阵元位置上可以随机放置，具有很强的灵活性，因此目前将MUSIC算法或者ESPRIT算法运用到非均匀线阵中的研究越来越多，同时也取得了较好的研究成果。但在这些基于非均匀线阵的DOA估计研究中，利用到的MUSIC算法或者ESPRIT算法都假定信号在介质中的传播速度是确定不变的，没有考虑速度对算法结果的影响，在水下环境复杂，速度变化较大、未知的情况下容易产生较大的误差，从而大大降低了目标定位的精确度。
发明内容
[0004] 本发明的主要目的在于克服
现有技术的缺点与不足，提供一种基于二维正交非均匀线阵的水下波达方向估计方法，为了克服ESPRIT算法在基于非均匀线阵的DOA估计中存在的缺点，本方法使用一个二维的正交非均匀线阵为接收阵列，实现对ESPRIT算法的改进，在算法过程中消去信号在介质中的传播速度，使最后的估计结果与速度无关，从而排除了速度对算法结果的影响，同时还利用了非均匀线阵的优势。因此，本方法在水下DOA估计中更具有实用性。
[0005] 本发明的另一目的在于提供一种基于二维正交非均匀线阵的水下波达方向估计装置。
[0006] 本发明的目的通过以下的技术方案实现：
[0007] 一种基于二维正交非均匀线阵的水下波达方向估计方法，包括以下步骤：
[0008] S1、推导水平非均匀线阵和竖直非均匀线阵的接收信号模型；
[0009] S2、结合接收信号模型推导水平非均匀线阵和竖直非均匀线阵的
相位补偿矩阵；
[0010] S3、求水平非均匀线阵和竖直非均匀线阵的初始相位补偿矩阵；
[0011] S4、对初始相位补偿后的水平接收信号阵列和竖直接收信号阵列进行
迭代，得到最终的水平旋转算子Φx和竖直旋转算子Φy的值；
[0012] S5、推导
声波在水平非均匀线阵上的方向角θx和在竖直非均匀线阵上的方向角θy之间的关系；
[0013] S6、根据声波在水平非均匀线阵上的方向角θx和在竖直非均匀线阵上的方向角θy之间的关系，以及最终的水平旋转算子Φx和竖直旋转算子Φy的值，求出方向角θx的最终大小。
[0014] 优选的，步骤S1具体步骤为：
[0015] 正交非均匀线阵中两条线阵所在直线的交点为原点O，窄带目标声源为S，中心
频率为f；水平和竖直非均匀线阵都有M个接收阵元，水平非均匀线阵第1阵元到第m个阵元相对于原点O的位置分别为x1,x2,…,xm,m＝1,2,…,M，其中 λ为声波的
波长，即水平非均匀线阵阵元之间的平均间距要小于声波信号的半波长；竖直非均匀线阵第1阵元到第m个阵元相对于原点O的位置分别为y1,y2,…,ym,m＝1,2,…,M，其中 即竖直非均匀线阵阵元之间的平均间距要小于声波信号的半波长；由于 而声波在探测路径上的速度v是未知的，因此取v为其范围中的最小值以确定λ的值，从而得到非均匀线阵阵元之间平均间距的范围；水平非均匀线阵各阵元之间的间距与竖直非均匀线阵各阵元之间的间距相等，即xm-xm-1＝ym-ym-1,m＝2,3,…,M；目标对应于水平线阵的方向角为θx，即从目标声源反射回来的声波与水平线阵的法线之间的夹角为θx；对应于竖直线阵的方向角为θy，即从目标声源反射回来的声波与竖直线阵的法线之间的夹角为θy；
[0016] 对于水平非均匀线阵，以处于位置x1处的第一个阵元为参考阵元，则第一个阵元接收的信号为：x1(t)＝sx(t)+nx1(t)，其中sx(t)表示t时刻到达水平非均匀线阵第一个阵元上的目标信号，即从目标反射回来的声波信号，nx1(t)表示水平非均匀线阵上第一个阵元上的噪声；
[0017] 接收信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个
相移；那么水平非均匀线阵上第m个阵元在同一时刻接收到的信号为：xm(t)＝sx(t)am(θx)+nxm(t),m＝1,2,…,M，其中 am(θx)中v表示声波在探测路径上的速度，nxm(t)表示水平非均匀线阵上第m个阵元上的噪声；
[0018] 于是，水平非均匀线阵的接收信号模型可表示为：
[0019] X(t)＝Axsx(t)+Nx(t) (1)
[0020] 其中， 为M×1的水平非均匀线阵的导向矢量矩阵，X(t)＝[x1(t),x2(t),…,xM(t)]T为M×1的水平非均匀线阵的接收信号矩阵，Nx(t)＝[nx1(t),nx2(t),…,nxM(t)]T为M×1的水平非均匀线阵的噪声矩阵；
[0021] 对于竖直非均匀线阵，以处于位置y1处的第一个阵元为参考阵元，则第一个阵元接收的信号为：y1(t)＝sy(t)+ny1(t)，其中sy(t)表示t时刻到达竖直非均匀线阵第一个阵元上的目标信号，即从目标反射回来的声波信号，其中ny1(t)表示竖直非均匀线阵上第一个阵元上的噪声；那么竖直非均匀线阵上第m个阵元在同一时刻接收到的信号为：ym(t)＝sy(t)am(θy)+nym(t),m＝1,2,…,M，其中 am(θy)中v表示声波在探测路径上的速度，nym(t)表示竖直非均匀线阵上第m个阵元上的噪声；
[0022] 于是，竖直非均匀线阵的接收信号模型可表示为：
[0023] Y(t)＝Aysy(t)+Ny(t) (2)
[0024] 其中， 为M×1的竖直非均匀线阵的导向矢量矩阵，Y(t)＝[y1(t),y2(t),…,yM(t)]T为M×1的竖直非均匀线阵的接收信号矩阵，Ny(t)＝[ny1(t),ny2(t),…,nyM(t)]T为M×1的竖直非均匀线阵的噪声矩阵。
[0025] 优选的，步骤S2具体步骤为：
[0026] 对于水平非均匀线阵，线阵X是实际水平非均匀的M元线阵，线阵 是水平非均匀线阵对应的参考水平均匀M元线阵，其中参考水平均匀线阵第一阵元到第m个阵元相对于原点的位置分别为 参考水平均匀线阵第一个阵元和最后一个阵元的位置与实际水平非均匀线阵第一个阵元和最后一个阵元的位置相同，即：
并且参考水平均匀线阵的阵元间距为 则实际水平非均匀线阵
和参考水平均匀线阵之间的位置差Δx可表示为：
[0027]
[0028] 其中，
[0029] 那么，参考水平均匀线阵的导向矩阵可以表示为
[0030] 实际水平非均匀线阵和参考水平均匀线阵阵元位置的偏差引起的相位偏移Px(θx)可以表示为：
[0031]
[0032] 从而可得实际水平非均匀线阵的导向矩阵与参考水平均匀线阵的导向矩阵之间的关系为：
[0033]
[0034] 定义θx方向上的水平相位补偿矩阵Bx(θx)为Px(θx)的逆矩阵，即水平相位补偿矩阵可表示为：
[0035]
[0036] 将式(6)代入式(5)，并
整理可得：
[0037]
[0038] 对于实际水平非均匀线阵接收到的数据，有X(t)＝Axsx(t)+Nx(t)，那么对其进行相位补偿后得到的参考均匀线阵接收到的数据有：
[0039]
[0040] 对于竖直非均匀线阵，线阵Y是实际竖直非均匀的M元线阵，线阵 是竖直非均匀线阵对应的参考竖直均匀M元线阵，其中参考竖直均匀线阵第一阵元到第m个阵元相对于原点的位置分别为 参考竖直均匀线阵第一个阵元和最后一个阵元的位置与实际竖直非均匀线阵第一个阵元和最后一个阵元的位置相同，即：
并且参考竖直均匀线阵的阵元间距 则实际竖直非均匀阵列和参考竖直均匀阵列之间的位置差可表示为：
[0041]
[0042] 其中，
[0043] 那么，参考均匀线阵的导向矩阵为
[0044] 实际竖直非均匀线阵和参考竖直均匀线阵阵元位置的偏差引起的相位偏移Py(θy)可以表示为：
[0045]
[0046] 从而可得实际竖直非均匀线阵的导向矩阵与参考竖直均匀线阵的导向矩阵之间的关系为：
[0047]
[0048] 定义θy方向上的竖直相位补偿矩阵By(θy)为Py(θy)的逆矩阵，即竖直相位补偿矩阵可表示为：
[0049]
[0050] 将式(12)代入式(11)，并整理可得：
[0051]
[0052] 对于实际竖直非均匀线阵接收到的数据，有Y(t)＝Aysy(t)+Ny(t)，那么对其进行相位补偿后得到的参考均匀线阵接收到的数据有：
[0053]
[0054] 优选的，步骤S3具体步骤为：
[0055] 对于实际水平非均匀线阵的接收信号矩阵X(t)，直接用ESPRIT算法对X(t)进行估计，可以得到水平旋转算子的值Φx，其中 从而可得到第一次估计后θx的值：
[0056]
[0057] 其中 θx1表示第一次估计后得到的θx值，Φx1表示第一次估计后得到的水平旋转算子Φx值；
[0058] 将式(15)代入式(4)，可得：
[0059]
[0060] 其中，Px1(θx1)表示第一估计后的Px(θx)值；
[0061] 因此，水平初始相位补偿矩阵为：
[0062] 于是，经过初始相位补偿后的参考水平均匀线阵的接收信号阵列为：
[0063]
[0064] 对于实际竖直非均匀线阵的接收信号矩阵Y(t)，直接用ESPRIT算法对Y(t)进行估计，可以得到竖直旋转算子的值Φy，其中 从而可得到第一次估计后θy的值：
[0065]
[0066] 其中 θy1表示第一次估计后得到的θy值，Φy1表示第一次估计后得到的竖直旋转算子Φy值；
[0067] 将式(17)代入式(10)，可得：
[0068]
[0069] 其中，Py1(θy1)表示第一估计后的Py(θy)值；
[0070] 因此，竖直初始相位补偿矩阵为：
[0071] 于是，经过初始相位补偿后的参考竖直均匀线阵的接收信号阵列为：
[0072]
[0073] 优选的，步骤S4具体步骤为：
[0074] 因为当相位补偿矩阵是基于每个接近于真实信号的角度时，所得角度将收敛于该真实信号的角度，而由公式(15)和公式(17)知，当角度确定时，旋转算子也将确定，因此在迭代过程中，每次估计的旋转算子的值也将收敛；因此，对于初始相位补偿后的水平均匀线阵的接收阵列 再次利用ESPRIT算法对其进行DOA估计，得到第二次估计后的水平旋转算子的值Φx2，判断 是否小于ε，其中ε是一个根据实际情况而定的
阈值，若小于ε，则Φx2就是真实角度θx对应的水平旋转算子Φx；若不小于ε，则利用Φx2继续根据步骤S3的步骤更新相位补偿矩阵并估计出第二次相位补偿后的均匀线阵接收信号阵列然后对其进行第三次估计得到水平旋转算子Φx3，判断其是否收敛，若不收敛，则继续迭代，直至当Φxt满足以下判断准则时:
[0075]
[0076] 其中Φxt表示第t次估计后得到的水平旋转算子的值，ε是一个根据实际情况而定的阈值；则有Φx＝Φxt，于是声波相对于水平非均匀线阵的方向角为：
[0077]
[0078] 对于初始相位补偿后的竖直均匀线阵的接收阵列 同样利用ESPRIT算法对其进行估计和迭代，当Φyt满足以下判断准则时:
[0079]
[0080] 其中Φyt表示第t次估计后得到的竖直旋转算子的值，ε是一个根据实际情况而定的阈值；则有Φy＝Φyt，于是声波相对于竖直非均匀线阵的方向角为：
[0081]
[0082] 优选的，步骤S5具体步骤为：
[0083] 以两个正交非均匀线阵所在的方向建立一个直角
坐标系，因为只是在水面的位置对水面下方的目标进行定位，因此可只考虑从实际水平非均匀线阵下方反射回来的声波；以反射回来的声波在实际水平非均匀线阵上第一个到达的阵元位置可将声波分为两种：一种是声波到达的第一个阵元是参考阵元，如声波①；另一种是声波到达的第一个阵元不是参考阵元，如声波②；设θxA和θyA分别为声波①与两个非均匀线阵法线的夹角，它们相加之和为 θxB和θyB分别为声波②与两个均匀线阵法线的夹角，它们的相加之和也为[0084] 当声波在水平非均匀线阵上到达的第一个阵元是参考阵元，即属于声波①的情形时，有θx＝θxA，θy＝θyA，由此可得θx和θy的关系为：
[0085]
[0086] 当声波在水平非均匀线阵上到达的第一个阵元不是参考阵元，即属于声波②的情形时，则参考阵元x1是最晚接收到信号的，因此时延参数τ是小于0的，又因为所以θx＜0，于是有θx＝-θxB，θy＝θyB，由此可得θx和θy的关系为：
[0087]
[0088] 根据公式(23)和公式(24)可以得到：
[0089]
[0090] 优选的，步骤S6具体步骤为：
[0091] 根据步骤五中得到的θx和θy之间的关系，即公式(25)，将它与公式(20)和公式(22)结合，又因为dx＝dy，从而可以得到θx的最终表达式：
[0092]
[0093] 当arg(Φx)大于0时，θx取正值，即目标声源在水平非均匀线阵参考阵元的左侧；当arg(Φy)小于0时，θx取负值，即目标声源在水平非均匀线阵参考阵元的右侧。
[0094] 一种基于非均匀正交线阵的水下波达方向估计装置，包括控
制模块、处理模块、AD转换模块、发射模块、接收模块、输出模块和电源模块；
[0095]
控制模块与处理模块相连，能够根据应用场景的不同而选择不同的参数，包括发送信号中心频率、脉冲长度，并将这些参数传送给处理模块，控制处理模块发送不同的信号；
[0096] 处理模块与AD转换模块、输出模块相连，能够根据控制模块输入的指令生成特定的
数字信号，该数字信号通过AD转换模块转换为
模拟信号；同时，它还能够对经过AD转换模块转换后的数字信号进行分析处理，从而得到声波信号的方向角；处理模块还能够将处理好的数据存储并且传送给输出模块，使数据传送到外部装置；
[0097] AD转换模块与发射模块、接收模块相连，A/D转换器将接收模块接收到的模拟信号转换为数字信号并传递给处理模块进行处理；D/A转换器将处理模块发出的数字信号转换为模拟信号并传送发射模块使其发射
指定的声波；
[0098] 电源模块与控制模块、处理模块、AD转换模块、发射模块、接收模块、输出模块相连，为这些模块供电。
[0099] 优选的，发射模块包括一个阻抗匹配
电路和一个
超声波发送
探头。
[0100] 优选的，接收模块包括两个以相互垂直方式放置的非均匀超声波接收探头线阵和阻抗匹配电路。
[0101] 优选的，输出模块包括有线/无线
接口和显示屏；有线/无线接口能够将处理模块中处理好的最终结果传送到外部设备或系统中去；显示屏作为直接的显示工具，能够将处理得到的最终结果在显示屏上进行显示。
[0102] 本发明与现有技术相比，具有如下优点和有益效果：
[0103] 1、本发明方法利用了非均匀线阵对于波达方向估计的优势，提高估计结果的角度分辨力，侧向精度和抗模糊性，同时还使超声波接收探头的摆放具有很强的灵活性。因此本发明算法有着很好的估计性能和很强的操作方便性。
[0104] 2、本发明相对于直接利用MUSIC算法或者ESPRIT算法进行水下波达方向估计的方法更具实用性，并且精确度也更高。在复杂的水声环境中，每一条探测路径上的声
波速度并不确定，即目标声源反射回来的声波速度是未知的，而传统MUSIC算法和ESPRIT算法假定声速是不变，在每条路径上都是相同的，没有考虑声速的不确定性对估计算法的影响，这使得在复杂水声环境下MUSIC算法和ESPRIT算法不具有实用性。本发明采用了两条相互正交的非均匀线阵，对ESPRIT算法进行了改进，在算法过程中消去了声速，使得最后的估计结果
和声速无关，从而排除了声速的不确定对波达方向估计的影响，可在复杂水声环境下进行测量，使得其实用性更高，并且目标定位的精度也更高。因此，本发明方法实用性强，且精度更高。
[0105] 3、本发明方法的运算量和复杂度相对来说都比较适中，从而保证了算法的可行性。
[0106] 4、本发明装置可行性强，安装简单。
附图说明
[0107] 图1为
实施例装置的模块连接组合示意图。
[0108] 图2为实施例具体装置连接组合示意图。
[0109] 图3为实施例正交非均匀线阵模型场景。
[0110] 图4为实施例水平非均匀线阵的接收信号示意图。
[0111] 图5为实施例水平非均匀线阵以及对应的参考水平均匀线阵示意图。
[0112] 图6为实施例竖直非均匀线阵以及对应的参考竖直均匀线阵示意图。
[0113] 图7为实施例正交非均匀线阵的直角坐标系示意图。
[0115] 图9为实施例装置具体工作步骤流程图。
具体实施方式
[0116] 下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
[0117] 实施例1
[0118] 一种基于二维正交非均匀线阵的水下波达方向估计方法，本方法的原理是先对非均匀线阵的导向矩阵进行相位补偿，得到非均匀线阵对应的均匀线阵的导向矩阵，然后利用ESPRIT算法对此均匀线阵的导向矩阵进行DOA估计迭代，得到两个非均匀线阵所对应的均匀线阵上的旋转算子，再利用从目标反射回来的声波在非均匀正交线阵上形成的方向角之间关系，推导出含有两个旋转算子且不含声速的波达方向角的表达式，最后将两个旋转算子的值代入表达式中从而求得最终方向角的大小。
[0119] 本方法的具体步骤包括以下几步：
[0120] 步骤一：推导水平非均匀线阵和竖直非均匀线阵的接收信号模型。
[0121] 正交非均匀线阵模型场景如图3所示，两条线阵所在直线的交点为原点O。窄带目标声源为S，中心频率为f。水平和竖直非均匀线阵都有M个接收阵元，水平非均匀线阵第1阵元到第m个阵元相对于原点O的位置分别为x1,x2,…,xm,m＝1,2,…,M，其中 λ为声波的波长，即水平非均匀线阵阵元之间的平均间距要小于声波信号的半波长；竖直非均匀线阵第1阵元到第m个阵元相对于原点O的位置分别为y1,y2,…,ym,m＝1,2,…,M，其中即竖直非均匀线阵阵元之间的平均间距要小于声波信号的半波长；由于而声波在探测路径上的速度v是未知的，因此取v为其范围中的最小值以确定λ的值，从而得到非均匀线阵阵元之间平均间距的范围。水平非均匀线阵各阵元之间的间距与竖直非均匀线阵各阵元之间的间距相等，即xm-xm-1＝ym-ym-1,m＝2,3,…,M。目标对应于水平线阵的方向角为θx，即从目标声源反射回来的声波与水平线阵的法线之间的夹角为θx；对应于竖直线阵的方向角为θy，即从目标声源反射回来的声波与竖直线阵的法线之间的夹角为θy。
[0122] 对于水平非均匀线阵，水平非均匀线阵接收信号的示意图如图4所示，以处于位置x1处的第一个阵元为参考阵元，则第一个阵元接收的信号为：x1(t)＝sx(t)+nx1(t)，其中sx(t)表示t时刻到达水平非均匀线阵第一个阵元上的目标信号，即从目标反射回来的声波信号，nx1(t)表示水平非均匀线阵上第一个阵元上的噪声。
[0123] 接收信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。那么水平非均匀线阵上第m个阵元在同一时刻接收到的信号为：xm(t)＝sx(t)am(θx)+nxm(t),m＝1,2,…,M，其中 am(θx)中v表示声波在探测路径上的速度，nxm(t)表示水平非均匀线阵上第m个阵元上的噪声。
[0124] 于是，水平非均匀线阵的接收信号模型可表示为：
[0125] X(t)＝Axsx(t)+Nx(t) (1)
[0126] 其中， 为M×1的水平非均匀线阵的导向矢量矩阵，X(t)＝[x1(t),x2(t),…,xM(t)]T为M×1的水平非均匀线阵的接收信号矩阵，Nx(t)＝[nx1(t),nx2(t),…,nxM(t)]T为M×1的水平非均匀线阵的噪声矩阵。
[0127] 对于竖直非均匀线阵，以处于位置y1处的第一个阵元为参考阵元，则第一个阵元接收的信号为：y1(t)＝sy(t)+ny1(t)，其中sy(t)表示t时刻到达竖直非均匀线阵第一个阵元上的目标信号，即从目标反射回来的声波信号，其中ny1(t)表示竖直非均匀线阵上第一个阵元上的噪声。那么竖直非均匀线阵上第m个阵元在同一时刻接收到的信号为：ym(t)＝sy(t)am(θy)+nym(t),m＝1,2,…,M，其中 am(θy)中v表示声波在探测路径上的速度，nym(t)表示竖直非均匀线阵上第m个阵元上的噪声。
[0128] 于是，竖直非均匀线阵的接收信号模型可表示为：
[0129] Y(t)＝Aysy(t)+Ny(t) (2)
[0130] 其中， 为M×1的竖直非均匀线阵的导向矢量矩阵，Y(t)＝[y1(t),y2(t),…,yM(t)]T为M×1的竖直非均匀线阵的接收信号矩阵，Ny(t)＝[ny1(t),ny2(t),…,nyM(t)]T为M×1的竖直非均匀线阵的噪声矩阵。
[0131] 步骤二：推导水平非均匀线阵和竖直非均匀线阵的相位补偿矩阵。
[0132] 对于水平非均匀线阵，如图5所示，线阵X是实际水平非均匀的M元线阵，线阵 是水平非均匀线阵对应的参考水平均匀M元线阵，其中参考水平均匀线阵第一阵元到第m个阵元相对于原点的位置分别为 参考水平均匀线阵第一个阵元和最后一个阵元的位置与实际水平非均匀线阵第一个阵元和最后一个阵元的位置相同，即：
并且参考水平均匀线阵的阵元间距为 则实际水平非均匀线阵
和参考水平均匀线阵之间的位置差Δx可表示为：
[0133]
[0134] 其中，
[0135] 那么，参考水平均匀线阵的导向矩阵可以表示为
[0136] 实际水平非均匀线阵和参考水平均匀线阵阵元位置的偏差引起的相位偏移Px(θx)可以表示为：
[0137]
[0138] 从而可得实际水平非均匀线阵的导向矩阵与参考水平均匀线阵的导向矩阵之间的关系为：
[0139]
[0140] 定义θx方向上的水平相位补偿矩阵Bx(θx)为Px(θx)的逆矩阵，即水平相位补偿矩阵可表示为：
[0141]
[0142] 将式(6)代入式(5)，并整理可得：
[0143]
[0144] 对于实际水平非均匀线阵接收到的数据，有X(t)＝Axsx(t)+Nx(t)，那么对其进行相位补偿后得到的参考均匀线阵接收到的数据有：
[0145]
[0146] 因为 是参考水平均匀线阵的接收信号矩阵，即接收数据，而参考均匀线阵具有移不变性，所以可以用ESPRIT算法对接收信号矩阵 进行DOA估计，从而可求出对应于水平线阵的方向角θx。
[0147] 对于竖直非均匀线阵，如图6所示，线阵Y是实际竖直非均匀的M元线阵，线阵 是竖直非均匀线阵对应的参考竖直均匀M元线阵，其中参考竖直均匀线阵第一阵元到第m个阵元相对于原点的位置分别为 参考竖直均匀线阵第一个阵元和最后一个阵元的位置与实际竖直非均匀线阵第一个阵元和最后一个阵元的位置相同，即：
并且参考竖直均匀线阵的阵元间距 则实际竖直非均匀阵列
和参考竖直均匀阵列之间的位置差可表示为：
[0148]
[0149] 其中，
[0150] 那么，参考均匀线阵的导向矩阵为
[0151] 实际竖直非均匀线阵和参考竖直均匀线阵阵元位置的偏差引起的相位偏移Py(θy)可以表示为：
[0152]
[0153] 从而可得实际竖直非均匀线阵的导向矩阵与参考竖直均匀线阵的导向矩阵之间的关系为：
[0154]
[0155] 定义θy方向上的竖直相位补偿矩阵By(θy)为Py(θy)的逆矩阵，即竖直相位补偿矩阵可表示为：
[0156]
[0157] 将式(12)代入式(11)，并整理可得：
[0158]
[0159] 对于实际竖直非均匀线阵接收到的数据，有Y(t)＝Aysy(t)+Ny(t)，那么对其进行相位补偿后得到的参考均匀线阵接收到的数据有：
[0160]
[0161] 因为 是参考竖直均匀线阵的接收信号矩阵，即接收数据，而参考均匀线阵具有移不变性，所以可以用ESPRIT算法对竖直均匀线阵的接收信号矩阵 进行DOA估计，同样可求出对应于竖直线阵的方向角θy。
[0162] 步骤三：求水平非均匀线阵和竖直非均匀线阵的初始相位补偿矩阵。
[0163] 因为相位补偿矩阵中含有入射波的方向角，而入射波的方向角一开始是未知的，所以需要先计算出方向角的一个近似值来构造初始相位补偿矩阵，从而得到近似的参考均匀线阵的接收信号矩阵。采用ESPRIT算法直接对实际的非均匀线阵进行DOA估计，可以得到真实方向角的一个近似值。
[0164] 对于实际水平非均匀线阵的接收信号矩阵X(t)，直接用ESPRIT算法对X(t)进行估计，可以得到水平旋转算子的值Φx，其中 从而可得到第一次估计后θx的值：
[0165]
[0166] 其中 θx1表示第一次估计后得到的θx值，Φx1表示第一次估计后得到的水平旋转算子Φx值。
[0167] 将式(15)代入式(4)，可得：
[0168]
[0169] 其中，Px1(θx1)表示第一估计后的Px(θx)值。
[0170] 因此，水平初始相位补偿矩阵为：
[0171] 于是，经过初始相位补偿后的参考水平均匀线阵的接收信号阵列为：
[0172]
[0173] 对于实际竖直非均匀线阵的接收信号矩阵Y(t)，直接用ESPRIT算法对Y(t)进行估计，可以得到竖直旋转算子的值Φy，其中 从而可得到第一次估计后θy的值：
[0174]
[0175] 其中 θy1表示第一次估计后得到的θy值，Φy1表示第一次估计后得到的竖直旋转算子Φy值。
[0176] 将式(17)代入式(10)，可得：
[0177]
[0178] 其中，Py1(θy1)表示第一估计后的Py(θy)值。
[0179] 因此，竖直初始相位补偿矩阵为：
[0180] 于是，经过初始相位补偿后的参考竖直均匀线阵的接收信号阵列为：
[0181]
[0182] 步骤四：对初始相位补偿后的水平接收信号阵列和竖直接收信号阵列进行迭代，得到最终的水平旋转算子Φx和竖直旋转算子Φy的值。
[0183] 因为当相位补偿矩阵是基于每个接近于真实信号的角度时，所得角度将收敛于该真实信号的角度，而由公式(15)和公式(17)知，当角度确定时，旋转算子也将确定，因此在迭代过程中，每次估计的旋转算子的值也将收敛。因此，对于初始相位补偿后的水平均匀线阵的接收阵列 再次利用ESPRIT算法对其进行DOA估计，得到第二次估计后的水平旋转算子的值Φx2，判断 是否小于ε，其中ε是一个根据实际情况而定的阈值，若小于ε，则Φx2就是真实角度θx对应的水平旋转算子Φx；若不小于ε，则利用Φx2继续根据步骤三的步骤更新相位补偿矩阵并估计出第二次相位补偿后的均匀线阵接收信号阵列 然后对其进行第三次估计得到水平旋转算子Φx3，判断其是否收敛，若不收敛，则继续迭代，直至当Φxt满足以下判断准则时:
[0184]
[0185] 其中Φxt表示第t次估计后得到的水平旋转算子的值，ε是一个根据实际情况而定的阈值。则有Φx＝Φxt，于是声波相对于水平非均匀线阵的方向角为：
[0186]
[0187] 对于初始相位补偿后的竖直均匀线阵的接收阵列 同样利用ESPRIT算法对其进行估计和迭代，当Φyt满足以下判断准则时:
[0188]
[0189] 其中Φyt表示第t次估计后得到的竖直旋转算子的值，ε是一个根据实际情况而定的阈值。
[0190] 则有Φy＝Φyt，于是声波相对于竖直非均匀线阵的方向角为：
[0191]
[0192] 步骤五：推导声波在水平非均匀线阵上的方向角θx和在竖直非均匀线阵上的方向角θy之间的关系。
[0193] 以两个正交非均匀线阵所在的方向建立一个直角坐标系，如图7所示。因为只是在水面的位置对水面下方的目标进行定位，因此可只考虑从实际水平非均匀线阵下方反射回来的声波。以反射回来的声波在实际水平非均匀线阵上第一个到达的阵元位置可将声波分为两种：一种是声波到达的第一个阵元是参考阵元，如图7中的声波①；另一种是声波到达的第一个阵元不是参考阵元，如图7中的声波②。设θxA和θyA分别为声波①与两个非均匀线阵法线的夹角，它们相加之和为 θxB和θyB分别为声波②与两个均匀线阵法线的夹角，它们的相加之和也为
[0194] 当声波在水平非均匀线阵上到达的第一个阵元是参考阵元，即属于声波①的情形时，有θx＝θxA，θy＝θyA，由此可得θx和θy的关系为：
[0195]
[0196] 当声波在水平非均匀线阵上到达的第一个阵元不是参考阵元，即属于声波②的情形时，则参考阵元x1是最晚接收到信号的，因此时延参数τ是小于0的，又因为所以θx＜0，于是有θx＝-θxB，θy＝θyB，由此可得θx和θy的关系为：
[0197]
[0198] 根据公式(23)和公式(24)可以得到：
[0199]
[0200] 步骤六：根据声波在水平非均匀线阵上的方向角θx和在竖直非均匀线阵上的方向角θy之间的关系求出方向角θx的最终大小。
[0201] 虽然在步骤四中求出θx和θy的表达式，即公式(20)和公式(22)。但由于公式中的速度是未知的，因此仅仅求出了参考均匀线阵的旋转算子Φx和Φy，并不能计算出最终的方向角θx和θy。
[0202] 根据步骤五中得到的θx和θy之间的关系，即公式(25)，将它与公式(20)和公式(22)结合，又因为dx＝dy，从而可以得到θx的最终表达式：
[0203]
[0204] 当arg(Φx)大于0时，θx取正值，即目标声源在水平非均匀线阵参考阵元的左侧；当arg(Φy)小于0时，θx取负值，即目标声源在水平非均匀线阵参考阵元的右侧。
[0205] 在最终的结果公式(26)中，没有出现声波在探测路径上的速度v，因此消除了由于声速变化过大、未知对结果带来的影响，进一步提高了精度。
[0206] 根据以上算法流程可知，本方法提出的改进算法不需要知道声速的大小就可以对波达方向角θx进行精确的估计，克服了传统ESPRIT算法的缺点，同时还利用了非均匀线阵的优势。
[0207] 以上的算法流程图可以由图8表示。
[0208] 实施例2
[0209] 一种基于非均匀正交线阵的水下波达方向估计装置如图1所示，包括控制模块、处理模块、AD转换模块、发射模块、接收模块、输出模块和电源模块。
[0210] 控制模块由一个PC机组成，并且与处理模块和电源模块相连，主要用于控制处理模块发送不同的信号。它能够根据应用场景的不同而选择不同的发送信号中心频率、脉冲长度等参数，并将这些参数传送给处理模块。
[0211] 处理模块由一个
数字信号处理器组成，并且与控制模块、AD转换模块、输出模块、电源模块相连。它能够根据控制模块输入的指令生成特定的数字信号，该数字信号通过AD转换模块转换为模拟信号，最后通过发射模块发射出指定频率的声波信号；同时，它还能够对从接收模块传送回来，并经过AD转换模块转换后的数字信号进行分析处理，从而得到声波信号的方向角。除此之外，处理模块还能够将处理好的数据存储并且传送给输出模块，使数据传送到外部装置。
[0212] AD转换模块由一个A/D转换器和D/A转换器组成，并且与处理模块、发射模块、接收模块和电源模块相连。A/D转换器将接收模块接收到的模拟信号转换为数字信号并传递给处理模块进行处理；D/A转换器将处理模块发出的数字信号转换为模拟信号并传送发射模块使其发射指定的声波。
[0213] 发射模块由一个阻抗匹配电路和一个超声波发送探头组成，并且与处理模块和电源模块相连。它能够根据处理模块发出，并经过AD转换模块转换后的模拟信号发射指定频率的超声波信号；
[0214] 接收模块由两个以相互垂直方式放置的非均匀超声波接收探头线阵组成，同时也包含阻抗匹配电路，并与AD转换模块和电源模块相连。它能够接收从目标反射回来的超声波信号，并通过AD转换模块将其转换为数字信号然后传送给处理模块进行分析处理。
[0215] 输出模块由有线/无线接口和LCD显示屏组成，并且与处理模块和电源模块相连。有线/无线接口能够将处理模块中处理好的最终结果传送到外部设备或系统中去；LCD显示屏作为直接的显示工具，能够将处理得到的最终结果在LCD上进行显示。
[0216] 电源模块由一个电源组成，并且与控制模块、处理模块、AD转换模块、发射模块、接收模块、输出模块相连。它能够为这些模块供电。
[0217] 本发明装置的主要工作流程如下：操作者通过控制模块输入指定的参数，使处理模块产生所需要的数字信号，该数字信号经过AD转换模块将其转换成模拟信号后，传送给发射模块，超声波发射探头将对其进行发射。接收模块中的非均匀超声波接收探头线阵接收到反射回来的声波信号并将其传送给AD转换模块，信号转换为数字信号后传送给处理模块，模块中的数字信号处理器根据实施例1提出的算法对这些数字信号进行处理，从而获得最终的波达方向角。最后，处理模块将分析处理后的数据结果通过输出模块传送给外部设备或系统，并将数据结果在LCD显示屏上显示出来。工作的全程中，电源模块为所有模块供电。
[0218] 图1为本装置的模块连接
框图，该装置主要包括控制模块101、处理模块102、AD转换模块103、发射模块104、接收模块105、输出模块106、电源模块107七个模块。具体的装置连接如图2所示，控制模块101由PC机110实现；处理模块102中的数字信号处理器201可用DSP芯片(如：TI公司的TMS320VC5509A型号的DSP芯片)或者ARM芯片等实现；AD转换模块103中的D/A转换器301和A/D转换器302可用A/D转换器实现；发射模块104由阻抗匹配电路401和一个超声波发射探头402组成；接收模块105由阻抗匹配电路501和两个相互垂直的非均匀超声波接收探头阵列502组成；输出模块106由有线/无线接口601和一块LCD显示屏602组成；电源模块107由电源组成。
[0219] 本装置的主要工作步骤如图9所示，具体如下：
[0220] 步骤1：按图2连接好具体装置，其中接收模块中的每个非均匀线阵中的阵元个数M定为8。利用操作模块给处理模块发送指令，控制超声发射探头发射单频超声信号s(t)，信号的频率为fs＝10kHz，脉冲长度5ms；
[0221] 步骤2：对超声接收探头线阵接收到的目标声源信号进行
采样；水平方向非均匀线阵接收到的信号为x1(t),x2(t),…,x8(t)，竖直方向非均匀线阵接收的信号为y1(t),y2(t),…,y8(t)。共采样接收200次，并将接收到的信号传递给处理模块进行分析处理。
[0222] 步骤3：信号在处理模块中的分析处理步骤具体如下：
[0223] 1)根据接收到的信号得出水平非均匀线阵的接收信号矩阵X(t)和竖直非均匀线阵的接收信号矩阵Y(t)，随后直接用EPSRIT算法对X(t)和Y(t)进行估计，得到第一次估计后的水平非均匀线阵的水平旋转算子Φx1和竖直非均匀线阵的竖直旋转算子Φy1。
[0224] 2)利用Φx1和Φy1分别得到水平初始相位补偿矩阵Bx1(θx1)和竖直初始相位补偿矩阵By1(θy1)，从而计算出经过初始相位补偿后的参考水平均匀线阵的接收信号矩阵 和参考竖直均匀线阵的接收信号矩阵
[0225] 3)利用ESPRIT算法对 和 进行估计，得到第二次估计后的水平旋转算子Φx2和竖直旋转算子Φy2。
[0226] 4)判断前后两次估计的旋转算子是否收敛，即是否满足 和其中，阈值ε的值设定为0.05，若满足，则Φx＝Φxt且Φy＝Φyt；若不满足，重复步骤2至步骤4。
[0227] 5)利用估计出的最终旋转算子Φx和Φy，计算出最终的波达方向角[0228] 步骤4：将计算出的波达方向角信息存储下来，并通过输出模块中的有线或者无线接口发送给外部设备或者系统，同时在LCD屏上显示出来。
[0229] 上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
