4.1 原理与算法

4.1.1 原理与估计准则

针对静止辐射源，二维测向定位基于多个方向线交汇进行位置估计，图4-3（a）是无测向误差定位情况，图4-3（b）是存在测向误差定位情况。三维空间测向定位基于多个方位面与俯仰圆锥构成的相交线交汇进行位置估计，可利用运动的单个观测站或多个分布在不同位置的观测站来测量辐射源方向信息，结合观测站自身位置，通过优化算法确定辐射源的位置。

以三维空间测向定位为例，假设有坐标为 x ₁ =[ x ₁ ， y ₁ ， z ₁ ] ^T ， x ₂ =[ x ₂ ， y ₂ ， z ₂ ] ^T ，…， x _L =[ x _L ， y _L ， z _L ] ^T 的观测站对某固定的目标辐射源 x _T =[ x _T ， y _T ， z _T ] ^T 进行定位。此时测得的方位角和俯仰角分别为 _l α 和 _l β ， l =1，2， … ， L ，则存在以下关系式：

测向定位的观测方程可表示为

式中，观测量 z = ，其中方位角 构成方位面，俯仰角 构成圆锥面；观测站的位置 x _O =[ x ₁ ， x ₂ ， … ， x _L ] ^T ， x _l =[ x _l ， y _l ， z _l ] ^T 为第 l 个观测站的位置或单个运动观测站在第 l 个观测时刻的位置， l =1，2， … ， L ； h （ x _T ， x _O ）= ， =arctan 2[（ y _T -y _l ），（ x _T -x _l ）]， =arctan ；观测噪声 ξ =[Δ α ₁ ，Δ β ₁ ，Δ α ₂ ，Δ β ₂ ， … ，Δ α _L ，Δ β _L ] ^T ，假设 = 0 ，协方差矩阵为 Σ _ξ = 。

可利用式（3-6）的广义最小二乘估计目标位置：

在一定条件下，方向测量的误差服从高斯分布 ^[1] ，此时式（4-3）也是辐射源位置的极大似然估计，且极大似然估计的 近似等于tr[CRLB（ x ）]。

4.1.2 优化算法

4.1.2.1 网格搜索法

对于式（4-3）的求解，在覆盖目标真实位置的区域范围内划分网格，通过比对网格点上的理论来波方向与实际测向结果的偏差值，遍历所有网格最终得到目标位置估计。除了基于理论来波方向进行比对，若采用干涉仪测向体制，则可利用理论相位差数据与实际相位差结果进行比对。下面介绍一种基于相位差数据的网格搜索法。

文献[2-4]提出最小相位误差的定位方法，即不测量辐射源的角度信息，利用多次测量相位差并使其在最小二乘意义下拟合误差最小化。另外，对不同基线相位差赋予不同的权值，可以避免某一基线的相位差误差较大时影响辐射源位置估计精度的情形。

设静止辐射源位置坐标为 x _T =[ x _T ， y _T ， z _T ] ^T ，考虑单个运动观测站进行 K 次观测，观测站第 k 次观测时的位置为 x _k =[ x _k ， y _k ， z _k ] ^T ，令观测站位置 x _k 至辐射源位置 x _T 的矢量为 r _k = x _T -x _k 。观测站的测向系统为 N 元天线阵， N 元天线阵最多可以组成 N （ N -1）/2条不同的基线，假设共选取 M （ M ≤ N （ N -1）/2）条不同的基线，则观测站第 k 次观测的第 m 条基线的相位差可以表示为

式中， a _{k ， m} 表示第 k 次观测的第 m 条基线两阵元位置构成的指向矢量（假设参考阵元为（0，0，0） ^T ，则阵元 x _m =（ x _m ， y _m ， z _m ） ^T 与参考阵元间的指向矢量 a _m = x _m ）， λ 表示辐射源信号波长。

当观测站的 M 条测向基线长度均小于 λ 2时，不同位置的辐射源信号在观测站构成的一组相位差也不相同，因此观测站可视区内不同位置的集合 Ω _x 与不同角度入射到观测站天线阵的相位差矢量集 Ψ 之间存在一一对应关系，通过 Ψ 到 Ω _x 的逆映射可以对辐射源进行定位，这就是最小相位误差无源定位的原理。当测向基线大于半波长并存在测向模糊的时候，不同位置的辐射源可能对应于同一组相位差矢量，出现定位模糊问题。对静止辐射源，如果观测站对辐射源进行了 K （ K >1）次相位差的测量，那么真实的辐射源位置都是这 K 次相位差矢量的原像，而虚假定位点则由于观测站运动产生的相位差非线性变化，并不总是这 K 次相位差矢量的原像，因此通过使辐射源入射信号的相位差测量值与理论值的误差平方和最小化可以去模糊定位，这就是最小相位误差单站无源定位去模糊的原理。

用 表示 φ _{k ， m} 的测量值， ξ _{k ， m} 表示测量误差。写成矩阵形式为

式中， z = ； h （ x _T ）=[ h _1，1 ， h _1，2 ， … ， h _{1， M} ， … ， h _{K ，1} ， h _{K ，2} ， … ， h _{K ， M} ] ^T ； ξ =[ ξ _1，1 ， ξ _1，2 ， … ， ξ _{1， M} ， … ， ξ _{K ，1} ， ξ _{K ，2} ， … ， ξ _{K ， M} ] ^T ， ，协方差矩阵为 Σ _ξ = 。

则利用最小二乘法估计辐射源位置 ，可以表示为

利用最小相位误差进行定位运算需要在目标可能的范围内进行搜索，结合观测站的航向角、俯仰角、横滚角，分别计算天线阵各测向基线的理论相位差，并与实际收到的相位差按式（4-6）寻优。可以先对目标可能存在的区域进行粗搜，然后在得到的粗定位结果附近减小步进进行精搜。

【例4.1】 双站测向定位二维目标

两个观测站坐标为 x ₁ =[0，0] ^T （km）， x ₂ =[40，0] ^T （km），目标辐射源位置为 x _T =[32，33] ^T （km）。理论方位角为 z =[ α ₁ ， α ₂ ] ^T =[45.88°，103.63°] ^T ，假设两个观测站的测向误差 σ _α =2°，得到的测向结果为 =[45.65°，101.81°] ^T 。基于相位差数据，利用2km的间隔对目标函数式（3-15）进行粗搜，其结果如图4-4（a）所示，得到目标位置粗搜估计值为 =[35，35] ^T （km），定位误差为3.61km。在粗搜的峰值附近再利用0.01km的间隔进行精搜，其结果如图4-4（b）所示。提取精搜的峰值得到目标位置估计值为 =[32.95，33.71] ^T （km），定位误差为1.19km。

4.1.2.2 G-N迭代法

针对式（4-2）中的非线性函数 h （ x ），可利用G-N迭代法求解 ^[5] 。基于一阶泰勒级数展开的数值迭代为

式中， = ， = + ， + + 。

令 δ = ，并重复上述过程，直到‖ δ ‖ < ε 为止， ε 为设定的一个较小的误差门限值。

利用G-N迭代法对例4.1求解，迭代初始值选为 =[5，15] ^T （km）， ε =1，则定位结果如图4-5所示，经过3次迭代得到目标位置的估计值为 =[32.95，33.71] ^T （km），定位误差为1.19km。

4.1.2.3 解析法

1.加权Stansfield法

文献[6-8]介绍了三维空间测向定位，提出了加权Stansfield 算法和辅助变量加权Stansfield算法。当测向误差服从高斯分布时，由式（4-2）构建最大似然估计的代价函数为

式（4-8）可近似拆分为两个代价函数之和：

式中， 为使 P ₁ （ x _T ， y _T ）达到最小值时的取值。

当测向误差较小时，上述两个代价函数可以近似表示为

使代价函数 达到最小值的参数 即二维测向定位 ^[9] 的结果：

式中，
H = ； Σ _r = ， = ， =（ x _T -x _l ） ² +（ y _T -y _l ） ² ； b = 。

根据 可计算代价函数 ，即

式中， = 。

对 关于 z _T 求偏导，并令其等于0可得

式中， g _l = 。

加权Stansfield法是通过线性化得到的解析解，其结果为有偏估计。用加权Stansfield法对例4.1求解可得目标位置的估计值为 =[32.95，33.71] ^T （km），定位误差为1.19km，定位结果如图4-6所示。

2.总体最小二乘法

基于式（4-3）结合极大似然估计可得到目标位置的无偏估计值，但前提是需要获得目标位置的粗略初始估计值，随后进行迭代计算或网格搜索。文献 ^[9] 将总体最小二乘算法应用于无源定位中，首先将非线性观测方程转化成伪线性方程，构造增广矩阵，并对该矩阵进行奇异值分解即可估计出目标位置，无须初始估计值与迭代计算，具有较高的定位实时性。

对于观测方程式（4-2），用三角函数关系进行伪线性化，得到如下线性方程组：

式中， H = ， b = 。

假设测向误差Δ _l α 与Δ _l β 较小，满足sin（Δ α _l ）≈Δ _l α ，cos（Δ α _l ）≈1，sin（Δ β _l ）≈Δ _l β ，cos（Δ β _l ）≈1，Δ α _l ·Δ β _l ≈0，则式中，

令 D =[Δ H ，-Δ b ]，此时定位解算问题可转化为如下最小二乘问题：

式中， 为矩阵的Frobenious范数。

由于考虑了误差矩阵Δ H 与Δ b ，该算法为总体最小二乘算法 ^[9] ，其基本思想是使Δ H 和Δ b 噪声扰动影响最小，即求得具有最小范数的矩阵 D 使得增广矩阵[ H +Δ H ， -b -Δ b ]是非满秩，利用奇异值分解可实现这一目标。将测量矩阵和观测参数矢量进行奇异值分解 ^[10] ：

式中，对于三维空间测向定位， U 为 L ×4的矩阵； Σ =diag（ σ ₁ ， σ ₂ ， σ ₃ ， σ ₄ ）为由按序排列的奇异值（ σ ₁ ≥ σ ₂ ≥ σ ₃ ≥ σ ₄ >0）构成的对角线矩阵； V =[ v ₁ ， v ₂ ， v ₃ ， v ₄ ]为4×4的矩阵。

当没有测量误差时，有rank{[ H ， -b ]}=3，此时 σ ₄ =0。当测量矩阵和观测参数矢量均存在误差扰动Δ H 与Δ b 时，rank{[ H +Δ H ， -b -Δ b ]}=4，且 σ ₄ ≠0。但当误差扰动较小时，满足 σ ₁ ≥ σ ₂ ≥ σ ₃ ≫σ ₄ ，则通过令 σ ₄ =0可以近似得到扰动的增广矩阵：

总体最小二乘定位结果 x _TLS 需满足

可以得到总体最小二乘算法的定位结果为

式中， v ₄ （ i ）为 v ₄ 的第 i 个分量； v ₄ 为矩阵[ H +Δ H ， -b -Δ b ]进行奇异值分解后，最小的奇异值所对应的右奇异向量。

总体最小二乘法涉及伪线性变换，因此该方法的使用需要满足较小的参数测量误差。利用总体最小二乘法对例4.1构建矩阵[ H +Δ H ， -b -Δ b ]进行奇异值分解，得到目标位置的估计值为 =[32.95，33.71] ^T （km），定位误差为1.19km。

3.带约束条件的总体最小二乘法

针对二维测向定位，文献 ^[11，12] 指出理论测量矩阵 H 中每个行向量的模为1，而总体最小二乘法中带有误差的测量矩阵 H +Δ H 中行向量不满足这一条件。因此，将矩阵 H +Δ H 中每个行向量元素相对于自身进行归一化处理，使其满足每个行向量的模为1。该算法在总体最小二乘法的基础上增加了非线性约束，因此称为约束总体最小二乘法，具体算法流程如下。

（1）令测量矩阵 H _k = H +Δ H ，对参数矢量 b _k = b +Δ b 构成的增广矩阵进行奇异值分解：

（2）得到三个奇异值 σ ₁ ≥ σ ₂ ≥ σ ₃ >0，计算误差值：

若误差值小于一个门限值，则终止迭代，利用公式 =[ v ₃ （1） /v ₃ （3）， v ₃ （2） /v ₃ （3）] ^T 得到约束总体最小二乘法的定位结果。否则进入步骤（3）。

（3）由于理论测量矩阵 H 中行向量的模为1，则对观测矩阵中的行向量进行自身归一化：

（4）对参数矢量 b _k 基于 H _k 中对应行向量进行归一化调整：

（5）更新 H _{k +1} 和 b _{k +1} ，并回到步骤（1）。

该算法相较于总体最小二乘法的单次解析具有更高的定位精度，但每次迭代运算均需要进行奇异值分解，算法过程较为复杂。

针对总体最小二乘法单次解析定位精度不高的问题，文献[13-15]提出三维空间测向定位的约束总体最小二乘法，将测向定位问题转化成一个约束总体最小二乘问题。假设目标相对于第 l 个观测站的真实方位角为 _l α = -Δ _l α ，真实俯仰角为 _l β = -Δ _l β ，则考虑如下三角函数的一阶泰勒展开：

由此可得带误差的测量矩阵 H _c 与参数矢量 b _c 表达式如下：

式中， H 和 b 具体表达式参见式（4-15）； F ₁ =diag[ f ₁₁ ， f ₁₂ ， … ， f _{1 L} ]， F ₂ =diag[ f ₂₁ ， f ₂₂ ， … ， f _{2 L} ]， F ₃ =diag[ f ₃₁ ， f ₃₂ ， … ， f _{3 L} ]， F ₄ =diag[ f ₄₁ ， f ₄₂ ， … ， f _{4 L} ]， f _{1 l} = ， f _{2 l} = ， f _{3 l} = ， f _{4 l} = ； E =[Δ α ₁ ，Δ β ₁ ，Δ α ₂ ，Δ β ₂ ， … ，Δ α _L ，Δ β _L ] ^T 。

由上式可以看出，矩阵 H _c 与 b _c 所受到的噪扰均来源于误差矢量 E ，此时定位解算问题可转化为如下带约束条件的总体最小二乘问题：

由文献 ^[16] 可知，式（4-22）带约束条件的总体最小二乘问题可转化为无约束的优化问题，其优化表达式如下：

式中，优化过程中目标位置 x _T =[ x ， y ， z ] ^T ， H _x = x F ₁ + y F ₂ + z F ₃ -F ₄ ， G = ， C =[ H _c ， b _c ]。

式（4-23）是关于目标位置 x _T 的非线性函数，可利用G-N迭代法进行求解：

式中， H _k =2（ H _c -B ₁ -B ₂ ） ^T （ H _c -B ₁ -B ₂ ）- 为 Hessian 矩阵， T _k =2（ U ^T H _c -U ^T B ₁ ） ^T 为梯度矢量， U = ， B ₁ = ， B ₂ = ， B ₃ = ； μ _k = μ ^k （ μ <1）为步长因子。

该算法收敛时可有效提高总体最小二乘法的单次解析定位精度，但需要较多的观测量信息，同时涉及梯度矢量与Hessian 矩阵的计算，以及G-N迭代求解，计算过程较为复杂。在此算法的基础上，还有学者提出将带约束条件的总体最小二乘法转化成结构总体最小二乘法，感兴趣的读者可参考文献 ^[17-19] 。

4.基于正则约束总体最小二乘法

上述介绍的优化算法均利用G-N迭代法进行求解，在足够多的测量信息和较小的测量误差前提下才具有较好的定位性能，否则迭代优化算法难以收敛。文献[20，21]提出基于正则约束总体最小二乘法定位，即在式（4-22）的约束条件中引入正则项 η ，可得基于正则约束总体最小二乘法的等价优化问题：

类似于式（4-23）的无约束优化问题，式（4-25）同样可等价于如下无约束优化问题：

对于式（4-26）可使用G-N迭代法求解，但是需要计算梯度矢量与Hessian矩阵，计算复杂度高，求解时间长。文献 ^[20] 提出对式（4-26）求导并忽略误差矢量 E 的高阶项，从而得到如下近似解：

当 = I 时，则近似解 退化为Tikhonov正则化法闭式解 ^[22] 。

该优化算法与带约束条件的总体最小二乘法相比，在观测站数目较少和测向误差标准差较大时具有较高的定位精度。但该算法是有偏估计，偏差随着正则化参数 η 的增大而增大，正则化参数 η 的确定需要在偏差和均方差之间取舍。此外，该算法前提条件是需要确切知道角度测量信息的误差协方差矩阵，同时需要对正则化参数 η 进行动态调整，这也就限制了该优化算法的工程易用性。 N+fXjgNwBWwgzSQyO38mGULSVzjrQ8mXmTLfsZeaOweQXwUQ1wd36Lvmlh8MdDJS