11、多目标模型中静态参数的近似似然估计

关注

原创最新推荐文章于 2026-05-23 14:02:01 发布 57 阅读 AI 写同款· GEO 优化 ›

多目标模型中静态参数的近似似然估计

1. 引言

目标跟踪问题旨在在线估计被监视目标随时间变化的状态向量。在单目标跟踪中，常用隐马尔可夫模型（HMM）描述目标状态和观测值的关系：
- (X_0 \sim \mu_{\theta}(\cdot))
- (X_n|(X_{n - 1} = x_{n - 1}) \sim f_{\theta}(\cdot|x_{n - 1}), n \geq 1)
- (Y_n|X_n = x_n \sim g_{\theta}(\cdot|x_n), n \geq 1)

其中，(\theta) 表示模型参数向量。而多目标跟踪问题更为复杂，每个时刻的状态 (X_n) 和观测 (Y_n) 是随机有限集：
- (X_n = {X_{n,1}, X_{n,2}, \ldots, X_{n,K_n}})
- (Y_n = {Y_{n,1}, Y_{n,2}, \ldots, Y_{n,M_n}}, n \geq 1)

目标数量 (K_n) 会随时间变化，还存在未被检测到的目标和虚假测量值。多目标跟踪的目标是根据截至当前时刻的所有观测集估计时变状态集，同时校准多目标跟踪模型的参数也是重要问题。

为了从数据中估计模型参数，本文采用近似似然函数并将其最大化的方法。当将未知数量的目标及其初始状态视为泊松点过程（PP）的实现时，(X_n) 的分布也为泊松分布。若能将 (X_n|Y_{1:n}) 近似为泊松 PP，就能大大简化估计目标数量和状态的问题。概率假设密度（PHD）滤波器可传播条件密度序列 ({p(x_n|y_{1:n})}_{n \geq 0}) 的强度，在多目标跟踪中展现出潜力。

2. 多目标模型

多目标统计模型中，泊松 PP 起着关键作用。未观测点 (X_n) 按以下方式演化到 (X_{n + 1})：
- 每个 (X_n) 中的点以概率 (p_{\theta}^S(x)) 存活，并根据马尔可夫转移密度 (f_{\theta}(x_{n + 1}|x_n)) 位移。
- 新点从强度函数为 (\gamma_{\theta}(x)) 的泊松过程中“诞生”。

(X_{n + 1}) 中的点通过以下模型被观测：
- 每个点以概率 (p_{\theta}^D(x)) 生成观测值，观测值通过密度 (g_{\theta}(y|x_{n + 1,j})) 产生。
- 除了这些检测到的点，还存在来自强度函数为 (\kappa_{\theta}(y)) 的独立泊松过程的虚假测量值。

3. PHD 滤波器回顾

3.1 部分观测泊松过程的推断

设 (Y) 的实现为 (y = {y_1, \ldots, y_m})，一般难以精确刻画 (X|Y = y) 的分布 (P_{X|y})。在泊松先验的情况下，文献 [21] 利用泊松 PP 的相关性质，给出了 (P_{X|y}) 的闭式表达式。

将 (X) 分解为 (bX)（产生观测值的点）和 (eX)（未观测到的点），(bX) 和 (eX) 均为独立泊松过程，强度分别为 (\alpha(x)p_D(x)) 和 (\alpha(x)(1 - p_D(x)))。

构造标记点过程 (Z)，其强度为 (\alpha(x)p_D(x)g(y|x)I_{E_2\times E_1}(y, x)dxdy + \kappa(y)I_{E_2\times{\Delta}}(y, x)\delta_{\Delta}(dx)dy)。根据定理 11.1，(P_{X|y}) 是强度为 (\alpha(x)(1 - p_D(x))dx) 的泊松 PP 和限制在 (E_1) 上的 (m) 点 PP 的叠加。

推论 11.2 给出了对于有界实值可测函数 (\phi) 在 (E_1) 上的期望表达式，同时得到 (P_{X|y}) 的强度为 (\alpha(x)(1 - p_D(x)) + \sum_{i = 1}^{m} \frac{\alpha(x)p_D(x)g(y_i|x)}{\int_{E_1} \alpha(x)p_D(x)g(y_i|x)dx + \kappa(y_i)})。

3.2 PHD 滤波器

PHD 滤波器基于以下两个事实：
- 若 (X_{n - 1}) 是强度为 (\alpha_{n - 1}(x_{n - 1})) 的泊松 PP，则 (X_n) 也是泊松 PP，强度为 (\alpha_n(x_n) = \int_{E_1} f_{\theta}(x_n|x_{n - 1})p_{\theta}^S(x_{n - 1})\alpha_{n - 1}(x_{n - 1})dx_{n - 1} + \gamma_{\theta}(x_n))。
- (X_1|y_1) 的最佳泊松近似是与 (P_{X_1|y}) 具有相同强度的泊松 PP，强度为 (\alpha_{1|1}(x_1) = \left[1 - p_{\theta}^D(x_1) + \sum_{y \in y_1} \frac{p_{\theta}^D(x_1)g_{\theta}(y|x_1)}{\int_{E_1} p_{\theta}^D(x)g_{\theta}(y|x)\alpha_{1|0}(x)dx + \kappa_{\theta}(y)}\right]\alpha_{1|0}(x_1))。

递归应用这两个事实得到 PHD 滤波器：
- (\alpha_{n|n - 1} = (\Phi\alpha_{n - 1|n - 1}) + \gamma_{\theta})
- (\alpha_{n|n} = (\Psi_n\alpha_{n|n - 1}))

4. 近似边际似然

根据模型，边际似然 (p(y_{1:n})) 可分解为 (p(y_{1:n}) = p(y_1) \prod_{k = 2}^{n} p(y_k|y_{1:k - 1}) = \int p(y_1|x_1)p(x_1)dx_1 \prod_{k = 2}^{n} \int p(y_k|x_k)p(x_k|y_{1:k - 1})dx_k)。

利用 PHD 滤波器对 (p(x_k|y_{1:k - 1})) 的泊松近似，可得预测似然 (p(y_k|y_{1:k - 1})) 也是泊松分布，强度函数为 (\int_{E_1} g_{\theta}(y|x_n)p_{\theta}^D(x_n)\alpha_{n|n - 1}(x_n)dx_n + \kappa_{\theta}(y))。

结合上述结果，观测数据的对数似然可近似为：
(\ell_{Po,n}(\theta) = - \sum_{k = 1}^{n} \int \left[\int g_{\theta}(y_k|x_k)p_{\theta}^D(x_k)\alpha_{k|k - 1}(x_k)dx_k + \kappa_{\theta}(y_k)\right] dy_k + \sum_{k = 1}^{n} \sum_{r = 1}^{m_k} \log \left[\int g_{\theta}(y_{k,r}|x_k)p_{\theta}^D(x_k)\alpha_{k|k - 1}(x_k)dx_k + \kappa_{\theta}(y_{k,r})\right])

5. PHD 滤波器及其梯度的 SMC 近似

多数情况下，无法精确计算近似对数似然 (\ell_{Po,n}(\theta))，需要数值方法近似其中的积分。文献 [27] 提出的顺序蒙特卡罗（SMC）方法可用于近似强度函数，其步骤如下：
- 预测步骤 ：
- 给定上一次迭代的粒子集 ({X_{n - 1}^{(i)}, W_{n - 1}^{(i)}} {i = 1}^{N})，从提议分布 (q_n(\cdot|X {n - 1}^{(i)})) 中抽取 (N) 个样本 ({X_n^{(i)}} {i = 1}^{N})，计算预测重要性权重 (W {n|n - 1}^{(i)} = \frac{f_{\theta}(X_n^{(i)}|X_{n - 1}^{(i)})p_{\theta}^S(X_{n - 1}^{(i)})}{q_n(X_n^{(i)}|X_{n - 1}^{(i)})}W_{n - 1}^{(i)})。
- 抽取 (L) 个额外样本 ({X_n^{(i)}} {i = N + 1}^{N + L}) 用于出生项，计算相应的重要性权重 (W {n|n - 1}^{(i)} = \frac{1}{L} \frac{\gamma_{\theta}(X_n^{(i)})}{p_n(X_n^{(i)})})。
- 近似预测强度 (\hat{\alpha} {n|n - 1}(dx_n) = \sum {i = 1}^{N + L} W_{n|n - 1}^{(i)}\delta_{X_n^{(i)}}(dx_n))。
- 更新步骤 ：
- 重新加权粒子，得到新的重要性权重 (W_n^{(i)} = \left[1 - p_{\theta}^D(X_n^{(i)}) + \sum_{r = 1}^{m_n} \frac{\psi_{n,r}(X_n^{(i)})}{\hat{Z} {n,r}}\right]W {n|n - 1}^{(i)})，其中 (\psi_{n,r}(x) = p_{\theta}^D(x)g_{\theta}(y_{n,r}|x))，(\hat{Z} {n,r} = \sum {i = 1}^{N + L} \psi_{n,r}(X_n^{(i)})W_{n|n - 1}^{(i)} + \kappa_{\theta}(y_{n,r}))。
- 近似更新强度 (\hat{\alpha} {n|n}(dx_n) = \sum {i = 1}^{N + L} W_n^{(i)}\delta_{X_n^{(i)}}(dx_n))。
- 归一化重要性权重并重新采样。

对于 PHD 滤波器的梯度，采用文献 [17, 18] 提出的方法进行近似。算法 11.1 总结了计算 (\hat{\ell} {Po,n}) 和 (\hat{\nabla}\ell {Po,n}) 的 SMC 方法，其计算成本与粒子数量 (N) 的平方成正比。

以下是 PHD 滤波器 SMC 近似的流程图：

graph TD;
    A[初始化粒子集] --> B[预测步骤];
    B --> C[计算预测强度和权重];
    C --> D[更新步骤];
    D --> E[重新加权粒子];
    E --> F[计算更新强度];
    F --> G[归一化和重新采样];
    G --> H[计算似然和似然梯度];
    H --> I[结束];

6. 参数估计

6.1 逐点梯度近似

使用梯度上升算法估计模型参数，可采用离线或在线方式。离线方式下，对于离散时间索引 (k)，通过 (\theta_{k + 1} = \theta_k + a_{k + 1} \nabla\ell_{Po,T}(\theta)| {\theta = \theta_k}) 更新参数，其中 ({a_k} {k \geq 1}) 是步长序列，需满足 (\sum_{k} a_k = \infty) 和 (\sum_{k} a_k^2 < \infty)。

在线方式下，可采用递归更新方案 (\theta_{n + 1} = \theta_n + a_{n + 1} \nabla\log p_{Po}(y_n|y_{1:n - 1})|_{\theta = \theta_n})。在实践中，可先采用恒定较小的步长，避免算法过早陷入局部最优。

6.2 同时扰动随机近似（SPSA）

也可使用 SMC 进行无显式梯度计算的似然最大化。SPSA 方法通过有限差分（FD）近似梯度，对于高维问题更为高效。在第 (k) 次迭代中，梯度估计为 (\hat{\nabla} p\ell {Po,T}(\theta_k) = \frac{\hat{\ell} {Po,T}(\theta_k + c_k\Delta_k) - \hat{\ell} {Po,T}(\theta_k - c_k\Delta_k)}{2c_k\Delta_{k,p}})，其中 (\Delta_k) 是随机扰动向量，({c_k}_{k \geq 0}) 是递减的小正数序列。

SPSA 方法的步骤总结在算法 11.2 中，其计算成本与粒子数量线性相关，但需要选择合适的序列 ((a_n)) 和 ((c_n))。

以下是参数估计方法的对比表格：
| 方法 | 优点 | 缺点 | 计算成本 |
| ---- | ---- | ---- | ---- |
| 逐点梯度近似 | 估计梯度方差低 | 计算梯度可能昂贵，需要调整步长 | 与粒子数量平方成正比 |
| SPSA | 实现简单，适用于高维问题 | 需要选择序列，可能需要手动调参 | 与粒子数量线性相关 |

7. 模拟研究

7.1 模型

采用多目标模型进行参数估计方法的评估，假设单个目标具有恒定速度模型，位置限制在 ([0, 100]×[0, 100]) 窗口内。目标状态由四维向量表示，状态演化方程为 (X_n = \begin{bmatrix}1 & 1 & 0 & 0 \ 0 & 1 & 0 & 0 \ 0 & 0 & 1 & 1 \ 0 & 0 & 0 & 1\end{bmatrix} X_{n - 1} + \begin{bmatrix}V_n(1) \ V_n(2) \ V_n(3) \ V_n(4)\end{bmatrix})，其中 (V_n(1)) 和 (V_n(3)) 是均值为 0、标准差为 (\sigma_{xs} = 0.01) 的独立高斯随机变量，(V_n(2)) 和 (V_n(4)) 是均值为 0、标准差为 (\sigma_{xv} = 0.25) 的独立高斯随机变量。

出生强度为 (\gamma_{\theta}(\cdot) = \Gamma N(\cdot; \mu_b, \Sigma_b))，观测值在加性各向同性高斯噪声中，噪声标准差为 (\sigma_y)，杂波强度 (\kappa(y)) 在 ([0, 100]×[0, 100]) 上均匀分布。待推断的模型参数为 (\theta = [\sigma_y, \kappa, \Gamma, \mu_{bx}, \mu_{by}, \sigma_{bs}, \sigma_{bv}, p_D]^T)，((\sigma_{xs}, \sigma_{xv}, p_S)) 假设已知。

7.2 逐点梯度近似

使用参数 (\theta^* = [5, 4×10^{-4}, 1, 50, 50, 5, 2, 0.9]^T) 生成观测记录，观测时间长度为 50 时，平均每个时间点有 9 个观测值，其中 4 个为虚假观测值，总共 50 个目标。使用 (N = 400) 和 (L = 400) 个粒子。

实验结果表明，估计的模型参数收敛到真实参数的附近，但存在一些差异。重复实验 50 次，不同观测记录长度下的结果显示，较长的观测记录可提高估计的准确性，但部分参数存在小偏差，可能是由于粒子数量不足或近似似然本身的问题。

7.3 SPSA

使用与逐点梯度方法相同的数据记录和初始条件运行 SPSA 方案，选择恒定步长 (a_k = a = 1×10^{-5}) 和 (c_k = c = 0.025)。经过 20000 次迭代后，得到的参数值与逐点梯度方法的结果相近，小差异归因于有限差分梯度近似的偏差。SPSA 方案实现更简单，但需要选择序列 ((a_n)) 和 ((c_n))，可能需要更多手动调参。

综上所述，通过近似边际似然和 SMC 方法，可以有效地估计多目标模型的静态参数。虽然结果令人鼓舞，但仍需进一步研究近似似然引入的偏差问题。

多目标模型中静态参数的近似似然估计

8. 实验结果分析

在模拟研究中，逐点梯度近似和 SPSA 两种方法都对多目标模型的参数估计进行了尝试。从逐点梯度近似的实验结果来看，估计参数能够收敛到真实参数附近，但存在一定偏差。例如，在估计杂波强度时出现了高估的情况，这可能是由于模拟中使用的粒子数量不足，无法充分跟踪强度函数的所有模式，导致目标产生的观测值被错误地归因于更高的杂波强度。

而 SPSA 方法在实现上更为简单，计算成本与粒子数量线性相关。在相同的数据记录和初始条件下，经过一定次数的迭代，得到的参数值与逐点梯度近似方法相近。不过，SPSA 方法需要选择合适的序列 ((a_n)) 和 ((c_n))，这可能需要更多的手动调参工作。

为了更直观地比较两种方法的性能，我们可以列出以下对比：
| 方法 | 收敛速度 | 参数估计准确性 | 实现复杂度 | 手动调参需求 |
| ---- | ---- | ---- | ---- | ---- |
| 逐点梯度近似 | 相对较慢，受粒子数量和步长影响 | 存在一定偏差，受粒子数量和近似似然影响 | 较高，需要计算梯度 | 较小，但步长选择需谨慎 |
| SPSA | 相对较快，受序列选择影响 | 与逐点梯度近似相近，有有限差分偏差 | 较低，无需显式计算梯度 | 较大，需要选择合适序列 |

9. 实际应用考虑

在实际应用中，选择合适的参数估计方法需要综合考虑多个因素。如果对计算效率要求较高，且数据维度较高，SPSA 方法可能是一个更好的选择。因为它无需显式计算梯度，实现简单，计算成本与粒子数量线性相关。但需要注意的是，SPSA 方法的性能很大程度上依赖于序列 ((a_n)) 和 ((c_n)) 的选择，需要进行一定的手动调参。

如果对参数估计的准确性要求较高，且计算资源充足，逐点梯度近似方法可能更合适。虽然它的计算成本与粒子数量的平方成正比，但通过合理选择粒子数量和步长，可以获得较为准确的参数估计。不过，该方法的实现复杂度较高，需要计算梯度，并且步长的选择需要谨慎，以避免算法陷入局部最优。

以下是实际应用中选择方法的决策流程图：

graph TD;
    A[开始] --> B{计算效率要求高?};
    B -- 是 --> C{数据维度高?};
    C -- 是 --> D[SPSA方法];
    C -- 否 --> E{参数估计准确性要求高?};
    E -- 是 --> F[逐点梯度近似方法];
    E -- 否 --> D;
    B -- 否 --> E;