1.   引言

短波通信是一种中远程超视距通信方式，不会受到网络枢纽节点和有源中继的制约，具有通信覆盖范围广、开通迅速、机动灵活、抗毁性能强、网络重构便捷等优点，是国防安全领域中的重要通信手段，在某些极端情况下甚至是唯一可用的通信手段。时至今日，短波通信依然是正向通信领域中的研究热点^[1]，因此研究针对短波辐射源的定位技术与方法具有重要意义。

在短波辐射源定位系统中最常使用的是多站角度交汇定位体制，该类系统需要每个观测站都安装天线阵列，并且通过阵列信号参数估计方法获得到达角度(Direction-Of-Arrival, DOA)参数，然后利用此观测信息估计目标辐射源位置坐标。时至今日，学者已经提出诸多不同类型的DOA定位方法^[2−7]，虽然这些方法的定位精度能渐近逼近克拉美罗界(Cramér-Rao Lower Bound, CRLB)，但大都是面向信号视距传播场景的定位方法。众所周知，短波信号主要通过天波散射传播，因此短波辐射源定位是在超视距场景下实现的，对此现有文献中的视距DOA观测模型已不再适用。在短波超视距定位场景中，信号到达方位角可以通过一个与观测站地理坐标相关的旋转矩阵来获得^[8,9]，其是一种线性变换，并且旋转矩阵是已知的，并不会对方位角观测模型带来本质影响。然而，信号到达仰角与电离层虚高观测紧密相关^[8−10]，因此仰角观测模型的非线性程度会得到显著增强。文献[8,9]针对短波超视距定位场景，通过推导一元二次方程的根以及引入辅助变量的方式得到仰角伪线性观测方程，并基于此提出两种迭代型DOA定位方法，他们的估计性能均可以逼近CRLB，但都需设置迭代初始值，因此存在局部收敛和迭代发散的风险。

另一种常用多站定位体制是时差(Time-Difference-Of-Arrival, TDOA)交汇定位，该类定位系统并不需要观测站安装天线阵列。针对视距场景下TDOA观测方程的非线性特性，现有TDOA定位方法大都将其转化为伪线性观测方程进行处理^[11−13]。然而，在短波超视距定位场景下，由于信号传播路径涉及电离层参数，此时TDOA观测方程具有更强的非线性特征，难以直接进行伪线性处理。近年来，短波超视距场景下的TDOA定位方法得到越来越多学者的关注。文献[14]研究短波信号TDOA估计方法，其实验结果表明，利用分布式传感器可以对2000 km远的短波信号提取TDOA参数，并具有较高精度。文献[15,16]利用准抛物线(Quasi-Parabolic, QP)电离层模型获得短波信号超视距传播的距离表达式，并基于此构建同时含有等式约束和不等式约束的非线性优化模型，针对该模型分别提出广义投影梯度下降算法、坐标下降算法等进行寻优计算。针对文献[15,16]中的定位方法易受迭代初始值影响而陷入局部最优解的问题，文献[17]提出一种与粒子群优化相结合的协作投影梯度方法，该方法能进行全局寻优，但需要更高计算复杂度。文献[18]针对天波超视距多输入多输出雷达，提出基于QP模型的到达时间(Time-Of-Arrival, TOA)定位方法，该方法将爬山算法与加权最小二乘估计相结合，通过3个计算阶段实现对短波辐射源的高精度定位。此外，文献[19]利用国际参考电离层(International Reference Ionosphere, ISI)模型，基于数值型三维射线追踪方法建立短波TDOA定位观测方程，推导相应的CRLB，并基于此进行理论性能分析。

QP模型和ISI模型都能以射线追踪方程为基础建立短波TDOA定位观测方程，但需要一些电离层参数(包括电离层离子密度、介电常数等)，在实际应用中难以准确获知，并且当电离层变化较复杂时，这些模型的定位精度极易受到影响。事实上，由于电离层时空变化特性无法精确获知，为了进行短波TDOA定位，在大多数情况下只能对电离层做近似处理。文献[20,21]提出一种能刻画信号传播距离的(单跳)电离层虚高模型，其中将电离层和地面近似为平行平面，定位观测方程相对简单，随着目标距离的增加，其定位精度会有所恶化。单跳模型可以覆盖3000 km远的通信距离，文献[22]提出一种电离层环境未知条件下的短波TDOA定位方法，该方法通过求解半正定规划模型获得短波辐射源位置参数，并且无需电离层先验信息，由于其对定位模型进行较大幅度的松弛，因此定位精度有时难以达到CRLB。文献[23]在假设短波辐射源信号到达各个观测站电离层虚高相等的条件下，提出一种基于网格选择的短波TDOA定位方法，该方法包含测试和搜索两个阶段，能扩展至电离层虚高未知的场景，获得了较优的定位精度，但需要较高的计算复杂度。

综上所述，短波多站DOA定位和短波多站TDOA定位是短波频段最具代表性的两种多站定位体制。他们各有优势，前者在信号参数估计环节无需站间数据交互，因此受同步误差影响相对较小；后者则对远距离目标具有更高定位精度。事实上，若能融合这两种体制进行协同定位，则会带来显著的协同增益。近些年来，学者虽然提出很多DOA/TDOA协同定位方法^[24−26]，但基本都是针对视距场景提出的，难以直接应用于短波超视距定位场景。文献[8]针对短波目标辐射源提出一种联合2D-DOA观测量和TDOA观测量的协同定位方法，取得了较高的定位精度，但该方法假设每个观测站仅能获得2D-DOA和TDOA中的一种参数。然而，利用现代阵列信号参数估计方法可以利用每个观测站的阵列天线联合估计2D-DOA和TDOA两种参数，在此基础上进行协同定位能获得更高协同增益。基于上述讨论，本文基于单跳电离层虚高模型，提出一种联合2D-DOA观测量和TDOA观测量的协同定位新方法。与文献[8]考虑的场景不同，本文假设每个观测站能联合估计信号2D-DOA/TDOA两域参数，因此需要每个观测站均安装天线阵列，减少了所需观测站个数下界。本文的主要工作和创新点包括：(1)针对短波超视距定位场景，利用多项式求根和拉格朗日乘子技术提出一种两阶段闭式解型2D-DOA/TDOA协同定位新方法，其无需迭代初始值，可有效避免局部收敛的影响；(2)利用约束误差扰动理论和正交投影矩阵的数学性质推导2D-DOA/TDOA协同定位新方法的估计均方误差理论表达式，并证明新方法的定位精度可以渐近逼近CRLB；(3)在短波辐射源高度信息误差存在下，基于约束误差扰动理论推导新方法的定位均方误差的理论表达式，并基于此提出能确保地球椭圆约束产生性能增益的短波辐射源高度信息误差的最大门限值。

2.   短波超视距2D-DOA/TDOA协同定位模型与克拉美罗界

2.1   定位几何模型

图1描绘了短波超视距定位几何模型，现在地球表面存在一个待定位的短波辐射源，其经度和纬度分别为${\omega _{\text{o}}}$和${\rho _{\text{o}}}$。根据文献[20−23]中的讨论可知，短波信号进入电离层和离开电离层的两条路径可以利用(单跳)平面反射模型来表征，因此存在一个等效电离层反射平面，反射面距离地面的高度即为电离层虚高。假设有N个观测站可以同时接收该短波信号，其中第n个观测站的经度和纬度分别为${\omega _{{\text{s,}}n}}$和${\rho _{{\text{s,}}n}}\;(1 \le n \le N)$，每个观测站均安装天线阵列，利用现代阵列信号参数估计技术可以联合估计辐射源信号2D-DOA/TDOA两域参数。

图  1  短波辐射源2D-DOA和TDOA几何模型示意图

Figure  1.  Schematic diagram of 2D-DOA and TDOA geometric models of shortwave radiation source

下载: 全尺寸图片 幻灯片

2.2   2D-DOA观测模型

将短波信号到达第$n\;(n \ge 1)$个观测站的方位角和仰角分别记为${\theta _n}$和${\beta _n}$，地面辐射源与第n个观测站在地心地固坐标系(Earth-Centered, Earth-Fixed, ECEF)下的位置向量分别如式(1)所示^[27,28]：

其中，$e = {\text{0}}{\text{.081819790992113}}$表示地球第一偏心率；${r_{\text{e}}} = {\text{6378}}{\text{.137 km}}$表示地球赤道半径。根据文献[8]中的讨论可知方位角${\theta _n}$的数学模型为

其中，${{\boldsymbol{S}}_n}$表示三维空间旋转矩阵，其表达式见文献[8]；${{\boldsymbol{s}}_{n{\text{,1}}}}$和${{\boldsymbol{s}}_{n{\text{,2}}}}$分别表示矩阵${\boldsymbol{S}}_n^{\text{T}}$中的第1列向量和第2列向量。定义方位角向量${\boldsymbol{\theta }} = {[{\theta _{\text{1}}}\;\;{\theta _{\text{2}}}\;\; \cdots \;\;{\theta _N}]^{\text{T}}}$，在实际中仅能得到其观测值${\boldsymbol{\hat \theta }}$，相应的观测模型如式(3)所示：

其中，${\boldsymbol{\theta }} = {{\boldsymbol{g}}_\theta }({{\boldsymbol{u}}_{\text{o}}})$表示由式(2)定义的非线性函数；向量${{\boldsymbol{\xi }}_\theta }$表示方位角观测误差，其服从零均值的高斯分布。

另外，将信号到达第n个观测站对应的电离层虚高记为${d_n}$，根据文献[8]中的讨论可知仰角${\beta _n}$的数学模型为

其中，$\phi_n=\arcsin(||\boldsymbol{u}_{\text{o}}-\boldsymbol{u}_{\text{s,}n}||_2/\left(2r_{\text{o}}\right))$；${r_{\text{o}}} = {\text{6371}}{\text{.393 km}}$表示地球平均半径。定义仰角向量${\boldsymbol{\beta }} = {[{\beta _{\text{1}}}\;\;{\beta _{\text{2}}}\;\; \cdots \;\;{\beta _N}]^{\text{T}}}$，在实际中仅能得到其观测值${\boldsymbol{\hat \beta }}$，相应的观测模型如式(5)所示：

其中，${\boldsymbol{d}} = {[{d_1}\;\;{d_2}\;\; \cdots \;\;{d_N}]^{\text{T}}}$表示电离层虚高向量；${\boldsymbol{\beta }} = {{\boldsymbol{g}}_\beta }({{\boldsymbol{u}}_{\text{o}}},{\boldsymbol{d}})$表示由式(4)定义的非线性函数；向量${{\boldsymbol{\xi }}_\beta }$表示仰角观测误差，其服从零均值的高斯分布。

注释1：由式(1)可知，短波辐射源位置向量${{\boldsymbol{u}}_{\text{o}}}$满足式(6)二次椭圆等式约束

显然，向量${{\boldsymbol{u}}_{\text{o}}}$的自由度应为2，不妨将其前面两个元素构成的向量记为${{\boldsymbol{\bar u}}_{\text{o}}}$，则第3个元素可以写成关于向量${{\boldsymbol{\bar u}}_{\text{o}}}$的函数，即${[{{\boldsymbol{u}}_{\text{o}}}]_3} = \pm \sqrt {(1 - {e^2})(r_{\text{e}}^2 - {\text{||}}{{{\boldsymbol{\bar u}}}_{\text{o}}}{\text{||}}_2^2)}$，这意味着向量${{\boldsymbol{u}}_{\text{o}}}$可以由向量${{\boldsymbol{\bar u}}_{\text{o}}}$来表示。

2.3   TDOA观测模型

假设第1个观测站为主观测站，其余均为辅助观测站，将短波信号到达第$n\;(n \ge 2)$个观测站与到达主观测站的TDOA记为${t_n}$，根据文献[8]中的讨论可知TDOA${t_n}$的数学模型为

其中，${\mathrm{c}} = {{3 \times 1}}{{\text{0}}^{\text{5}}}\;{\text{km/s}}$表示信号传播速度；${a_n} = 2r_{\text{o}}^{\text{2}} + 2{r_{\text{o}}}{d_n} + d_n^2$；${b_n} = 2{r_{\text{o}}}({r_{\text{o}}} + {d_n})$。定义TDOA向量${\boldsymbol{t}} = {[{t_{\text{2}}}\;\;{t_{\text{3}}}\;\; \cdots \;\;{t_N}]^{\text{T}}}$，在实际中仅能得到其观测值${\boldsymbol{\hat t}}$，相应的观测模型如式(8)所示：

其中，${\boldsymbol{t}} = {{\boldsymbol{g}}_t}({{\boldsymbol{u}}_{\text{o}}},{\boldsymbol{d}})$表示由式(7)定义的非线性函数；向量${{\boldsymbol{\xi }}_t}$表示TDOA观测误差，其服从零均值的高斯分布。

2.4   2D-DOA/TDOA协同定位观测模型

将式(3)、式(5)以及式(8)进行合并可以得到短波超视距2D-DOA/TDOA协同定位观测模型，如式(9)所示：

其中

其中，向量${{\boldsymbol{\xi }}_z}$表示2D-DOA/TDOA观测误差，其服从零均值的高斯分布，协方差矩阵为${{\boldsymbol{\varSigma }}_z} = {\bf{E}}[{{\boldsymbol{\xi }}_z}{\boldsymbol{\xi }}_z^{\text{T}}]$。

式(9)中的短波超视距2D-DOA/TDOA协同定位观测模型与电离层虚高向量d紧密相关，实际中可以通过垂直探测、后向散射反演等技术得到其先验观测值。电离层虚高观测值与其真实值之间通常存在观测误差，因此可以将向量d的观测模型表示为

其中，向量${{\boldsymbol{\xi }}_d}$表示电离层虚高观测误差，其服从零均值的高斯分布，协方差矩阵为${{\boldsymbol{\varSigma }}_d} = {\bf{E}}[{{\boldsymbol{\xi }}_d}{\boldsymbol{\xi }}_d^{\text{T}}]$。需要指出的是，尽管电离层具有实时变化性，但是在短波频段完成单次定位所需要的信号采集与参数解算仅为毫秒级时间，因此可以假设在定位计算过程中电离层虚高观测量${\boldsymbol{\hat d}}$保持不变。

2.5   克拉美罗界

虽然文中考虑的短波超视距协同定位场景与文献[8]中的协同定位场景并不完全一致，但可以在相同框架下推导两种CRLB，他们也具有相似的形式。两者的主要区别在于文献[8]中的2D-DOA观测和TDOA观测源自不同观测站，因此协方差矩阵${{\boldsymbol{\varSigma }}_z}$具有块状对角结构，而本文假设每个观测站联合估计信号2D-DOA/TDOA两域参数，此时协方差矩阵${{\boldsymbol{\varSigma }}_z}$不再具有块状对角结构。

实际中仅关心短波辐射源位置向量${{\boldsymbol{u}}_{\text{o}}}$的CRLB矩阵，基于文献[8]中的结论可得

其中，${\bf{CRL}}{{\bf{B}}_{{\text{un}}}}({{\boldsymbol{u}}_{\text{o}}})$表示没有等式约束条件下向量${{\boldsymbol{u}}_{\text{o}}}$的CRLB矩阵，其表达式为

其中，$\boldsymbol{G}_{z\text{,}u} (\boldsymbol{u}_{\text{o}},\boldsymbol{d}) = \partial\boldsymbol{g}_z(\boldsymbol{u}_{\text{o}},\boldsymbol{d}) / \partial\boldsymbol{u}_{\text{o}}^{\text{T}}$和${{\boldsymbol{G}}_{z{\text{,}}d}}({{\boldsymbol{u}}_{\text{o}}},{\boldsymbol{d}}) = {{\partial {{\boldsymbol{g}}_z}({{\boldsymbol{u}}_{\text{o}}},{\boldsymbol{d}})}}/{{\partial {{\boldsymbol{d}}^{\text{T}}}}}$。另一方面，如果定义矩阵$\boldsymbol{J}(\boldsymbol{u}_{\text{o}})= \left[\boldsymbol{I}_2\; \; (\mathrm{\mathit{e}}^2-1)\boldsymbol{\overline{u}}_{\text{o}}/[\boldsymbol{u}_{\text{o}}]_3\right]^{\text{T}}$，则还能得到矩阵${\bf{CRLB}}({{\boldsymbol{u}}_{\text{o}}})$的另一种形式，如式(14)所示：

为了便于后续理论性能分析，这里还需要获得关于向量${{\boldsymbol{\bar u}}_{\text{o}}}$的CRLB矩阵，如式(15)所示：

结合式(14)和式(15)可得

3.   短波超视距2D-DOA/TDOA协同定位新方法

相比视距场景下的2D-DOA和TDOA观测模型，在短波超视距场景下其定位观测模型的非线性程度显著增强，此时对迭代初始值也会更加敏感，因此本文提出一种闭式解型定位方法，其中包括两个计算阶段，阶段1需要将2D-DOA/TDOA非线性观测方程转化为伪线性观测方程，并通过求解一元六次多项式的根获得目标位置向量的闭式解，阶段2则通过构建等式约束优化模型对阶段1的估计误差进行改良，并且利用拉格朗日乘子技术得到精度更高的定位结果。下面描述新方法的基本原理和计算步骤。

3.1   阶段1的计算原理与方法

3.1.1   伪线性观测方程

首先给出方位角伪线性观测方程，其可以直接由式(2)获得^[8]，如式(17)所示：

其中

接着推导仰角伪线性观测方程。需要指出的是，文献[8]利用一元二次方程的根获得仰角伪线性观测方程，其计算过程相对复杂，本文则从几何关系着手进行分析，可以省略一元二次方程求根环节。在图1$\bigtriangleup {{ABC}}$中利用正弦定理可知：

结合式(4)中${\phi _n}$的定义和式(19)可得

不妨将标量$||{{\boldsymbol{u}}_{\text{o}}}{\text{||}}_{\text{2}}^{\text{2}}$作为辅助变量，此时综合式(20)中的N个方程可以将仰角伪线性观测方程写为

其中

最后推导TDOA伪线性观测方程。文献[8]假设每个观测站仅能获得2D-DOA和TDOA中的一种观测量，由于短波超视距TDOA观测模型的高度非线性特征，难以直接对其进行伪线性化处理，需要借助中间定位结果。然而，本文假设每个观测站能联合估计信号的2D-DOA/TDOA两域参数，此时在2D-DOA观测量的协同下可以直接得到TDOA伪线性观测方程。在图1$\bigtriangleup {{ABC}}$中利用几何关系可知：

将式(23)代入式(7)中可得

仍然将标量$||{{\boldsymbol{u}}_{\text{o}}}{\text{||}}_{\text{2}}^{\text{2}}$视为辅助变量，此时综合式(24)中的$N - 1$个方程可以将TDOA伪线性观测方程写为

其中

为了进行短波超视距2D-DOA/TDOA协同定位，需要合并式(17)、式(21)以及式(25)得到矩阵形式为

其中

其满足下面两个二次椭圆等式约束

注释2：类似于文献[8]中的分析，将${\boldsymbol{z}} = {{\boldsymbol{g}}_z}({{\boldsymbol{u}}_{\text{o}}},{\boldsymbol{d}})$代入式(27)中可以得到关于向量${{\boldsymbol{u}}_{\text{o}}}$和d的恒等式，将该恒等式两边分别对向量${{\boldsymbol{u}}_{\text{o}}}$和d求导可得

其中

式(29)对于后续定位方法的推导以及理论性能分析均能起到重要作用。

3.1.2   估计准则

基于式(27)可以定义如下误差向量

其中，$\Delta {\boldsymbol{p}} = {\boldsymbol{p}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}) - {\boldsymbol{p}}({\boldsymbol{z}},{\boldsymbol{d}})$, $\Delta {\boldsymbol{Q}} = {\boldsymbol{Q}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}) - {\boldsymbol{Q}}({\boldsymbol{z}},{\boldsymbol{d}})$。运用一阶误差分析方法可知：

将式(32)和式(33)代入式(31)中可得${\boldsymbol{e}} \approx {{\boldsymbol{C}}_z}{{\boldsymbol{\xi }}_z} + {{\boldsymbol{C}}_d}{{\boldsymbol{\xi }}_d}$，由该式可知，伪线性观测误差向量e渐近服从零均值高斯分布，并且其协方差矩阵为

与经典两步加权最小二乘定位方法类似，阶段1同样忽略了辅助变量$||{{\boldsymbol{u}}_{\text{o}}}{\text{||}}_{\text{2}}^{\text{2}}$与向量${{\boldsymbol{u}}_{\text{o}}}$之间的代数关系，仅考虑式(28)中的第1个等式约束，从而建立如下含有二次椭圆等式约束的优化模型：

由于加权矩阵${({\bf{COV}}({\boldsymbol{e}}))^{ - 1}}$中涉及2D-DOA/TDOA观测误差协方差矩阵${{\boldsymbol{\varSigma }}_z} = {\bf{E}}[{{\boldsymbol{\xi }}_z}{\boldsymbol{\xi }}_z^{\text{T}}]$和电离层虚高观测误差协方差矩阵${{\boldsymbol{\varSigma }}_d} = {\bf{E}}[{{\boldsymbol{\xi }}_d}{\boldsymbol{\xi }}_d^{\text{T}}]$，因此本文新方法需要将这两个矩阵作为先验知识。矩阵${{\boldsymbol{\varSigma }}_z}$中的元素结构与估计信号2D-DOA/TDOA参数时所处理的数据相关性有关，这里参照文献[24]设置该矩阵，其中2D-DOA观测误差协方差矩阵可设为对角矩阵，TDOA观测误差协方差矩阵可设为对角线元素为1、其余元素均为0.5的非对角矩阵乘以某个标量，此外还需考虑2D-DOA观测误差与TDOA观测误差之间的相关性，因此矩阵${{\boldsymbol{\varSigma }}_z}$通常并不具备块状对角结构。矩阵${{\boldsymbol{\varSigma }}_d}$可以通过放置信号源并进行大量的离线观测获得，可将其设为对角矩阵。

3.1.3   求解算法

由于${\nabla ^2}{f_{\text{f}}}{\text{(}}{{\boldsymbol{v}}_{\text{o}}}{\text{)}} = {({\boldsymbol{Q}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{\text{T}}}{({\bf{COV}}({\boldsymbol{e}}))^{ - 1}}{\boldsymbol{Q}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})$是正定矩阵，${{\boldsymbol{\varOmega }}_{\text{1}}}$仅是半正定矩阵，根据文献[29]中的定理2.1可知，式(35)一定存在全局最优解。不妨将其全局最优解记为${{\boldsymbol{\hat v}}_{{\text{o,f}}}}$，利用文献[29]中的定理3.2可得

其中，$\lambda$表示待求解的拉格朗日乘子。式(36)前面两个等式即为KKT (Karush-Kuhn-Tucker)条件。

利用式(36)第1个等式可得

其中

将式(37)代入式(36)第2个等式可知

式(39)是关于$\lambda$的非线性方程，下面将其转化为关于$\lambda$的多项式形式，以使得$\lambda$的求解问题变成一元多项式求根问题。

首先在附录A中证明矩阵${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}{{\boldsymbol{\varOmega }}_1}$可以进行如下对角化：

其中，${{\boldsymbol{\varLambda }}_{\text{w}}} = {\text{diag}}\{ {\kappa _1}\;,\;{\kappa _2}\;,\;{\kappa _3}\;,\;{\kappa _4}\}$表示由特征值构成的对角矩阵，由于${\mathrm{rank}}[{({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}{{\boldsymbol{\varOmega }}_1}] = {\mathrm{rank}}[{{\boldsymbol{\varOmega }}_1}] = 3$，于是有${\kappa _4} = 0$；${{\boldsymbol{H}}_{\text{w}}}$表示由特征向量构成的特征矩阵。利用式(40)可知：

将式(41)代入式(39)中，并且结合式(40)可得

接着定义两个向量为

并将式(43)代入式(42)中可知：

其中，${h_j} = {h_{1,j}}{h_{2,j}}$，其第2个等号是由于${\kappa _4} = 0$。基于式(44)可以得到关于拉格朗日乘子$\lambda$的一元六次多项式方程，如式(45)所示：

其中

利用一元多项式求根算法求解式(45)的根，并将其代入式(37)中即可得到估计值${{\boldsymbol{\hat v}}_{{\text{o,f}}}}$。利用该估计值可以获得辐射源位置向量${{\boldsymbol{u}}_{\text{o}}}$在阶段1的估计结果，不妨将其记为${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$，于是有

其中，${{\boldsymbol{\bar I}}_3} = [{{\boldsymbol{I}}_3}\;\;{{\bf{0}}_3}]$。若将估计值${{\boldsymbol{\hat v}}_{{\text{o,f}}}}$和${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$中的误差向量分别记为$\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}$和$\Delta {{\boldsymbol{u}}_{{\text{o,f}}}}$，则有

注释3：式(45)中的多项式共包含6个根，因此需要排除其中5个假根，可以依次通过以下方式进行排除。首先，将成共轭对的复数根直接排除；然后，将矩阵束$[{{\boldsymbol{\varOmega }}_{\text{1}}}\;,\;{\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})]$的最大广义特征值记为${\kappa _{\max }}$，由式(36)中的第3式可知，$\lambda$的取值范围应为$\lambda \in ( - 1/{\kappa _{\max }}\;,\; + \infty )$；最后，将剩余未被排除的根记为${\{ {\lambda _j}\} _{1 \le j \le J}}$(其中J表示总个数)，并将由每个根得到的辐射源位置向量的估计值记为${\{ {{\boldsymbol{\hat u}}_{{\text{o,f}}}}({\lambda _j})\} _{1 \le j \le J}}$，进而可以利用下面的最大似然估计(Maximum Likelihood Estimation, MLE)准则选定唯一正确的根，如式(53)所示：

众所周知，门限效应对于所有非线性问题都是不可避免的，并且在相同条件下检测发生门限效应时的观测误差阈值通常高于参数估计发生门限效应时的观测误差阈值，这意味着如果标准差增加到使误判发生，那么即使选用正确的根进行定位，其定位精度也将大概率超过CRLB，即参数估计的门限效应也会发生。由于这里选根的最后一步是通过MLE准则进行评判的，因此导致误判发生的观测误差阈值相对较高。第5节的仿真实验结果表明，在大多数情况下仅利用门限值$- 1/{\kappa _{\max }}$就已经能够挑选出唯一正确的根。

3.1.4   理论性能分析

若将向量${{\boldsymbol{\hat v}}_{{\text{o,f}}}}$的估计误差记为$\Delta {{\boldsymbol{v}}_{{\text{o,f}}}} = {{\boldsymbol{\hat v}}_{{\text{o,f}}}} - {\boldsymbol{v}}$，运用约束误差扰动理论可知，误差向量$\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}$是如下优化问题的最优解

式(54)的最优闭式解为

从式(55)中可以看出，误差向量$\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}$渐近服从零均值高斯分布，因此向量${{\boldsymbol{\hat v}}_{{\text{o,f}}}}$是渐近无偏估计值，并且其均方误差矩阵等于

其中，${{\boldsymbol{\varPi }}^ \bot }[ \cdot ]$表示矩阵列补空间上的正交投影矩阵。综合式(52)和式(56)可知，误差向量$\Delta {{\boldsymbol{u}}_{{\text{o,f}}}}$渐近服从零均值高斯分布，并且估计值${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$的均方误差矩阵为

注意到向量${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$虽然是渐近无偏估计值，但却不是渐近统计最优估计值，具体可见如下命题。

命题1：在一阶误差分析理论框架下满足${\text{tr}}[{\bf{MSE}}({{\boldsymbol{\hat u}}_{{\text{o,f}}}})] \ge {\text{tr}}[{\bf{CRLB}}({{\boldsymbol{u}}_{\text{o}}})]$。

证明：首先将式(29)代入式(13)中可得

式(58)中第1个等号和第2个等号分别利用式(34)和式(38)。然后将式(58)代入式(14)中可得

其中，${\boldsymbol{\varPi }}[ \cdot ]$表示矩阵列空间上的正交投影矩阵。另一方面，可以验证

由式(60)可知：

式(61)的证明见附录B。最后，综合式(57)、式(59)以及式(61)可得${\text{tr}}[{\bf{MSE}}({{\boldsymbol{\hat u}}_{{\text{o,f}}}})] \ge {\text{tr}}[{\bf{CRLB}} ({{\boldsymbol{u}}_{\text{o}}})]$。证毕。

估计值${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$不具有渐近统计最优性的原因在于，向量${{\boldsymbol{v}}_{\text{o}}}$中的第4个元素(即辅助变量$||{{\boldsymbol{u}}_{\text{o}}}{\text{||}}_{\text{2}}^{\text{2}}$)与其他元素之间存在相关性，但是该相关性尚未在阶段1中得到利用。另一方面，根据注释1中的讨论可知，向量${{\boldsymbol{v}}_{\text{o}}}$中的第3个元素可以由其前面两个元素表征，为此可以直接舍去该元素，而不会丢失任何有效信息，于是不妨定义一个新的向量${{\boldsymbol{\bar v}}_{\text{o}}} = {[{[{{\boldsymbol{v}}_{\text{o}}}]_1}\;\;{[{{\boldsymbol{v}}_{\text{o}}}]_2}\;\;{[{{\boldsymbol{v}}_{\text{o}}}]_4}]^{\text{T}}}$，其在阶段1的估计值及其估计误差分别为${{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}} = {[{[{{\boldsymbol{\hat v}}_{{\text{o,f}}}}]_1}\;\;{[{{\boldsymbol{\hat v}}_{{\text{o,f}}}}]_2}\;\;{[{{\boldsymbol{\hat v}}_{{\text{o,f}}}}]_4}]^{\text{T}}}$和$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}} = {[{[\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}]_1}\;\;{[\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}]_2}\;\;{[\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}]_4}]^{\text{T}}}$，于是有

其中，${{\boldsymbol{\tilde I}}_3} = \left[\begin{array}{c:c:cc} {{{\boldsymbol{I}}_2}} & {{{\bf{0}}_2}} & {{{\bf{0}}_2}} \\ \hdashline \\[-7pt] {{\bf{0}}_2^{\text{T}}}& 0 & 1 \end{array} \right]$。注意到辅助变量$||{{\boldsymbol{u}}_{\text{o}}}{\text{||}}_{\text{2}}^{\text{2}}$的引入会使误差向量$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$服从某个等式约束，阶段2正是利用该等式约束进一步提高估计精度。

3.2   阶段2的计算原理与方法

阶段2需要首先推导误差向量$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$满足的等式约束，并由此获得向量$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$的估计值，进而得到向量${{\boldsymbol{\bar u}}_{\text{o}}}$的估计值，最终获得精度更高的短波辐射源定位结果。根据向量${{\boldsymbol{\bar v}}_{\text{o}}}$中的第3个元素的定义可知

将式(63)第2个等号右边第1项在点${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$处进行一阶泰勒级数展开可得

其中

基于式(54)中的等式约束可知，误差向量$\Delta {{\boldsymbol{v}}_{{\text{o,f}}}}$与$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$之间满足式(66)关系式：

将式(66)代入式(64)中可得

式(67)即为误差向量$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$所满足的等式约束。基于误差向量$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$的高斯分布特性可以建立估计误差向量$\Delta {{\boldsymbol{\bar v}}_{{\text{o,f}}}}$的等式约束优化模型，如式(68)所示：

最后利用拉格朗日乘子技术可得式(68)的最优闭式解为

基于式(69)可以获得向量${{\boldsymbol{\bar u}}_{\text{o}}}$的估计值，如式(70)所示：

其中，${{\boldsymbol{\bar I}}_2} = [{{\boldsymbol{I}}_2}\;\;{{\bf{0}}_2}]$。根据注释1中的讨论可以将辐射源位置向量${{\boldsymbol{u}}_{\text{o}}}$在阶段2中的估计值表示为

式(71)中第3个元素的符号与阶段1的估计结果$< {{\boldsymbol{\hat u}}_{{\text{o,f}}}}{ > _3}$的符号保持一致。

3.3   协同定位新方法的计算步骤与讨论

综合前面的讨论，表1总结了本文新方法的计算步骤，并且给出了每个步骤所需的乘法次数。

表  1  新方法的计算步骤及其所需的乘法次数

Table  1.  Calculation steps of the new method and the number of multiplications required

序号	计算内容	乘法次数
步骤1	利用式(18)计算矩阵${{\boldsymbol{Q}}_\theta }$和向量${{\boldsymbol{p}}_\theta }$；利用式(22)计算矩阵${{\boldsymbol{Q}}_\beta }$和向量${{\boldsymbol{p}}_\beta }$； 利用式(26)计算矩阵${{\boldsymbol{Q}}_t}$和向量${{\boldsymbol{p}}_t}$	$27N - 8$
步骤2	利用式(27)构造矩阵${\boldsymbol{Q}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})$和向量${\boldsymbol{p}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})$	无实质乘法计算
步骤3	利用式(30)计算矩阵${{\boldsymbol{C}}_z}$和${{\boldsymbol{C}}_d}$	$48{N^2} - 4N - 13$
步骤4	利用式(34)计算协方差矩阵${\bf{COV}}({\boldsymbol{e}})$及其逆矩阵${({\bf{COV}}({\boldsymbol{e}}))^{ - 1}}$	$\begin{gathered} O[{(3N - 1)^3}] + 66{N^3} \\[-3pt] - 61{N^2} + 19N - 2 \end{gathered}$
步骤5	利用式(38)计算矩阵${\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})$及其逆矩阵${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}$	$\begin{gathered} O(64) + 36{N^2} \\[-3pt] + 24N - 12 \end{gathered}$
步骤6	对矩阵${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}{{\boldsymbol{\varOmega }}_1}$进行特征值分解得到特征值对角矩阵${{\boldsymbol{\varLambda }}_{\text{w}}}$和特征矩阵${{\boldsymbol{H}}_{\text{w}}}$	$O(64) + 4$
步骤7	利用式(43)计算向量${{\boldsymbol{h}}_1}$和${{\boldsymbol{h}}_2}$	$\begin{gathered} O(64) + 9{N^2} \\ + 6N + 45 \\ \end{gathered}$
步骤8	利用式(46)—式(50)计算多项式系数${\{ {\alpha _j}\} _{0 \le j \le 6}}$	$82$
步骤9	求解式(45)中的一元六次多项式的根，并利用式(53)选择正确的根	$\begin{gathered} J \cdot O[{(3N - 1)^3}] + O(6) \\[-3pt] + J(12{N^3} + 2{N^2} - 2N) \end{gathered}$
步骤10	利用式(37)计算向量${{\boldsymbol{\hat v}}_{{\text{o,f}}}}$，并进而构造向量${{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}}$和${{\boldsymbol{\hat u}}_{{\text{o,f}}}}$	$O(64) + 16$
步骤11	利用式(62)计算均方误差矩阵${\bf{MSE}}({{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}})$及其逆矩阵${({\bf{MSE}}({{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}}))^{ - 1}}$	$O(64) + O(27) + 145$
步骤12	利用式(65)计算标量$\hat \mu$和向量${\boldsymbol{\hat \eta }}$	6
步骤13	利用式(66)构造矩阵${\boldsymbol{\bar J}}({{\boldsymbol{\hat u}}_{{\text{o,f}}}})$	5
步骤14	利用式(69)计算向量$\Delta {{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}}$	$28$
步骤15	利用式(70)计算向量${{{\hat{ \bar {\boldsymbol{u}}}}}_{{\text{o,s}}}}$	无实质乘法计算
步骤16	利用式(71)计算最终定位结果${{\boldsymbol{\hat u}}_{{\text{o,s}}}}$	3

下载: 导出CSV 

| 显示表格

关于上述定位方法，下面给出两点注释：

注释4：为确保本文新方法的计算能够顺利进行，需要式(38)中的矩阵${\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})$是可逆的，这要求式(27)中的$(3N - 1) \times 4$矩阵${\boldsymbol{Q}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}})$具有列满秩性，由此可知观测站个数N应满足$3N - 1 \ge 4 \Rightarrow N \ge 2$，所以新方法仅需两个观测站即可实现定位。该结论并不难以理解，事实上在只有2D-DOA观测量的条件下也仅需两个观测站即可实现交汇定位，而新方法在融合TDOA观测量的条件下自然不会增加观测量个数N的下界。另一方面，由于文献[8]中的场景假设观测站仅获得2D-DOA观测量和TDOA观测量中的一种，采用类似的分析可知其至少需要5个观测站才能实现定位。此外，两种方法的定位精度都会随着观测站个数的增加而得到提升。

注释5：新方法无需迭代运算，其中阶段1所需的乘法次数为

阶段2所需的乘法次数为$O(64) + O(27) + 187$，因此总的乘法次数为

显然，阶段1的乘法次数是${N^3}$阶，而阶段2的乘法次数与N无关，是常数项。

4.   理论性能分析

本节运用约束误差扰动理论和正交投影矩阵的数学性质推导定位结果${{\boldsymbol{\hat u}}_{{\text{o,s}}}}$的统计特性，主要推导其估计均方误差，并将所得到的结果与CRLB进行比较，以证明新方法的渐近统计最优性。虽然本文新方法与文献[8]中的定位方法都需要两个计算阶段，但两种方法的计算框架有着显著差异，因此两者理论性能的推导思路也截然不同。此外，本节还新增分析短波辐射源高度信息误差的影响，并提出相应定位均方误差的理论表达式。

4.1   估计均方误差

首先将式(68)中的等式约束代入式(69)中可得

其中，

由式(74)可知，估计值$\Delta {{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}}$渐近服从零均值高斯分布，并且其均方误差矩阵等于

结合式(70)和式(75)可知，估计值${{{\hat{ \bar {\boldsymbol{u}}}}}_{{\text{o,s}}}}$渐近服从零均值高斯分布，并且其均方误差矩阵为

因此最终定位结果${{\boldsymbol{\hat u}}_{{\text{o,s}}}}$也渐近服从零均值高斯分布，并且其均方误差矩阵等于

4.2   渐近统计最优性分析

下面证明定位结果${{\boldsymbol{\hat u}}_{{\text{o,s}}}}$具有渐近统计最优性，为此需要证明估计值${{{\hat {\bar {\boldsymbol{u}}}}}_{{\text{o,s}}}}$也具有渐近统计最优性。为了完成证明，下面首先推导均方误差矩阵${\bf{MSE}}({{{\hat{ \bar {\boldsymbol{u}}}}}_{{\text{o,s}}}})$的另一种表达式，具体可见命题2。

命题2：均方误差矩阵${\bf{MSE}}({{{\hat{ \bar {\boldsymbol{u}}}}}_{{\text{o,s}}}})$可以表示为

证明：首先利用正交投影矩阵的定义可以将式(76)重新表示为

考虑矩阵${({\bf{MSE}}({{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}}))^{ - 1/2}} \dfrac{{\partial {{{\boldsymbol{\bar v}}}_{\text{o}}}}}{{\partial {\boldsymbol{\bar u}}_{\text{o}}^{\text{T}}}} \in {\mathbb{R}^{3 \times 2}}$和向量${({\bf{MSE}}({{{\hat {\bar {\boldsymbol{v}}}}}_{{\text{o,f}}}}))^{1/2}}\,{({\boldsymbol{\bar J}}({{\boldsymbol{u}}_{\text{o}}}))^{\text{T}}}{\boldsymbol{\eta }} \in {\mathbb{R}^{3 \times 1}}$，他们的列数之和等于3，并且满足

于是有

此时利用正交投影矩阵的性质可知：

接着将式(82)代入式(79)中可得

不难验证${{\boldsymbol{\bar I}}_2} \dfrac{{\partial {{{\boldsymbol{\bar v}}}_{\text{o}}}}}{{\partial {\boldsymbol{\bar u}}_{\text{o}}^{\text{T}}}} = {{\boldsymbol{I}}_2}$，最后将该式代入式(83)中可知式(78)成立。证毕。

命题3：在一阶误差分析理论框架下满足${\bf{MSE}}({{{\hat {\bar {\boldsymbol{u}}}}}_{{\text{o,s}}}}) = {\bf{CRLB}}({{\boldsymbol{\bar u}}_{\text{o}}})$。

证明：首先将式(56)和式(62)代入式(78)中，并且结合关系式$\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{I}}_2}} \\ {2{e^2}{\boldsymbol{\bar u}}_{\text{o}}^{\text{T}}} \end{array}} \right] = {{\boldsymbol{\tilde I}}_3}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{J}}({{\boldsymbol{u}}_{\text{o}}})} \\ {2{e^2}{\boldsymbol{\bar u}}_{\text{o}}^{\text{T}}} \end{array}} \right]$可知

然后将式(58)代入式(15)中可得

式(85)中第2个等号利用了关系式$\dfrac{{\partial {{\boldsymbol{v}}_{\text{o}}}}}{{\partial {\boldsymbol{u}}_{\text{o}}^{\text{T}}}}{\boldsymbol{J}}({{\boldsymbol{u}}_{\text{o}}}) = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{J}}({{\boldsymbol{u}}_{\text{o}}})} \\ {2{e^2}{\boldsymbol{\bar u}}_{\text{o}}^{\text{T}}} \end{array}} \right]$。结合式(84)和式(85)可知：

利用矩阵和的求逆公式可得

其中，${\boldsymbol{\tilde W}}({\boldsymbol{z}},{\boldsymbol{d}}) = {\boldsymbol{\tilde I}}_3^{\text{T}}{[{{\boldsymbol{\tilde I}}_3}{({\boldsymbol{W}}({\boldsymbol{z}},{\boldsymbol{d}}))^{ - 1}}{\boldsymbol{\tilde I}}_3^{\text{T}}]^{ - 1}}{{\boldsymbol{\tilde I}}_3}$。接着将式(87)代入式(86)中可知

附录C中证明了矩阵等式式(89)：

综合式(88)和式(89)可知${({\bf{MSE}}({{{\hat {\bar {\boldsymbol{u}}}}}_{{\text{o,s}}}}))^{ - 1}} = {({\bf{CRLB}}({{\boldsymbol{\bar u}}_{\text{o}}}))^{ - 1}}$，由此可得${\bf{MSE}}({{{\hat {\bar {\boldsymbol{u}}}}}_{{\text{o,s}}}}) = {\bf{CRLB}}({{\boldsymbol{\bar u}}_{\text{o}}})$。证毕。

最后结合命题3、式(77)以及式(16)可知

式(90)证明了定位结果${{\boldsymbol{\hat u}}_{{\text{o,s}}}}$的渐近统计最优性。

4.3   短波辐射源高度信息误差的影响

本文定位新方法假设短波辐射源位于地球表面(即高度为零)，如果短波辐射源高度并不为零，但是先验已知，此时基于文献[30]中的讨论同样可以建立类似于式(6)的二次椭圆等式约束，因此对新方法的计算过程并无本质影响。然而，如果短波辐射源高度无法先验已知，则被忽略的高度必然会对新方法的定位精度产生负面影响，本小节仍然运用约束误差扰动理论对此影响进行定量分析。

假设短波辐射源高度为${h_{\text{e}}}$，此时式(1)中的位置向量应修正为^[30]

其中，${N_{\text{e}}} = {r_{\text{e}}}/\sqrt {1 - {e^2}{{(\sin ({\rho _{\text{o}}}))}^2}}$。当高度${h_{\text{e}}}$较小时，利用一阶近似可得

其中，${\varphi _{\text{o}}} = 2{N_{\text{e}}} = 2{r_{\text{e}}}/\sqrt {1 - {e^2}{{(\sin ({\rho _{\text{o}}}))}^2}}$。

由于新方法所求解的等式约束优化问题可以表示为

若在高度${h_{\text{e}}}$非零的情况下将向量${{\boldsymbol{v}}_{\text{o}}}$和${{\boldsymbol{u}}_{\text{o}}}$的估计值分别记为${{\boldsymbol{\hat v}}_{{\text{o,h}}}}$和${{\boldsymbol{\hat u}}_{{\text{o,h}}}} = {{\boldsymbol{\bar I}}_3}{{\boldsymbol{\hat v}}_{{\text{o,h}}}}$，相应的估计误差分别记为$\Delta {{\boldsymbol{v}}_{{\text{o,h}}}} = {{\boldsymbol{\hat v}}_{{\text{o,h}}}} - {{\boldsymbol{v}}_{\text{o}}}$和$\Delta {{\boldsymbol{u}}_{{\text{o,h}}}} = {{\boldsymbol{\hat u}}_{{\text{o,h}}}} - {{\boldsymbol{u}}_{\text{o}}} = {{\boldsymbol{\bar I}}_3}\Delta {{\boldsymbol{v}}_{{\text{o,h}}}}$，则误差向量$\Delta {{\boldsymbol{v}}_{{\text{o,h}}}}$应是式(94)等式约束优化问题的最优解

若令${\boldsymbol{\bar \varOmega }} = \left[ {\begin{array}{*{20}{c}} {{{({{\boldsymbol{\varOmega }}_{\text{1}}}{{\boldsymbol{v}}_{\text{o}}})}^{\text{T}}}} \\ {{{({{\boldsymbol{\varOmega }}_{\text{2}}}{{\boldsymbol{v}}_{\text{o}}} + {{\boldsymbol{\gamma }}_2}{\text{/}}2)}^{\text{T}}}} \end{array}} \right]$，利用拉格朗日乘子法可以得到误差向量$\Delta {{\boldsymbol{v}}_{{\text{o,h}}}}$的表达式为

基于式(95)可得

因此，当高度${h_{\text{e}}}$非零时，新方法是有偏估计器，即存在定位偏置。另一方面，估计值${{\boldsymbol{\hat u}}_{{\text{o,h}}}}$的均方误差矩阵等于

附录D和附录E中分别证明了式(98)，式(99)两个矩阵等式

最后将式(98)和式(99)代入式(97)中可得

注释6：首先根据式(100)可知${\text{tr}}[{\bf{MSE}}({{\boldsymbol{\hat u}}_{{\text{o,h}}}})] \ge {\text{tr}}[{\bf{CRLB}}({{\boldsymbol{u}}_{\text{o}}})]$，即被忽略的辐射源高度${h_{\text{e}}}$必然会增加定位均方误差，从而降低定位精度。此外，为了使二次椭圆等式约束式(6)产生增益(即${\text{tr}}[{\bf{MSE}}({{\boldsymbol{\hat u}}_{{\text{o,h}}}})] \le {\text{tr}}[{{\bf{CRLB}}_{{\text{un}}}}({{\boldsymbol{u}}_{\text{o}}})]$)，由式(100)可知高度${h_{\text{e}}}$应满足${h_{\text{e}}} \le 2\sqrt {{{({{\boldsymbol{\varOmega }}_{\text{0}}}{{\boldsymbol{u}}_{\text{o}}})}^{\text{T}}}{{\bf{CRLB}}_{{\text{un}}}}({{\boldsymbol{u}}_{\text{o}}}){{\boldsymbol{\varOmega }}_{\text{0}}}{{\boldsymbol{u}}_{\text{o}}}} /{\varphi _{\text{o}}}$。该式表明，在没有二次椭圆等式约束式(6)的条件下，定位精度越高(即${\text{tr}}[{\bf{CRL}}{{\bf{B}}_{{\text{un}}}}({{\boldsymbol{u}}_{\text{o}}})]$较小)，所允许的高度${h_{\text{e}}}$的范围就越小。

注释7：本节仅从理论上分析了短波辐射源高度信息误差的影响，在实际定位中还应考虑如何有效抑制该误差，这里提供两种策略：第1种是将辐射源高度信息也作为未知参数进行处理，并且为了融入辐射源高度先验信息，可以通过贝叶斯估计框架完成定位；第2种是将式(28)中的第1个二次椭圆等式约束松弛为不等式约束形式${\boldsymbol{v}}_{\text{o}}^{\text{T}}{{\boldsymbol{\varOmega }}_{\text{1}}}{{\boldsymbol{v}}_{\text{o}}} \le {({r_{\text{e}}} + \Delta r)^2}$，其中$\Delta r$表示半径增量，然后再通过有效集方法、凸优化方法等进行寻优计算。

5.   仿真分析

在未特殊说明的情况下本节的基础实验参数为：短波辐射源经度和纬度分别为128.14°和32.20°，其在ECEF坐标系下的位置向量为${{\boldsymbol{u}}_{\text{o}}} = {[ - 3336.4\;\;4248.9\;\;3379.2]^{\text{T}}}$ km，现有5个观测站的阵列天线接收该短波信号，并且能联合估计信号的2D-DOA/TDOA两域参数，他们的经纬度与信号到达观测站对应的电离层虚高数值见表2，将2D-DOA观测误差协方差矩阵设为${{\boldsymbol{\varSigma }}_1} = \sigma _1^{\text{2}}{{\boldsymbol{I}}_{10}}$，TDOA观测误差协方差矩阵设为${{\boldsymbol{\varSigma }}_2} = \sigma _2^{\text{2}}({{\boldsymbol{I}}_4} + {{\bf{1}}_{4 \times 4}})/2$，2D-DOA观测误差与TDOA观测误差之间的协方差矩阵设为${{\boldsymbol{\varSigma }}_{1,2}} = {c_{{\text{DT}}}}{\sigma _1}{\sigma _2}\left( {{{\bf{1}}_2} \otimes \left[ {\begin{array}{*{20}{c}} {{\bf{1}}_4^{\text{T}}/2} \\ {{{\boldsymbol{I}}_4}/2} \end{array}} \right]} \right)$，其中${\sigma _1}$和${\sigma _2}$分别表示2D-DOA观测误差标准差和TDOA观测误差标准差，${c_{{\text{DT}}}}$表示2D-DOA/TDOA两种观测误差的相关系数，将电离层虚高观测误差协方差矩阵设为${{\boldsymbol{\varSigma }}_d} = \sigma _d^2{{\boldsymbol{I}}_5}$，其中${\sigma _d}$表示电离层虚高观测误差标准差。

表  2  观测站的经纬度与电离层虚高数值

Table  2.  Longitude, latitude and ionospheric virtual height of shortwave observer

观测站	经度(°)	纬度(°)	电离层虚高(km)
1	123.04	40.94	273
2	115.66	40.27	316
3	114.61	34.98	358
4	113.47	30.09	347
5	115.75	26.22	255

下载: 导出CSV 

| 显示表格

5.1   验证新方法的有效性

将2D-DOA, TDOA以及电离层虚高观测误差标准差分别设为${\sigma _1} = 0.5^\circ$, ${\sigma _2} = 5{\text{/}}{\mathrm{c}} \;{\mathrm{s}}$以及${\sigma _d} = {\text{5 km}}$，2D-DOA/TDOA两种观测误差的相关系数设为${c_{{\text{DT}}}} = 0.1$，并且进行5000次蒙特卡罗实验，图2给出了新方法在ECEF坐标系下的X-Y平面中的定位结果散布图与定位误差椭圆曲线。

图  2  在ECEF坐标系下的X-Y平面中的定位结果散布图与定位误差椭圆曲线

Figure  2.  Scatter plot of positioning results and elliptical curve of positioning errors in the X-Y plane in the ECEF coordinate system

下载: 全尺寸图片 幻灯片

从图2中可以看出，本文新方法的定位误差椭圆的形状与定位结果散点图的形状相匹配，并且概率值越小对应的定位误差椭圆面积越小，概率值越大对应的定位误差椭圆面积越大，由此验证了新方法的有效性。

下面给出选根的实验结果。表3给出了在单次实验中一元六次多项式的6个根以及他们对应的MLE目标函数值(见式(53))，注意到由最大广义特征值得到的门限值为$- 1/{\kappa _{\max }} = - {\text{4}}{\text{.02}} \times {\text{1}}{{\text{0}}^{ - {\text{5}}}}$。从表3中可以看出，利用门限值$- 1/{\kappa _{\max }}$可以选出唯一正确的根(对应序号6)，而利用MLE目标函数值也能将序号6对应的根挑选出来。事实上，通过大量实验发现在大多数情况下仅通过门限值$- 1/{\kappa _{\max }}$就已经能把唯一正确的根挑选出来，该结论可通过图3得到验证，图3中给出了5000次蒙特卡罗实验中正确的根、与正确的根最邻近的假根以及门限值$- 1/{\kappa _{\max }}$。在少数出现多个根大于门限值$- 1/{\kappa _{\max }}$的情形下，利用MLE准则也总能选出唯一的根。当然，随着各类观测误差标准差的不断增大，终究会出现误判的情形，但由于最后一步是通过MLE准则进行评判，因此导致误判发生的观测误差阈值是相对较高的，并且高于后续定位发生门限效应时所对应的观测误差阈值。

表  3  多项式方程的根及其对应的MLE目标函数值

Table  3.  Roots of polynomial equation and their corresponding MLE objective function values

序号	一元六次多项式方程的根	与门限值–1/${\kappa _{\max }}$= –4.02×10–5的关系	由式(53)给出的MLE目标函数值
1	–8.12×10–2	小于	3.34×106
2	–7.98×10–2	小于	3.22×106
3	–1.81×10–2	小于	7.56×105
4	–1.73×10–2	小于	6.90×105
5	–8.08×10–5	小于	6.60×103
6	2.89×10–7	唯一大于	1.23×101 唯一的最小值

下载: 导出CSV 

| 显示表格

图  3  正确的根、与正确的根最邻近的假根以及门限值$- 1/{\kappa _{\max }}$

Figure  3.  Correct root, the closest false root to the correct root and threshold value $- 1/{\kappa _{\max }}$

下载: 全尺寸图片 幻灯片

5.2   验证协同定位的优势

本节将新方法与文献[9]中的短波2D-DOA定位方法和文献[23]中的短波TDOA定位方法进行比较，以验证新方法能获得较高的协同增益。文献[9]中的定位方法属迭代型方法，这里将其迭代初始值设为真实值，文献[23]中的定位方法属网格选择方法，可通过调整参数使其性能达到渐近统计最优，将2D-DOA/TDOA两种观测误差的相关系数设为${c_{{\text{DT}}}} = 0.1$。首先将TDOA和电离层虚高观测误差标准差分别设为${\sigma _2} = 3{\text{/}}{\mathrm{c}} \;{\mathrm{s}}$和${\sigma _d} = 3{\text{ km}}$，2D-DOA标准差设为${\sigma _1} = 0.{\text{05}} \times {\delta _1}^\circ$，图4给出了定位均方根误差随着参数${\delta _1}$的变化曲线；接着将2D-DOA和电离层虚高观测误差标准差分别设为${\sigma _1} = 0.8^\circ$和${\sigma _d} = 3{\text{ km}}$，TDOA标准差设为${\sigma _2} = 0.3 \times {\delta _2}{\text{/}}{\mathrm{c}}{\text{ s}}$，图5给出了定位均方根误差随着参数${\delta _2}$的变化曲线；最后将2D-DOA和TDOA标准差分别设为${\sigma _1} = 0.8^\circ$和${\sigma _d} = 3{\text{ km}}$，电离层虚高观测误差标准差设为${\sigma _d} = 0.{\text{3}} \times {\delta _3}{\text{ km}}$，图6给出了定位均方根误差随着参数${\delta _3}$的变化曲线。

图  4  定位均方根误差随着参数δ₁的变化曲线

Figure  4.  Localization RMSE as a function of parameter δ₁

下载: 全尺寸图片 幻灯片

图  5  定位均方根误差随着参数δ₂的变化曲线

Figure  5.  Localization RMSE as a function of parameter δ₂

下载: 全尺寸图片 幻灯片

图  6  定位均方根误差随着参数δ₃的变化曲线

Figure  6.  Localization RMSE as a function of parameter δ₃

下载: 全尺寸图片 幻灯片

从图4至图6中可以看出：(1)本文新方法的定位均方根误差总能渐近逼近CRLB，从而验证了第4节理论性能分析的有效性；(2)文献[9]中的短波2D-DOA定位方法和文献[23]中的短波TDOA定位方法均具有渐近统计最优性，但这两种方法的定位均方根误差都明显高于新方法的定位均方根误差，从而说明了新方法能获得较高的协同增益。

5.3   与其他短波辐射源2D-DOA/TDOA协同定位方法的性能比较

由于文献[8]中的短波辐射源2D-DOA/TDOA协同定位方法可直接应用于本文场景，因此本节将新方法与文献[8]中的定位方法进行比较。注意到文献[8]中的定位方法假设每个观测站仅获得2D-DOA/TDOA中的一种观测量，因此该方法认为2D-DOA/TDOA两种观测误差统计独立。然而，本文假设每个观测站同时获得2D-DOA/TDOA两域参数，此时这两种观测误差统计相关，但文献[8]中的定位方法并未考虑此相关性，因此其定位精度无法达到本文场景下的CRLB，下面的实验验证该结论。将2D-DOA, TDOA以及电离层虚高观测误差标准差分别设为${\sigma _1} = 0.5^\circ$, ${\sigma _2} = 5{\text{/}}{\mathrm{c}}{\text{ s}}$以及${\sigma _d} = {\text{5 km}}$，改变2D-DOA/TDOA两种观测误差的相关系数${c_{{\text{DT}}}}$，使其从0增加至0.95，图7给出了两种方法的定位均方根误差随着相关系数${c_{{\text{DT}}}}$的变化曲线。

图  7  定位均方根误差随着相关系数c_DT的变化曲线

Figure  7.  Localization RMSE as a function of correlation coefficient c_DT

下载: 全尺寸图片 幻灯片

从图7可以看出：(1)随着2D-DOA/TDOA两种观测误差相关系数${c_{{\text{DT}}}}$的增加，两种方法的定位精度都得到了提升，这是因为相关系数${c_{{\text{DT}}}}$的增加预示着观测误差总体不确定性的降低，从而使定位误差逐渐减小；(2)在本文场景下，新方法的精度高于文献[8]中的定位方法的精度，并且两者的性能差异随着相关系数${c_{{\text{DT}}}}$的增加而增大，从而验证了前面的结论；(3)即使当相关系数${c_{{\text{DT}}}} = 0$时，文献[8]中的定位方法的估计均方根误差仍然高出了CRLB达0.262 km，因为该方法不仅假设2D-DOA/TDOA两种观测误差统计独立，还假设测DOA观测站与测TDOA观测站各自对应的电离层虚高观测误差也相互统计独立，但是在本文场景中，2D-DOA/TDOA两域参数对应相同的观测站，所以他们各自对应的电离层虚高观测误差完全相同(并不统计独立)，而文献[8]中的定位方法未考虑此信息，从而损失了0.262 km的定位精度。

总体而言，两种方法在各自考虑的场景下均具有渐近统计最优性，但如果将文献[8]中的定位方法直接应用于本文场景，由于一些观测误差相关性的丢失，其定位性能无法达到本文场景下的CRLB。

5.4   验证短波辐射源高度信息误差的影响

将2D-DOA, TDOA以及电离层虚高观测误差标准差分别设为${\sigma _1} = 0.3^\circ$, ${\sigma _2} = 3{\text{/}}{\mathrm{c}}{\text{ s}}$以及${\sigma _d} = {\text{3 km}}$，2D-DOA/TDOA两种观测误差的相关系数设为${c_{{\text{DT}}}} = 0.1$，辐射源高度${h_{\text{e}}}$从5 km变化至100 km。假设定位方法忽略了此高度，图8给出了本文新方法的定位偏置随着被忽略的辐射源高度${h_{\text{e}}}$的变化曲线，图8中同时包含了ECEF坐标系下X轴、Y轴以及Z轴方向上的定位偏置，并且将蒙特卡罗实验值与式(96)给出的理论值进行比较。图9给出了新方法的定位均方根误差随着被忽略的辐射源高度${h_{\text{e}}}$的变化曲线，图9中给出了没有等式约束条件下的CRLB，并且将蒙特卡罗实验值与式(100)给出的理论值进行比较。

图  8  定位偏置随着高度h_e的变化曲线

Figure  8.  Localization bias as a function of height h_e

下载: 全尺寸图片 幻灯片

图  9  定位均方根误差随着高度h_e的变化曲线

Figure  9.  Localization RMSE as a function of height h_e

下载: 全尺寸图片 幻灯片

从图8和图9中可以看出：(1)当辐射源高度${h_{\text{e}}}$被忽略时，将导致定位偏置，并且随着高度${h_{\text{e}}}$的增加，其定位偏置逐渐增大；(2)定位均方根误差也随着辐射源高度${h_{\text{e}}}$的增加而增大，当${h_{\text{e}}}$超过门限值后，定位均方根误差将高于没有二次椭圆等式约束条件下的CRLB，此时意味着二次椭圆等式约束将不再起作用；(3)定位偏置与式(96)给出的理论值基本吻合，定位均方根误差与式(100)给出的理论值基本吻合，从而验证了4.3节理论性能分析的有效性。

6.   结语

本文基于单跳电离层虚高模型，在观测站同时获得2D-DOA/TDOA两域参数的场景下提出一种协同这两种观测量的定位新方法。新方法可以有效避免迭代运算，并利用两个计算阶段获得目标位置向量闭式解。此外，该文利用约束误差扰动理论和正交投影矩阵的数学性质对新方法的估计性能进行理论分析，证明其具有渐近统计最优性，同时定量分析短波辐射源高度信息误差对定位精度产生的影响。仿真实验结果表明，相比于短波2D-DOA定位方法和短波TDOA定位方法，新方法能获得较高的协同增益；相比于文献[8]中的短波辐射源2D-DOA/TDOA协同定位方法，新方法在两类观测误差具有统计相关性的条件下具有更高定位精度，并且随着相关性的增加其优势更加显著。最后需要指出的是，在电离层影响下短波信号传播易发生多径效应，并且电离层参数具有较强的时变性，如何将本文的新方法应用于多径传播、电离层参数时变等复杂环境和边界场景值得后续进一步深入研究。

附录A　证明矩阵${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}{{\boldsymbol{\varOmega }}_1}$可以进行对角化

附录A首先证明如下引理。

引理1：假设${{\boldsymbol{A}}_1} \in {\mathbb{R}^{n \times n}}$和${{\boldsymbol{A}}_2} \in {\mathbb{R}^{n \times n}}$均为正定矩阵，则矩阵${{\boldsymbol{A}}_1}{{\boldsymbol{A}}_2} \in {\mathbb{R}^{n \times n}}$可以进行对角化。

证明：首先有${{\boldsymbol{A}}_1}{{\boldsymbol{A}}_2} = {\boldsymbol{A}}_2^{ - 1/2}({\boldsymbol{A}}_2^{1/2}{{\boldsymbol{A}}_1}{\boldsymbol{A}}_2^{1/2}){\boldsymbol{A}}_2^{1/2}$，由于${\boldsymbol{A}}_2^{1/2}$是正定矩阵，${\boldsymbol{A}}_2^{1/2}{{\boldsymbol{A}}_1}{\boldsymbol{A}}_2^{1/2}$也是正定矩阵，因此该矩阵可以进行对角化，即有${\boldsymbol{A}}_2^{1/2}{{\boldsymbol{A}}_1}{\boldsymbol{A}}_2^{1/2} = {{\boldsymbol{H}}_{\text{a}}}{{\boldsymbol{\varLambda }}_{\text{a}}}{\boldsymbol{H}}_{\text{a}}^{ - 1}$，其中，${{\boldsymbol{\varLambda }}_{\text{a}}}$表示由特征值构成的对角矩阵；${{\boldsymbol{H}}_{\text{a}}}$表示由特征向量构成的特征矩阵。至此可得

于是矩阵${{\boldsymbol{A}}_1}{{\boldsymbol{A}}_2}$也可以进行对角化。证毕。

下面基于引理1证明矩阵${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}{{\boldsymbol{\varOmega }}_1}$可以进行对角化。由于${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}$是正定矩阵，可将其分块表示为

其中左上角子矩阵${{\boldsymbol{\hat W}}_1}$是正定矩阵。利用矩阵${{\boldsymbol{\varOmega }}_{\text{1}}}$的定义可得

由于${{\boldsymbol{\hat W}}_1}$和${{\boldsymbol{\varOmega }}_{\text{0}}}$均为正定矩阵，根据引理1可知，${{\boldsymbol{\hat W}}_1}{{\boldsymbol{\varOmega }}_{\text{0}}}$可以进行对角化，即有${{\boldsymbol{\hat W}}_1}{{\boldsymbol{\varOmega }}_{\text{0}}} = {{\boldsymbol{\bar H}}_{\text{w}}}{{\boldsymbol{\bar \varLambda }}_{\text{w}}}{\boldsymbol{\bar H}}_{\text{w}}^{ - 1}$，其中，${{\boldsymbol{\bar \varLambda }}_{\text{w}}}$表示由特征值构成的对角矩阵；${{\boldsymbol{\bar H}}_{\text{w}}}$表示由特征向量构成的特征矩阵。将该式代入式(A-3)中可得

其中，${{\boldsymbol{\bar h}}_{\text{w}}} = {\boldsymbol{\bar \varLambda }}_{\text{w}}^{ - 1}{\boldsymbol{\bar H}}_{\text{w}}^{\text{T}}{{\boldsymbol{\varOmega }}_{\text{0}}}{{\boldsymbol{\hat w}}_{\text{2}}}$。由于

将式(A-5)代入式(A-4)中可得

式中，${{\boldsymbol{\varLambda }}_{\text{w}}} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\bar \varLambda }}}_{\text{w}}}}&{{{\bf{0}}_3}} \\ {{\bf{0}}_3^{\text{T}}}&0 \end{array}} \right]$表示由特征值构成的对角矩阵；${{\boldsymbol{H}}_{\text{w}}} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\bar H}}}_{\text{w}}}}&{{{\bf{0}}_3}} \\ {{\boldsymbol{\bar h}}_{\text{w}}^{\text{T}}}&1 \end{array}} \right]$表示由特征向量构成的特征矩阵。综上可知，矩阵${({\boldsymbol{W}}({\boldsymbol{\hat z}},{\boldsymbol{\hat d}}))^{ - 1}}{{\boldsymbol{\varOmega }}_1}$可以进行对角化。

附录B　证明式(61)

附录B首先证明引理2。

引理2：若矩阵${\boldsymbol{A}} \in {\mathbb{R}^{m \times n}}\;(m > n)$列满秩，非零向量${\boldsymbol{a}} \in {\mathbb{R}^{m \times 1}}$满足${\text{range}}[{\boldsymbol{a}}] \bot {\text{range}}[{\boldsymbol{A}}]$，则有${\boldsymbol{\varPi }}[\begin{array}{c:ccc} {\boldsymbol{A}}& {\boldsymbol{a}} \end{array}] \ge {\boldsymbol{\varPi }}[{\boldsymbol{A}}]$。

证明：根据正交投影矩阵的定义可得

由于式(B-1)中最后一个等号右边第2项为半正定矩阵，于是有${\boldsymbol{\varPi }}\left[\begin{array}{c:ccc} {}&{}\\[-8pt] {\boldsymbol{A}}& {\boldsymbol{a}}\\ {}&{}\\[-8pt]\end{array}\right] \ge {\boldsymbol{\varPi }}[{\boldsymbol{A}}]$。证毕。

下面基于引理2证明式(61)。首先由式(60)可得

于是一定存在非零向量${{\boldsymbol{w}}_v} \in {\mathbb{R}^{4 \times 1}}$满足

然后结合式(B-3)和式(60)可知

此时根据正交投影矩阵的性质可知

最后利用引理2中的结论可知式(61)成立。

附录C　证明式(89)

附录C首先证明引理3。

引理3：假设${{\boldsymbol{A}}_1} \in {\mathbb{R}^{n \times n}}$为正定矩阵，${{\boldsymbol{A}}_2} \in {\mathbb{R}^{n \times (n - 1)}}$为列满秩矩阵，矩阵${{\boldsymbol{A}}_3} \in {\mathbb{R}^{n \times m}} \;(m < n)$与向量${\boldsymbol{a}} \in {\mathbb{R}^{n \times 1}}$满足${{\boldsymbol{a}}^{\text{T}}}{{\boldsymbol{A}}_3} = {\bf{0}}_m^{\text{T}}$，则有矩阵等式为

证明：首先利用关系式${{\boldsymbol{a}}^{\text{T}}}{{\boldsymbol{A}}_3} = {\bf{0}}_m^{\text{T}}$可得

然后将式(C-2)代入式(C-1)中的等号左侧可知

考虑矩阵${\boldsymbol{A}}_1^{ - 1/2}{{\boldsymbol{A}}_2} \in {\mathbb{R}^{n \times (n - 1)}}$和向量${{\boldsymbol{\varPi }}^ \bot }[{\boldsymbol{A}}_1^{ - 1/2}{{\boldsymbol{A}}_2}]{\boldsymbol{A}}_1^{ - 1/2}{\boldsymbol{a}} \in {\mathbb{R}^{n \times 1}}$，他们的列数之和等于n，并且满足

于是有

此时利用正交投影矩阵的性质可知${\boldsymbol{\varPi }}[{{\boldsymbol{\varPi }}^ \bot }[{\boldsymbol{A}}_1^{ - 1/2}{{\boldsymbol{A}}_2}]{\boldsymbol{A}}_1^{ - 1/2}{\boldsymbol{a}}] = {{\boldsymbol{\varPi }}^ \bot }[{\boldsymbol{A}}_1^{ - 1/2}{{\boldsymbol{A}}_2}]$，最后结合该式与式(C-3)可知式(C-1)成立。证毕。

下面基于引理3证明式(89)。由于${\boldsymbol{W}}({\boldsymbol{z}},{\boldsymbol{d}}) \in {\mathbb{R}^{4 \times 4}}$是正定矩阵，${\boldsymbol{\tilde I}}_3^{\text{T}} \in {\mathbb{R}^{4 \times 3}}$是列满秩矩阵，并且根据矩阵${\boldsymbol{J}}({{\boldsymbol{u}}_{\text{o}}})$的定义式可以验证(C-6)关系式

此时利用引理3可知式(89)成立。

附录D　证明式(98)

首先根据正交投影矩阵的定义可知

然后结合矩阵${\boldsymbol{\bar \varOmega }}$的定义可得

另一方面，可以验证

以及${\mathrm{rank}}\left[ {{{({\boldsymbol{W}}({\boldsymbol{z}},{\boldsymbol{d}}))}^{1/2}}\dfrac{{\partial {{\boldsymbol{v}}_{\text{o}}}}}{{\partial {\boldsymbol{u}}_{\text{o}}^{\text{T}}}}} \right] = {\mathrm{rank}}\left[ {\dfrac{{\partial {{\boldsymbol{v}}_{\text{o}}}}}{{\partial {\boldsymbol{u}}_{\text{o}}^{\text{T}}}}} \right] = 3$，于是有

此时根据正交投影矩阵的性质可知

将式(D-2)和式(D-5)代入式(D-1)中，并且利用关系式${\boldsymbol{\bar I}}_3^{\text{T}}{{\boldsymbol{\varOmega }}_{\text{0}}}{{\boldsymbol{u}}_{\text{o}}} = {{\boldsymbol{\varOmega }}_{\text{1}}}{{\boldsymbol{v}}_{\text{o}}}$可得

将关系式${{\boldsymbol{\bar I}}_3}\dfrac{{\partial {{\boldsymbol{v}}_{\text{o}}}}}{{\partial {\boldsymbol{u}}_{\text{o}}^{\text{T}}}} = {{\boldsymbol{I}}_3}$和式(58)代入式(D-6)中可知式(98)成立。

附录E　证明式(99)

首先结合分块矩阵求逆公式和矩阵${\boldsymbol{\bar \varOmega }}$的定义可知

然后结合式(D-5)和式(58)可得

将式(E-2)和关系式${\boldsymbol{\bar I}}_3^{\text{T}}{{\boldsymbol{\varOmega }}_{\text{0}}}{{\boldsymbol{u}}_{\text{o}}} = {{\boldsymbol{\varOmega }}_{\text{1}}}{{\boldsymbol{v}}_{\text{o}}}$代入式(E-1)中可知

最后由式(E-3)可得

至此式(99)得证。