旋转长基线干涉仪“维特比-卡尔曼”滤波定位算法

2021-08-12刘元华

航天电子对抗 2021年3期

马爽，蔡凯，刘元华

（北京遥感信息研究所，北京 100011）

0 引言

旋转长基线干涉仪（RLBI）通过不断旋转形成时序上的多条虚拟基线来完成相位差解模糊和定位。与传统的多通道干涉仪和阵列测向等方法相比，旋转基线体制具有设备简单、通道一致性要求低等优点，但是，长基线干涉仪只能测量得到模糊相位差，如何利用模糊相位差高效准确地定位，是旋转干涉仪定位体制需要解决的首要问题。

要解决该问题，最直接的思路是首先将模糊相位差恢复成真实相位差，再进行定位解算。文献[1-4]分别研究了圆阵、GPS测向、导弹导引头测向和近场测向等条件下相位差解模糊的方法，实现了针对具体应用场景的定位解算。但是，该类方法对信号测量间隔和精度要求高，一般需要满足相对无模糊条件才能实现。为了扩展适用范围，文献[5-7]探索了利用模糊相位差直接定位的方法，在定位解空间利用各种优化算法直接搜索最优解。该类方法不再需要恢复真实相位差，但是，计算量随着可能定位区域的扩大迅速增加，且定位精度与搜索粒度密切相关，在具体应用中受到一定限制。

本文在研究旋转长基线干涉仪相位差和模糊数变化特点的基础上，提出了相位差“条带”和“路径”的概念，将需要搜索的解空间限制在模糊数序列空间中，与直接在定位空间中搜索相比，大大降低了解空间的维度。用隐马尔可夫模型（HMM）对“条带”转移过程和相位差观测过程进行建模，通过结合“卡尔曼”滤波过程和“维特比”状态解码过程，同步实现了相位差解模糊和滤波定位。该算法能够对数据进行序贯处理，实时完成解模糊和定位，

定位性能与使用真实相位差进行定位时的“卡尔曼”滤波算法一致。

1 旋转长基线干涉仪相位差观测过程HMM模型

旋转基线通过不断旋转来获取目标在不同基线转角下的相位差测量值。在旋转过程中，目标与基线的相对位置不断变化，对应的模糊数和相位差测量值也随之改变，这一过程可以用HMM[8]来建模。

1.1 模糊相位差、“条带”和“路径”

旋转基线干涉仪采用单根基线测量信号相位差，真实相位差为：

式中，d为基线长度，λ是信号波长，θ为信号到达方向与基线法平面的夹角。

干涉仪输出的测量相位差φ′∈(-π,π]，当d＞λ时，存在相位差模糊，测量相位差φ′与真实相位差φ之间满足如下关系：

式中，mod2π(·)表示按照2π取模，整数n称为模糊数。

模糊数n的取值范围为[-N,N]，N为不超过dλ的整数。例如，当信号频率为3 GHz、基线长度为1 m时，模糊数的取值范围为[-10,10]。

假设干涉仪基线在定位区域上方水平放置，定位区域内每个点对应的真实相位差φ、模糊相位差φ′和模糊数n具有如下特点：

1）真实相位差φ连续分布，并沿基线矢量方向单调变化。φ的取值范围由基线波长比d/λ决定，d/λ越大，取值范围也越大；基线法平面与定位平面交线对应的相位差为0°。

2）模糊相位差φ′是真实相位差φ向(-π,π]区间内塌缩的结果。由φ＝0的区域向外，每当φ超出(-π,π]范围时，就加上或减去2π的整数倍，使其值保持在(-π,π]区间中。

3）定位区域被划分为一系列具有相同模糊数n的区域，区域的边界由φ′＝π对应等相位差线确定，这些区域按照-N到N的顺序沿基线矢量方向依次排列，随着基线的旋转不断变化。

将定位区域上具有相同模糊数的区域称为“条带”。随着基线的旋转，目标位置所在条带不断变化。将每个测量时刻目标所在的条带n连接起来，就构成了一条“路径”。图1显示一条可能的路径在条带空间的变化情况。

图1 每个测量时刻目标所在条带构成的路径

1.2 观测过程的HMM模型

目标所在条带与其位置和基线转角有关，将条带作为隐变量，旋转基线相位差序列的观测过程可用HMM来描述。

HMM是一种时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程，其中隐藏的随机序列称为状态序列,观测产生的随机序列称为观测序列。

设Q＝{q1,q2,q3,…,q K}代表所有状态的集合，V＝{v1,v2,v3,…,v M}代表所有可能观测的集合，HMM可由如下三元组来描述：

式中，A＝[aij]K×K为状态转移概率矩阵，aij＝P(st＋1＝qj|st＝qi)表示从时刻t到时刻t＋1由状态qi转移到状态qj的概率；B＝[bij]K×M为观测概率矩阵，bij＝P(ot＝v j|st＝qi)表示时刻t在qi状态下观测到v j的概率；Π＝[πi]K×1为初始状态概率向量，πi＝P(s1＝qi)代表t＝1时刻处于状态qi的概率。

具体到旋转基线干涉仪定位的情况，目标所在条带n是不可观测的隐藏状态，状态集合Q＝{-N,-(N-1),…,N-1,N}。观测集合对应相位差测量值，是(-π,π]之间的连续变量。状态转移矩阵A用于描述基线旋转导致的目标在各“条带”之间的转移概率。观测概率矩阵B表示目标位于某一“条带”时观测到测量相位差的概率。Π代表目标位置最初所在条带的概率。

由上文可知，任意一条“路径”对应一个状态序列s1s2s3…sT，如果能够由测量相位差序列o1o2o3…o T推断出目标真实位置对应的状态序列，就实现了相位差解模糊。

接下来，通过将“卡尔曼”滤波(Kalman filtering)过程嵌入到“维特比”算法(Viterbi algorithm)中，同步实现相位差动态解模糊和滤波定位。

2 “维特比-卡尔曼”滤波定位算法

2.1 假设真实相位差已知情况下的“卡尔曼”滤波定位过程

“卡尔曼”滤波[9]是一种经典的递归滤波器，典型应用是从一组包含噪声的观测序列中预测物体的运动状态。在旋转干涉仪定位中，假设真实相位差已知，卡尔曼滤波定位过程如下。

1）状态方程和观测方程

为了讨论方便，假定目标固定不动，对应的状态方程为：

状态S＝rT，rT＝[x y z]T是目标的直角坐标矢量。

假设能够测量到真实相位差φ，观测方程为：φ＝

式中，d为干涉仪基线长度，b＝[x b yb zb]T为单位基线矢量，r1＝[x1y1z1]T是干涉仪所在位置的直角坐标矢量。

2）卡尔曼滤波过程

根据以上的状态方程和观测方程，卡尔曼滤波的具体过程如下。

预测：

最小预测均方误差(MSE)矩阵：

卡尔曼增益：

修正：

最小MSE矩阵：

式中，σ2[n]为真实相位差φ的测量噪声方差，H[n]为观测方程φ＝h(S)的雅克比矩阵

滤波过程中v(n)＝φ[n]-h(Sˆ[n|n-1]称为“新息”，表示观测值与预测值的差。

令ˆ[-1|-1]等于目标的初始位置[x0y0z0]T，M[-1|-1]等于目标初始位置方差矩阵为相位差φ的观测方差，通过递归计算，就可以得到对系统状态S的估计。

3）雅克比矩阵H[n]的推导

雅克比矩阵H[n]是真实相位差φ关于目标位置rT＝[x y z]T的导数：

式中，

令：

得：

以上是假设真实相位差φ已知情况下的卡尔曼滤波过程，但实际中干涉仪只能输出测量相位差φ′，要完成滤波定位过程，就需要实时解出目标所在的条带n，接下来，通过引入“维特比”算法来解决这一问题。

2.2 “维特比”算法解码最大概率路径

“维特比”算法是HMM模型中根据观测序列解码状态序列的标准算法。给定观测序列o1o2o3…o T，要找出最可能的状态序列s1s2s3…sT，使式（16）中的条件概率最大：

根据马尔可夫齐次性假设和观测独立性假设，等价于:

如果直接计算式（17），计算复杂度将是O(K T)，式中,K为状态数，T为序列长度，计算量非常大。“维特比”算法通过在“篱笆网络（Lattice）”中应用动态规划技术，将计算复杂度降低为O(TK2)，有效解决了状态解码问题。

图2给出了篱笆网络的结构。路径在网络中从左向右传播，由于所有路径必然经过网络中的每一列，且如果最大概率路径经过节点X ij，那么该路径上从起始时刻到第i个时刻的部分路径也必然是该时段所有路径中的最大概率路径。因此，在任一时刻t，只要考虑K条最大概率路径即可[10]。

图2 篱笆网络（Lattice）的结构

基于以上性质，篱笆网络中最大概率的路径的计算过程如下：

1）初始化。计算第一列各状态节点的初始概率δ1，并记录每个节点前一时刻的状态ψ1：

式中，πk为状态k的初始概率，bk(o1)为状态k条件下观测到o1的概率。由于第一列没有前一时刻，因此ψ1(k)统一记为0。

2）递推。对于t＝2,3,…,T，根据下式递推计算篱笆网络中每一列的K个节点的概率δt，并记录每个节点前一个时刻的状态ψt：

式中，ajk为状态j转移到状态k的概率，bk(ot)为状态k条件下观测到ot的概率。

3）终止。当递推到达时刻T后，递推终止，此时，按照下式计算最优路径概率P*和最优路径的最后一个状态k*T：

4）回溯。找到最优路径的最后一个状态k*T后，根据下式回溯出最大概率路径：

最后求得最大概率路径I*：

接下来，通过将“卡尔曼”滤波步骤嵌入到“维特比”算法过程中，扩展上述递归过程，完成状态转移概率A＝[aij]K×K和观测概率B t＝[bi(ot)]K×1的计算，同步实现旋转基线干涉仪相位差解模糊和滤波定位。

2.3 “维特比-卡尔曼”滤波过程

2.3.1 初始化篱笆网络和滤波状态

根据式（18）、式（19），初始化篱笆网络的状态，由于目标的位置未知，信号到达的初始时刻基线的旋转角度也是随机的，因此，可以简单地认为目标位于各个条带的概率是相等的，即πi＝1/k,i＝1,2,…,K。同样，由于目标在条带中的具体位置未知，观测概率的最初状态概率也认为是相等的，bk(o1)＝1/k,i＝1,2,…,K，因此：

在篱笆网络中，初始时刻的每个状态代表目标所在的条带，也就是可能的模糊数n∈[-N,N]。根据测量相位差φ′和模糊数n可以恢复出K＝2N＋1个可能的真实相位差φ(k)，利用φ(k)根据式（6）-（10）执行一步滤波定位，得到初始时刻每个状态下的一次定位位置S1(k)和最小均方误差矩阵M1(k)。

2.3.2 递推计算滤波位置和最大概率路径

注意到“卡尔曼”滤波和“维特比”最大概率路径计算过程都采用递归的方式计算，将2种过程融合在一起，用篱笆网络的状态k将测量相位差φ′转换为真实相位差φ(k)，用其进行一步滤波定位，同时，用该滤波定位结果对应的估计相位差φˆ(k)与真实相位差φ(k)之间的相似程度来计算观测概率，完成一步递归计算，具体步骤如下：

1）执行一步滤波定位

对于篱笆网络中的节点X tk，利用该节点对应的模糊数n将测量相位差φ′t转换为真实相位差φt(k)：

分别以t-1时刻的K个节点对应的滤波位置St-1(k)，k＝1,2,…,K和最小均方误差矩阵M t-1(k),k＝1,2,…,K为状态，φt(k)为观测值，根据式（6）-（10）执行一步滤波定位，得到节点X tk下的K个滤波位置Stk(j),j＝1,2,…,K和M tk(j),j＝1,2,…,K。

2）计算状态转移概率

状态转移概率矩阵A＝[aij]K×K，其中aij表示随着基线的旋转，在相邻时刻目标在从条带i转移到条带j的概率。根据在时刻t对系统状态了解程度的不同，按照以下几种情况来计算状态转移概率。

①无任何先验知识的情况

在目标位置无任何先验知识，且时刻t-1到t之间的时间间隔随机产生的情况下，随着基线的旋转，目标位置可能在任意2个条带之间转移，此时，可以简单认为转移概率是平均分布的，即：

②目标粗略位置已知的情况

假设已经获取了目标的粗略位置，例如在滤波定位过程中不断接近目标真实位置，那么，在t-1时刻，根据位置利用式（2）和式（5）可以估计出目标可能的条带范围，例如连续的3个条带(i-1,i,i＋1)，基线旋转到t时刻，同样可以计算出另一个可能的条带范围，例如(j-1,j,j＋1)，这样，从t-1时刻到t时刻，转移概率可以表示为：

具体的条带范围需要根据估计位置与真实位置的接近程度来设置，越接近，则可能的条带范围越小，否则，需要给多个可能的条带赋予非零概率。

③相对相位差无模糊的情况

当t-1和t时刻的真实相位差满足|φt-φt-1|≤π时，称为相位差相对无模糊。此时，可以恢复出测量相位差φ′t＋1和φ′t之间真实的差值。文献[1]中证明，当|Δω|≤2arcsin(λ(4d))时，干涉仪覆盖范围内任意位置的相位差相对无模糊。此时，转移概率可以表示为：

当φ′t-φ′t-1≤-π时，

当φ′t-φ′t-1＞π时，

当-π＜φ′t-φ′t-1≤π时，

3）计算观测概率

观测概率bk(ot)表示t时刻在状态k下观测到ot的概率。相位差是由位置决定的，经过一次滤波后，在节点X tk下有K个滤波位置Stk(j)，j＝1,2,…,K，根据式（5）可计算出对应的K个相位差：

将φtk(j)与真实相位差φt(i)之间的相似程度作为计算观测概率bk(ot)的依据。首先利用sigmoid函数将φtk(j)-φt(i)变换到(-1,1)的范围内：

然后按照标准正态分布计算观测概率：

4）保存最大概率状态

对于进入节点X tk的K条路径，选取当前最大的部分路径概率对应的一条作为X tk的状态。根据下式计算篱笆网络中时刻t的K个节点的概率δt，并记录每个节点前一个时刻的状态ψt：

记录最大概率路径对应的滤波位置和最小均方误差矩阵：

5）递归终止和路径回溯

当递推到达时刻T后，递推终止，此时，按照下式计算最优路径概率P*、最优路径的最后一个状态k*T和滤波定位位置S*T：

找到最优路径的最后一个状态k*T后，根据下式回溯出最大概率路径：

最后求得最大概率路径I*：

至此，完成了“维特比-卡尔曼”滤波定位的全过程。该算法利用“维特比”算法根据测量相位差序列解码最优路径，相当于解出每次测量时对应的模糊数，从而利用真实相位差完成滤波定位。定位性能方面与“卡尔曼”滤波一致，此处不再赘述。

3 仿真分析

仿真场景设置如下：目标位置的经纬度为[123.117 5°,23.553 1°]，旋转基线位置的经纬度为[125°,28°]，在500 km高度水平放置，基线长度为2 m，按照360°/s的速度旋转。辐射源信号脉冲重复间隔在[1 000µs,10 000µs]区间内随机选取，脉冲列长度为100。滤波初始位置经纬度随机选取为[124.132 7°,22.080 9°]。

采用“维特比-卡尔曼”滤波算法进行处理，经过100步滤波后得到定位位置的经纬度为[123.114 7°，23.551 3°]，定位误差为351 m，与利用真实相位差进行卡尔曼滤波后的定位结果一致。滤波定位每一步的结果如图3所示。其中，图3(a)给出了滤波定位过程中估计位置不断接近真实位置的情况；图3(b)是滤波过程中最大概率路径对应的相位差“新息”的变化情况，可见，随着滤波的进行，相位差“新息”在不断减小，说明估计相位差与真实相位差越来越接近；图3(c)给出了篱笆网格中由最终状态回溯的所有路径（13条）的每一步的定位结果，可见，非真实路径对应的定位结果始终是发散的，而真实路径对应的结果逐渐收敛到真实位置（红色星号）；图3(d)给出了在维特比递推过程中篱笆网络中每一列最大概率的变化情况，可见，随着递归的进行，真实路径上的概率比重不断增加，而其它路径上的概率由于采用了错误的相位差进行滤波，无法收敛到真实位置，导致概率比重逐渐减小，通过选取最大概率的状态并回溯后，找出真实路径，从而实现相位差解模糊和定位。