使用广义互相关进行源定位

打开实时脚本

这个例子展示了如何使用广义互相关(GCC)和三角测量来确定宽带信号源的位置。为了简单起见，这个示例仅限于一个二维场景，由一个源和两个接收传感器阵列组成。您可以将这种方法扩展到两个以上的传感器或传感器阵列，并扩展到三维。

简介

源定位不同于到达方向(DOA)估计。DOA估计仅寻求从传感器确定源的方向。源定位决定其位置。在本例中，源定位包括两个步骤，第一步是DOA估计。

使用DOA估计算法估计每个传感器阵列的源方向。对于宽带信号，许多著名的到达方向估计算法，如Capon方法或MUSIC，都不能应用，因为它们利用了元素之间的相位差，使得它们只适用于窄带信号。在宽带情况下，您可以使用元素之间信号到达时间的差异而不是相位信息。为了计算到达时间差异，本示例使用相变的广义互相关(GCC-PHAT)算法。根据到达时间的差异，可以计算出DOA。(有关窄带DOA估计算法的另一个示例，请参见高分辨率到达方向估计）.
通过三角剖分计算源位置。首先，从数组沿着到达方向画直线。然后，计算这两条直线的交点。这是源位置。源定位需要了解接收传感器或传感器阵列的位置和方向。

三角公式

三角剖分算法基于简单的三角公式。假设传感器阵列位于二维坐标中(0,0)而且(L, 0)而未知源的位置是(x, y)．通过了解传感器阵列的位置和到达阵列的两个方向， $θ_{1}$ 而且 $θ_{2}$ ，可以计算(x, y)坐标从

$l ＝ y 棕褐色 θ_{1} + y 棕褐色 θ_{2}$

你能解出哪个y

$y ＝ l / （棕褐色 θ_{1} + 棕褐色 θ_{2} ）$

然后对于x

$x ＝ y 棕褐色 θ_{1}$

本例的其余部分将展示如何使用相控阵系统工具箱™的函数和系统对象来计算源位置。

源和传感器几何

设置两个接收4元ula，沿x-轴，间隔50米。第一个阵列的相心为(0, 0, 0)．第二阵列相心为(50, 0, 0)．源位于(30100)米。如图所示，接收阵列的增益点在+ y方向。源在可能是方向。

指定传感器阵列之间的基线。

L = 50;

创建一个全向麦克风的4元素接收器ULA。你可以用同样的分阶段。齿龈的系统对象™分阶段。WidebandCollector而且分阶段。GCCEstimator两个数组的系统对象。

N = 4;rxULA =相控的。齿龈(“元素”,分阶段。OmnidirectionalMicrophoneElement,.．.“NumElements”N);

指定第一个传感器阵列的位置和方向。创建ULA时，数组元素会自动沿着y设在。必须将数组的局部轴旋转90°以使元素沿x全球坐标系的-轴。

Rxpos1 = [0;0;0];Rxvel1 = [0;0;0];Rxax1 = azelaxes(90,0);

指定第二个传感器阵列的位置和方向。选择第二个数组的局部轴与第一个数组的局部轴对齐。

rxpos2 = [L;0;0];Rxvel2 = [0;0;0];Rxax2 = rxax1;

指定信号源为单个全向换能器。

Srcpos = [30;100;0];Srcvel = [0;0;0];Srcax =杜鹃花(-90,0);srcULA = phase . omnidirectionalmicrophoneelement;

定义波形

选择源信号为宽带LFM波形。假设系统的工作频率为300 kHz，并将信号的带宽设置为100 kHz。假设最大操作距离为150米。然后，可以设置脉冲重复间隔PRI (pulse repeat interval)和脉冲重复频率PRF (pulse repeat frequency)。假设占空比为10%，并设置脉冲宽度。最后，在水下通道中使用音速1500米/秒。

设置LFM波形参数并创建分阶段。linearFMWaveform系统对象。

Fc = 300e3;% 300 kHzC = 1500;% 1500 m/sDmax = 150;% 150 mPri = (2*dmax)/c;PRF = 1/pri;Bw = 100.0e3;% 100 kHzFs = 2*bw;波形=相控。LinearFMWaveform (“SampleRate”fs,“SweepBandwidth”bw,.．.脉冲重复频率的脉冲重复频率,“脉冲宽度”革命制度党/ 10);

发射信号可以生成为

信号=波形();

辐射、传播和收集信号

宽带系统的辐射和传播建模比窄带系统的建模更复杂。例如，衰减取决于频率。多普勒频移以及因信号入射方向而引起的元件间的相移也随频率而变化。因此，在处理宽带信号时，对这些行为建模是至关重要的。本例使用子带方法。

设置子带数为128。

NFFT = 128;

指定源散热器和传感器阵列收集器。

散热器=阶段性。WidebandRadiator (“传感器”srcULA,.．.“PropagationSpeed”c“SampleRate”fs,.．.“CarrierFrequency”足球俱乐部,“NumSubbands”, nfft);Collector1 = phase。WidebandCollector (“传感器”rxULA,.．.“PropagationSpeed”c“SampleRate”fs,.．.“CarrierFrequency”足球俱乐部,“NumSubbands”, nfft);Collector2 = phase。WidebandCollector (“传感器”rxULA,.．.“PropagationSpeed”c“SampleRate”fs,.．.“CarrierFrequency”足球俱乐部,“NumSubbands”, nfft);

为从源到两个传感器阵列的路径创建宽带信号传播器。

通道1 =分阶段。WidebandFreeSpace (“PropagationSpeed”c.．.“SampleRate”fs,“OperatingFrequency”足球俱乐部,“NumSubbands”, nfft);Channel2 =分阶段。WidebandFreeSpace (“PropagationSpeed”c.．.“SampleRate”fs,“OperatingFrequency”足球俱乐部,“NumSubbands”, nfft);

确定从源到传感器阵列的传播方向。传播方向是相对于源的局部坐标系的。

[~，ang1t] = rangeangle(rxpos1,srcpos,srcax);[~，ang2t] = rangeangle(rxpos2,srcpos,srcax);

从信号源向传感器阵列的方向辐射信号。

Sigt =散热器(信号，[ang1t ang2t]);

然后，将信号传播到传感器阵列。

Sigp1 = channel1(sigt(:，1)，srcpos,rxpos1,srcvel,rxvel1);Sigp2 = channel2(sigt(:，2)，srcpos,rxpos2,srcvel,rxvel2);

计算传播信号在传感器阵列处的到达方向。由于收集器响应是传感器阵列局部坐标系中到达方向的函数，因此将局部坐标轴矩阵传递给rangeangle函数。

[~，ang1r] = rangeangle(srcpos,rxpos1,rxax1);[~，ang2r] = rangeangle(srcpos,rxpos2,rxax2);

在接收传感器阵列处收集信号。

Sigr1 = collector1(sigp1,ang1r);Sigr2 = collector2(sigp2,ang2r);

GCC估计和三角测量

创建GCC-PHAT估计器。

Doa1 =阶段性。GCCEstimator (“SensorArray”rxULA,“SampleRate”fs,.．.“PropagationSpeed”c);Doa2 =阶段性。GCCEstimator (“SensorArray”rxULA,“SampleRate”fs,.．.“PropagationSpeed”c);

估计到达的方向。

Angest1 = doa1(sigr1);Angest2 = doa2(sigr2);

使用前面建立的公式对源位置进行三角测量。因为这个场景局限于x - y平面，设z坐标为0。

= L/(abs(t (angest1)) + abs(t (angest2)));Xest = est*abs(tand(angest1));热情= 0;Srcpos_est = [xest;yest;zest]

srcpos_est =3×129.9881 100.5743 0

估计的震源位置与真实位置匹配，误差不超过30厘米。

总结

这个例子展示了如何使用三角测量来执行源定位。特别地，该示例演示了如何模拟、传播和处理宽带信号。GCC-PHAT算法用于估计宽带信号的到达方向。