地球的大部分区域被海洋覆盖，探索海洋具有重大意义。在实现对海洋的有效监测方面，声呐设备是一种有效的工具^[1-3]，根据声呐的工作模式，可将声呐设备分为被动声呐和主动声呐。被动声呐只被动地接收信号，具有很好的隐蔽性，但随着消声瓦和减振技术的应用，水下目标的辐射噪声和反射强度显著降低，这导致被动声呐的有效探测范围不断减小。主动声呐通过发射声波并接收目标反射的回波信号进行检测和参数估计来探测目标，虽然不依赖目标的噪声辐射水平，但是自发自收的工作模式存在隐蔽性差、灵活性低和有效探测距离受限等问题。为了弥补被动声呐和主动声呐存在的不足，收发分置的双基地声呐系统被提出，并成为当前热门的水下应用平台之一。

水下目标定位是双基地声呐实现水下安全航行、快速打捞救生、水下目标探测及海洋资源开发等应用的基础^[4-6]。目标定位基于获得的各种测量物理量构建相应的定位模型，然后使用不同类型的方法得到目标位置。常用的物理量包括到达角(Angle of arrival，AoA)、到达时间(Time of arrival，ToA)、接收信号强度(Received signal strength，RSS)、到达时间差(Time difference of arrival，TDoA)和多普勒因子等。由于声信号在水下传播的速度较慢，基站与目标的相对运动会引起显著的多普勒效应。基于此，ToA和多普勒因子已被广泛应用于水下目标定位^[6-8]。在联合ToA和多普勒因子进行水下目标定位时，以上两个物理量的获取有两种方式：实测和仿真；由于水下实测数据获取困难，因此在验证定位算法性能时通常采用仿真数据。这两个物理量的精度是决定目标定位精度的关键因素，而定位精度通常与协方差矩阵有关。为了得到测量协方差矩阵，传统方法假设不同测量物理量之间是独立的，将计算测量协方差矩阵简化为计算各测量物理量的方差，然后利用经验法得到各测量物理量的方差^[6-8]。然而，传统方法存在两个问题：首先，从同一探测信号中提取的ToA和多普勒信息实际上是高度相关的；其次，经验方法无法提供测量物理量的精确方差^[9-11]。这些问题导致使用传统方法得到的协方差矩阵精度受限。为提高协方差矩阵的精度，研究者提出了诸如短时傅里叶变换(Short-time Fourier transform，STFT)^[12]、Wigner分布^[13]和分数阶傅里叶变换(Fractional Fourier transform，FrFT)^[14]等方法。虽然这些方法能够较为准确地估计ToA和多普勒因子各自的方差，但它们仍然未能有效捕捉两者之间的相关性。

本文针对ToA与多普勒因子之间高度相关导致精确协方差矩阵难以计算的问题，提出了一种基于时延-多普勒联合估计的协方差矩阵计算方法。与传统方法不同，该方法直接从信号模型出发推导协方差矩阵。考虑到水下双基地声呐系统所用探测信号通常具有大时宽带宽积，传统的窄带模型已不适用^[15-17]，本方法首先构建宽带水声信道模型以准确表征基站接收的目标回波信号；进而基于此模型推导信道参数估计的克拉美罗下界(Cramer-Rao lower bound，CRLB)；最后，利用信道参数与时延、多普勒因子之间的数学关系，推导出二者联合估计的CRLB，并将其作为协方差矩阵。仿真结果表明，所提方法有效提升了时延和多普勒协方差矩阵的精度，验证了该计算方法的有效性。

1 测量模型

本文使用双基地声呐测量时延和多普勒信息实现对水下运动目标的定位。在k时刻的测量场景如图 1所示，目标位于未知位置s_k且运动速度为v_k，发射基站在已知位置t发射探测信号并接收来自目标的回波信号，同时在已知位置u的接收基站被动接收来自目标的回波信号。

探测信号从发射基站发射到接收的双程距离R_{t, k}为

$ R_{t, k}=2 d_{t, k}=2\left\|\mathit{\boldsymbol{\boldsymbol{s}}}_k-\mathit{\boldsymbol{\boldsymbol{t}}}\right\| 。$

(1)

式中：‖·‖表示欧几里得范数; d_{t, k}=‖s_k－t‖代表发射基站到目标的距离。相应地，传播时延τ_{t, k}可表示为

$ \tau_{t, k}=R_{t, k} / c \text { 。} $

(2)

式中c为水声信号传播速度。另外，在信道的时间相干性较好，大于信号时间展宽后的处理时宽时，探测信号从发射基站经目标反射到接收基站的距离R_{u, k}为

$ R_{u, k}=d_{t, k}+d_{u, k} \text { 。} $

(3)

式中d_{u, k}=‖s_k－u‖为目标到接收基站的距离，相应的传播时延τ_{u, k}表示为

$ \tau_{u, k}=R_{u, k} / c \text { 。} $

(4)

式(2)和(4)定义了理想的传播时延。由于测量噪声的存在，时延的测量值$\tilde{\tau}_{t, k}，\tilde{\tau}_{u, k}$可表示为

$ \tilde{\tau}_{t, k}=\tau_{t, k}+n_{\tau_{t, k}}，$

(5)

和

$ \tilde{\tau}_{u, k}=\tau_{u, k}+n_{\tau_{u, k}} 。$

(6)

式中n_{τt, k}和n_{τu, k}是零均值高斯噪声且方差分别为δ_τ²_{t, k}和δ_τ²_{u, k}。

定义目标相对发射基站和接收基站接近时径向速度为正，因此目标相对发射基站和接收基站的径向速度可分别表示为

$ v_{t, k}=\frac{\left(t-s_k\right)^{\mathrm{T}} \boldsymbol{v}_k}{\left\|t-s_k\right\|}，$

(7)

和

$ v_{u, k}=\frac{\left(\boldsymbol{u}-\boldsymbol{s}_k\right)^{\mathrm{T}} \boldsymbol{v}_k}{\left\|\boldsymbol{u}-\boldsymbol{s}_k\right\|} 。$

(8)

根据以上径向速度，发射基站和接收基站的多普勒因子可表示为

$ \alpha_{t, k}=\frac{c+v_{t, k}}{c-v_{t, k}}，$

(9)

和

$ \alpha_{u, k}=\frac{c+v_{t, k}}{c-v_{u, k}} 。$

(10)

同样地，由于测量噪声的存在，发射基站和接收基站的多普勒因子的测量值可表示为

$ \tilde{\alpha}_{t, k}=\alpha_{t, k}+n_{\alpha_{t, k}}，$

(11)

和

$ \tilde{\alpha}_{u, k}=\alpha_{u, k}+n_{\alpha_{u, k}} 。$

(12)

式中n_{αt, k}和n_{αu, k}是零均值高斯噪声且方差分别为δ_α²_{t, k}和δ_α²_{u, k}。为了表示方便，将时延和多普勒因子写成向量形式m_k=[τ_{t, k}, α_{t, k}, τ_{u, k}, α_{u, k}]^T，相应的测量向量为$\tilde{\boldsymbol{m}}_k=\left[\tilde{\tau}_{t, k}, \tilde{\alpha}_{t, k}, \tilde{\tau}_{u, k}, \tilde{\alpha}_{u, k}\right]^{\mathrm{T}} 。$时延和多普勒因子的测量值和真实值间满足以下关系

$ \tilde{\boldsymbol{m}}_k=\boldsymbol{m}_k+\boldsymbol{n}_k 。$

(13)

式中n_k是协方差为Q_k的零均值高斯测量噪声向量，满足$\boldsymbol{n}_k \sim \mathcal{N}\left(\mathbf{0}, \boldsymbol{Q}_k\right) \text { 。}$

2 问题描述

双基地声呐对运动目标的定位可通过式(13)构建相关定位模型来估计目标的位置s_k和速度v_k。通常CRLB作为评估定位性能的基准，将待求的目标位置和速度写成向量形式x_k=[s_k; v_k]，基于模型(13)的目标位置和速度估计的CRLB (用$C_{\widehat{x}_k}$表示)满足

$ \mathbb{E}\left\{\left(\tilde{\boldsymbol{x}}_k-\boldsymbol{x}_k\right)\left(\tilde{\boldsymbol{x}}_k-\boldsymbol{x}_k\right)^{\mathrm{T}}\right\} \succeq \boldsymbol{C}_{\widehat{x_k}}=\boldsymbol{F}\left(\boldsymbol{x}_k\right)^{-1} 。$

(14)

式中：$\mathbb{E}\{\cdot\}$表示期望；F(x_k)∈R^4×4是费雪信息矩阵(Fisher information matrix，FIM)，可表示为

$ \boldsymbol{F}\left(\boldsymbol{x}_k\right)=\mathrm{E}\left\{\frac{\ln p\left(\tilde{\boldsymbol{m}}_k \mid \boldsymbol{x}_k\right)}{\partial \boldsymbol{x}_k} \cdot \frac{\ln p\left(\tilde{\boldsymbol{m}}_k \mid \boldsymbol{x}_k\right)}{\partial \boldsymbol{x}_k^{\mathrm{T}}}\right\} 。$

(15)

式中$\ln p\left(\tilde{m}_k \mid x_k\right)$是关于x_k的对数条件概率密度函数，可表示为

$ \ln p\left(\tilde{\boldsymbol{m}}_k \mid \boldsymbol{x}_k\right)=-\frac{1}{2}\left(\tilde{\boldsymbol{m}}_k-\boldsymbol{m}_k\right)^{\mathrm{T}} \boldsymbol{Q}_k^{-1}\left(\tilde{\boldsymbol{m}}-\boldsymbol{m}_k\right) \text { 。} $

(16)

根据式(16)可推导出

$ \frac{\partial \ln p\left(\tilde{\boldsymbol{m}}_k \mid \boldsymbol{x}_k\right)}{\partial \boldsymbol{x}_k^{\mathrm{T}}}=\left(\tilde{\boldsymbol{m}}_k-\boldsymbol{m}_k\right)^{\mathrm{T}} \boldsymbol{Q}_k^{-1} \frac{\partial \boldsymbol{m}_k}{\partial \boldsymbol{x}_k^{\mathrm{T}}} 。$

(17)

将式(17)代入式(15)，FIM可进一步表示为

$ \boldsymbol{F}\left(\boldsymbol{x}_k\right)=\frac{\partial \boldsymbol{m}_k}{\partial \boldsymbol{x}_k} \boldsymbol{Q}_k^{-1} \frac{\partial \boldsymbol{m}_k}{\partial \boldsymbol{x}_k^{\mathrm{T}}}，$

(18)

其中

$ \frac{\partial \boldsymbol{m}_k}{\partial \boldsymbol{x}_k^{\mathrm{T}}}=\left[\frac{\partial \tau_{t, k}}{\partial \boldsymbol{x}_k^{\mathrm{T}}} ; \frac{\partial \alpha_{t, k}}{\partial \boldsymbol{x}_k^{\mathrm{T}}} ; \frac{\partial \tau_{u, k}}{\partial \boldsymbol{x}_k^{\mathrm{T}}} ; \frac{\partial \alpha_{u, k}}{\partial \boldsymbol{x}_k^{\mathrm{T}}}\right] 。$

(19)

式(19)中各偏导数的推导详见附录。

显然，根据式(14)和(18)可以看到，双基地定位性能与测量协方差矩阵Q_k密切相关。由于不同基站的测量是独立的^[18]，测量协方差矩阵Q_k可表示为

$ \boldsymbol{Q}_k=\left[\begin{array}{cc} \boldsymbol{R}_t & \mathbf{0} \\ \mathbf{0} & \boldsymbol{R}_u \end{array}\right] 。$

(20)

式中R_t和R_u分别是发射基站和接收基站获得时延和多普勒因子的测量协方差矩阵，可分别表示为

$ \boldsymbol{R}_t=\left[\begin{array}{cc} \operatorname{var}\left(\tilde{\tau}_{t, k}\right) & \operatorname{cov}\left(\tilde{\tau}_{t, k}, \tilde{\alpha}_{t, k}\right) \\ \operatorname{cov}\left(\tilde{\tau}_{t, k}, \tilde{\alpha}_{t, k}\right) & \operatorname{var}\left(\tilde{\alpha}_{t, k}\right) \end{array}\right]，$

(21)

和

$ \boldsymbol{R}_u=\left[\begin{array}{cc} \operatorname{var}\left(\tilde{\tau}_{u, k}\right) & \operatorname{cov}\left(\tilde{\tau}_{u, k}, \tilde{\alpha}_{u, k}\right) \\ \operatorname{cov}\left(\tilde{\tau}_{u, k}, \tilde{\alpha}_{u, k}\right) & \operatorname{var}\left(\tilde{\alpha}_{u, k}\right) \end{array}\right] 。$

(22)

式中var(·)和cov(·, ·)分别表示求方差和协方差。

现有计算时延和多普勒因子测量协方差矩阵的方法为经验法^{[8, 10]}，该方法首先假设时延和多普勒因子测量不相关，即

$ \begin{aligned} & \operatorname{cov}\left(\tilde{\tau}_{t, k}, \tilde{\alpha}_{t, k}\right)=0，\\ & \operatorname{cov}\left(\tilde{\tau}_{u, k}, \tilde{\alpha}_{u, k}\right)=0 。\end{aligned} $

(23)

然后依赖经验方法构造两者测量方差的线性关系，即

$ \begin{aligned} & \operatorname{var}\left(\tilde{\alpha}_{t, k}\right)=\beta_t \cdot \operatorname{var}\left(\tilde{\tau}_{t, k}\right)，\\ & \operatorname{var}\left(\tilde{\alpha}_{u, k}\right)=\beta_u \cdot \operatorname{var}\left(\tilde{\tau}_{u, k}\right) 。\end{aligned} $

(24)

式中β_t和β_u是经验常数值。显然经验法忽略了测量值间的强相关性且经验值β_t和β_u通常不是精确值，这导致计算的测量协方差矩阵Q_k不精确。

3 时延和多普勒因子协方差矩阵计算

本节首先构建宽带水声信道模型得到接收信号模型，然后基于接收信号推导信道参数估计的CRLB，最后根据信道参数估计的CRLB计算时延和多普勒因子协方差矩阵。

3.1 宽带水声信道模型

探测信号通过水声信道传播，常用的宽带水声信道的冲激响应可建模为

$ h(t, \tau)=A(t) \delta(\tau-\tau(t)) \text { 。} $

(25)

式中：δ(·)表示狄拉克函数；A(t)和τ(t)分别是时变幅度和时变时延。时延的变化主要是由发射基站、接收基站和目标的相对运动引起，而路径幅度的变化是由时变的海洋环境导致的；在相对短的时间内，可以假设路径的幅度A(t)为常数，即式(25)中的A(t)变成与时间无关的参数A。时延τ(t)可以近似表示为

$ \tau(t)=a-\mu t \text { 。} $

(26)

式中a是τ(t)的零阶泰勒展开和μ描述了传播路径时延变化快慢。将式(26)代入式(25)，可以重写信道冲激响应为

$ h(t, \tau)=A \delta(\tau-(a-\mu t)) \text {。} $

(27)

当发射的探测信号为s(t)时，不考虑噪声的接收信号r(t)为

$ \begin{gathered} r(t)=\int_{-\infty}^{\infty} h(t, \tau) s(t-\tau) \mathrm{d} \tau= \\ A \int_{-\infty}^{\infty} \delta(\tau-(a-\mu t)) s(t-\tau) \mathrm{d} \tau= \\ A s((1+\mu) t-a) 。\end{gathered} $

(28)

定义α=1+μ和τ′=a/(1+μ)，接收信号可表示为

$ r(t)=A s\left(\alpha\left(t-\tau^{\prime}\right)\right) \text { 。} $

(29)

式中α和τ′分别是多普勒因子和时延。

3.2 信道参数估计CRLB

由于线性调频(Linear frequency modulation，LFM)信号在时域和频域分辨率之间实现了很好的平衡，因此在水下声呐探测系统中常被使用，并可以表示为^[19]

$ s(t)=P \mathrm{e}^{\mathrm{j}\left(2 \pi f_0 t+\gamma \pi t^2+\varphi_0\right)}, t \in[0, T] 。$

(30)

式中P、f₀、γ和φ₀分别是信号幅度、起始频率、调频斜率和初始相位。将式(30)代入式(29)得到接收信号r(t)为

$ r(t)=\underbrace{\widetilde{A} e^{j\left(2 \pi \widetilde{J}_t t-\widetilde{R} \pi t^2+\widetilde{\varphi}\right)}}_{g(t)}+\omega(t) \text { 。} $

(31)

式中ω(t)是加性高斯噪声，且

$ \left\{\begin{array}{l} \tilde{f}=\alpha f_0-\gamma \alpha^2 \tau^{\prime} \\ \widetilde{k}=\gamma \alpha^2 \\ \widetilde{\varphi}=\gamma \pi \alpha^2 \tau^{\prime 2}-2 \pi f_0 \alpha \tau^{\prime}+\varphi_0 \\ \widetilde{A}=A P \end{array} 。\right. $

(32)

由式(32)可知接收信号r(t)包含了四个未知参数，汇集所有参数得到向量形式$\boldsymbol{\theta}=[\tilde{f}, \tilde{k}, \tilde{\varphi}, \tilde{A}]^{\mathrm{T}} 。$假设采样频率为f_s，即采样周期为T_s=1/f_s，信号采样序列为

$ r[n]=g[n]+w[n], n=0, 1, \cdots, N-1 。$

(33)

式中：N是总采样数；w[n]是高斯白噪声。将式(33)写成向量形式

$ r=g+w，$

(34)

其中

$ \boldsymbol{g}=(g[0], g[1], \cdots, g[N-1])^{\mathrm{T}} \text { 。} $

(35)

式中：r=r[0], r[1], …, r[N－1]^T是接收信号的采样向量；w=w[0], w[1], …, w[N－1]^T是噪声向量，且满足$w \sim \mathscr{C} \mathcal{N}\left(\mathbf{0}, \sigma^2 \mathbf{I}\right) \text { 。}$根据式(34)可知关于θ的似然函数为

$ \lambda_\theta=\frac{1}{\left(\pi \sigma^2\right)^N} \exp \left\{-\frac{(\boldsymbol{r}-\boldsymbol{g})^{\mathrm{H}}(\boldsymbol{r}-\boldsymbol{g})}{\sigma^2}\right\} 。$

(36)

随后，相应的对数似然函数可表示为

$ \psi_\theta=\ln \left(\lambda_\theta\right)=-N \ln \left(\pi \sigma^2\right)-\frac{(\boldsymbol{r}-g)^{\mathrm{H}}(\boldsymbol{r}-g)}{\sigma^2} 。$

(37)

然后，关于θ的FIM可表示为

$ \boldsymbol{F}_\theta=\mathbb{E}\left\{\boldsymbol{\nabla}_\theta \psi_\theta \boldsymbol{\nabla}_\theta^{\mathrm{T}} \psi_\theta\right\}=\frac{2}{\sigma^2} \sum\limits_{n=0}^{N-1} \Re\left\{\mathbf{∇}_\theta g[n] \mathbf{∇}_\theta^{\mathrm{H}} g[n]\right\}，$

(38)

其中

$ \boldsymbol{\nabla}_\theta g[n]=g[n] \cdot\left[\begin{array}{llll} \mathrm{j} 2 \pi n T_s & \mathrm{j} \pi n^2 T_s^2 & \mathrm{j} & 1 / \tilde{A} \end{array}\right]^{\mathrm{T}} \text { 。} $

(39)

将式(39)代入式(38)得到

$ \boldsymbol{F}_\theta=\frac{2}{\sigma^2} \sum\limits_{n=0}^{N-1}\left\{\Re\left\{g[n] g^*[n]\right\} \cdot \boldsymbol{V}\right\}，$

(40)

其中

$ \boldsymbol{V}=\left[\begin{array}{cccc} 4 \pi^2 n^2 T_s^2 & 2 \pi^2 n^3 T_s^3 & 2 \pi n T_s & 0 \\ 2 \pi^2 n^3 T_s^3 & \pi^2 n^4 T_s^4 & \pi n^2 T_s^2 & 0 \\ 2 \pi n T_s & \pi n^2 T_s^2 & 1 & 0 \\ 0 & 0 & 0 & 1 / \widetilde{A}^2 \end{array}\right] 。$

由于采样总数N是足够大的整数，因此有如下的近似关系

$ \begin{aligned} & \sum\limits_{n=0}^{N-1} n \sim N^2 / 2, \sum\limits_{n=0}^{N-1} n^2 \sim N^3 / 3，\\ & \sum\limits_{n=0}^{N=1} n^4 \sim N^4 / 2, \sum\limits_{n=0}^{N-1} n^4 \sim N^5 / 5 。\end{aligned} $

(41)

将式(41)代入式(40)且采样时间满足T=NT_s，信道参数向量θ的FIM可近似表示为

$ \begin{aligned} &\boldsymbol{F}_{\boldsymbol{0}} \sim \frac{2 \bar{A}^2 f_A}{\sigma^2} U=\\ &\frac{2 \widetilde{A}^2 f_0}{\sigma^2}\left[\begin{array}{cccc} 4 \pi^2 T^3 / 3 & \pi^2 T^4 / 2 & \pi T^2 & 0 \\ \pi^2 T^4 / 2 & \pi^2 T^5 / 5 & \pi T^3 / 3 & 0 \\ \pi T^2 & \pi T^3 / 3 & 1 & 0 \\ 0 & 0 & 1 & 1 / \tilde{A}^2 \end{array}\right] \end{aligned} $

因此可得到信道参数估计的CRLB为

$ R_0=F_0^{-1} 。$

(42)

3.3 时延-多普勒因子协方差估计

我们仅对目标位置和速度相关的参数感兴趣，即$\tilde{f}$和$\tilde{k}，$下面将推导参数向量$\boldsymbol{o}=[\widetilde{f}, \widetilde{k}]^{\mathrm{T}}$的CRLB，用符号R_o表示。为了表示方便，将矩阵U分块表示为

$ \boldsymbol{U}=\left[\begin{array}{cc} \boldsymbol{G} & \boldsymbol{B} \\ \boldsymbol{B}^{\mathrm{T}} & \boldsymbol{C} \end{array}\right]，$

(43)

其中

$ \begin{gathered} \boldsymbol{G}=\left[\begin{array}{cc} 4 \pi^2 T^3 / 3 & \pi^2 T^4 / 2 \\ \pi^2 T^4 / 2 & \pi^2 T^5 / 5 \end{array}\right]，\\ \boldsymbol{B}=\left[\begin{array}{cc} \pi T^2 & 0 \\ \pi T^3 / 3 & 0 \end{array}\right]，\\ \boldsymbol{C}=\left[\begin{array}{cc} 1 & 0 \\ 0 & 1 / \widetilde{A}^2 \end{array}\right] 。\end{gathered} $

使用分块矩阵求逆定理^[20]可得R_o为

$ \begin{gathered} \boldsymbol{R}_o=\boldsymbol{F}_{\boldsymbol{\theta}}^{-1}[1: 2, 1: 2]= \\ \frac{\sigma^2}{2 N \widetilde{A}^2}\left(\boldsymbol{G}-\boldsymbol{B C}^{-1} \boldsymbol{B}^{\mathrm{T}}\right)-1= \\ \frac{\sigma^2}{2 N \widetilde{A}^2} \cdot \frac{6}{\pi^2 T^2}\left[\begin{array}{cc} 8 & -15 / T \\ -15 / T & 30 / T^2 \end{array}\right] 。\end{gathered} $

(44)

其中F_θ[1∶2, 1∶2]表示取矩阵F_θ中第1—2行及第1—2列的元素。在获得CRLB(即R_o)后，可以根据链式法则计算需要的时延τ′和多普勒因子α的CRLB。为了表示方便，我们定义向量参数η=[τ′, α]^T。关于向量参数η的CRLB为

$ \boldsymbol{R}_\eta=\nabla_\eta^{\mathrm{T}} \boldsymbol{o} \boldsymbol{R}_o \nabla_\eta \boldsymbol{o}，$

(45)

其中

$ \boldsymbol{\nabla}_{\boldsymbol{\eta}} \boldsymbol{o}=\left[\begin{array}{cc} -\gamma_\alpha^2 & f_0-2 \gamma_{\alpha \tau} \\ 0 & 2 \gamma_\alpha \end{array}\right] 。$

我们知道有效估计器得到的估计值精度能达到CRLB。因此，可以合理地使用R_η得到发射基站和接收基站测量时延和多普勒因子的协方差矩阵R_t和R_u。

式(45)验证了同一基站得到的时延和多普勒因子之间具有强相关性，并给出了时延和多普勒因子的估计方差和协方差。因此，使用式(23)和(24)所示的经验法不能得到精确的时延和多普勒因子估计协方差矩阵。

4 仿真结果及性能分析

本节对目标不同运动的情形进行数值仿真实验，通过对比FrFT方法、短时傅里叶变换方法和Wigner分布方法的结果，验证提出的时延和多普勒因子协方差估计方法的有效性。

4.1 发射波形和接收波形生成

LFM探测信号的参数设置如下：带宽为B=5 000 Hz，脉宽为T=30 ms，发射周期为1 s，起始频率为f₀=5 000 Hz，采样频率为f_s=8B Hz。根据以上探测信号参数设置，图 2展示了200 ms时长的发射信号时域波形。

首先考虑如下的运动目标定位场景，发射基站和接收基站的位置分别为t=[100, 100]^Tm和u=[3 000, 400]^Tm。图 3和4展示了时刻k目标位置和速度分别为s_k=[1 500, 1 500]^Tm和v_k=[3,5]^Tm/s时，进行1次仿真得到的发射基站和接收基站在不同信噪比(Signal-to-noise ratio，SNR)下基于发射波形接收的回波时域波形；图 5和6展示了目标速度为v_k=[－5, －3]^Tm/s时，进行1次仿真得到的发射基站和接收基站在不同信噪比下接收的回波时域波形。

4.2 时延和多普勒协方差矩阵估计性能

基于5 000次蒙特卡洛仿真得到发射基站和接收基站接收的回波，使用本文提出的方法、FrFT方法、短时傅里叶变换和Wigner分布方法分别得到时延和多普勒因子估计方差。图 7和8显示了在目标速度为v_k=[3,5]^Tm/s，不同信噪比下使用本文提出的方法、FrFT方法、短时傅里叶变换和Wigner分布方法得到的时延和多普勒因子估计方差。从图 7、8中可以看到，相比其他方法，本文提出的方法估计的时延方差和多普勒因子方差更小，这说明所提方法的准确性更高。

为了进一步验证所提方法的普适性，图 9和图 10显示了在目标速度为v_k=[－5, －3]^Tm/s时，不同信噪比下使用本文提出的方法、FrFT方法、短时傅里叶变换和Wigner分布方法得到的时延和多普勒因子估计方差，从图 9、10中可看出结果与目标速度v_k=[3,5]^Tm/s时一致。

对于时延和多普勒因子估计的协方差矩阵，除了计算时延和多普勒因子的估计方差，还需要计算时延和多普勒因子之间的协方差。图 11和图 12分别显示了目标速度为v_k=[3,5]^Tm/s和v_k=[－5, －3]^Tm/s时，不同信噪比下使用本文提出的方法和传统方法使用经验法得到的时延和多普勒因子估计的协方差，可以看到时延和多普勒因子间是强相关的，而经验法假设时延和多普勒因子是独立的，即

$ \operatorname{cov}\left(\tilde{\tau}_{t, k}, \tilde{\alpha}_{t, k}\right)=0，$

和

$ \operatorname{cov}\left(\tilde{\tau}_{u, k}, \tilde{\alpha}_{u, k}\right)=0 。$

这种忽略时延和多普勒因子间相关性的经验法将导致得到的时延和多普勒因子估计协方差矩阵不精确。

4.3 双基地定位性能

为了进一步验证所提方法的有效性，接下来从双基地定位性能的角度进行评估。由于发射基站和接收基站在接收回波信号时的噪声水平不同，因此设置发射基站和接收基站的信噪比从-25 dB到20 dB变化。上节结果显示FrFT方法相比另外两种传统方法(STFT和Wigner分布方法)具有最小的时延和多普勒因子的估计方差，在不丢失一般性的情况下，为了更清晰地在图中展现所提方法与传统方法的定位性能，下面将基于本文提出的方法与FrFT方法使用经验法得到的时延和多普勒因子的协方差矩阵比较双基地定位精度。

图 13和图 14展示了目标速度为v_k=[3,5]^Tm/s时不同信噪比下使用本文提出的方法和FrFT方法估计双基地目标位置和速度的精度。可以看到，相比忽略时延和多普勒因子间强相关性的传统方法，本文提出的方法具有更高的目标位置和速度估计精度。

此外，对于不同的目标速度(v_k=[－5, －3]^Tm/s)，使用本文提出的方法和FrFT方法估计双基地目标位置和速度的精度。从图 15和图 16中可以看到本文提出的方法具有更高的目标位置和速度估计精度，但相较于目标速度为v_k=[3,5]^Tm/s时的结果目标位置和速度估计误差更大，这是由于目标速度的变化造成时延和多普勒因子的估计方差和协方差变大导致的(见图 7—12)。

综合以上仿真结果可知, 本文提出的方法能计算出精确的时延和多普勒因子的协方差矩阵。

5 结语

针对时延和多普勒之间的强相关性导致计算精确的协方差矩阵难的问题，本文研究了基于时延和多普勒联合估计的协方差计算方法。该方法通过构建宽带水声信道模型得到了接收信号，并基于得到的接收信号推导了估计信道参数的理论CRLB；考虑到时延-多普勒因子与信道参数是非线性相关的，因此基于链式法则推导了时延和多普勒因子的CRLB；结合有效估计器的估计误差接近CRLB的特性，可以根据时延和多普勒因子的CRLB得到其协方差矩阵。仿真实验表明：与FrFT方法、短时傅里叶变换和Wigner分布方法相比，本文所提方法估计的时延和多普勒因子的方差与之接近，说明所提方法的有效性。与经验法相比，本文提出的方法考虑了时延和多普勒因子间强相关性，得到了精度更高的协方差矩阵；另外提出的方法也具有更高的目标位置和速度估计精度。

本论文从理论上证实了所提方法的有效性，但缺乏与真实试验数据的对比结果，因此在后续工作中将开展双基地声呐系统试验，根据获取的真实试验数据进一步验证所提方法的有效性。

附录

式(19)中各偏导数如下：

$ \begin{gathered} \frac{\partial \tau_{t, k}}{\partial x_k^{\mathrm{T}}}=\left[2 \boldsymbol{\rho}_{t-s_k} ; \mathbf{0}_{2 \times 1}\right]^{\mathrm{T}}，\\ \frac{\partial \tau_{u, k}}{\partial x_k^{\mathrm{T}}}=\left[2 \boldsymbol{\rho}_{t-s_k} ; \mathbf{0}_{2 \times 1}\right]^{\mathrm{T}}，\\ \frac{\partial \alpha_{t, k}}{\partial x_k^{\mathrm{T}}}=\left[\frac{2 c \mathcal{P}_{t-s_k} v_k}{\left(c-v_{t, k}\right)^2\left\|t-s_k\right\|} ; \frac{2 c \rho_{t-s_k}}{\left(c-v_{t, k}\right)^2}\right]^{\mathrm{T}}，\\ \frac{\partial \alpha_{n, k}}{\partial x_k^{\mathrm{T}}}=\left[\begin{array}{c} \frac{\mathcal{P} \frac{1}{t-s_k} v_k}{\left(c-v_{u, k}\right)\left\|t-s_k\right\|}+\frac{\left(c+v_{t, k}\right) \mathcal{P}_{n-s_k}^{\perp} v_k}{\left(c-v_{n, k}\right)^2\left\|u-s_k\right\|} \\ \frac{t-s_k}{\left(c-v_{u, k}\right)\left\|t-s_k\right\|}+\frac{\left.\left(c+v_{n, k}\right) u-s_k\right)}{\left(c-v_{n, k}\right)^2\left\|u-s_k\right\|} \end{array}\right]^{\mathrm{T}} 。\end{gathered} $

式中：$\boldsymbol{\rho}_a=\frac{a}{\|a\|} ；$$\mathcal{P}_{\boldsymbol{a}}^{\perp}=\mathbf{I}_{2 \times 2}-\frac{\boldsymbol{a} \boldsymbol{a}^{\mathrm{T}}}{\|\boldsymbol{a}\|^2} ；$a表示t－s_k或u－s_k。