Processing math: 0%

      一种参考源修正的短波辐射源时差定位方法

      于家傲, 敬东, 夏冬玉, 夏楠, 杨洋, 张玉

      于家傲,敬东,夏冬玉,等. 一种参考源修正的短波辐射源时差定位方法[J]. 电波科学学报,2024,39(3):455-464. DOI: 10.12265/j.cjors.2023175
      引用格式: 于家傲,敬东,夏冬玉,等. 一种参考源修正的短波辐射源时差定位方法[J]. 电波科学学报,2024,39(3):455-464. DOI: 10.12265/j.cjors.2023175
      YU J A, JING D, XIA D Y, et al. A short-wave target time difference of arrival localization method corrected by calibration emitters[J]. Chinese journal of radio science,2024,39(3):455-464. (in Chinese). DOI: 10.12265/j.cjors.2023175
      Reference format: YU J A, JING D, XIA D Y, et al. A short-wave target time difference of arrival localization method corrected by calibration emitters[J]. Chinese journal of radio science,2024,39(3):455-464. (in Chinese). DOI: 10.12265/j.cjors.2023175

      一种参考源修正的短波辐射源时差定位方法

      详细信息
        作者简介:

        于家傲: (1989—),男,黑龙江齐齐哈尔人,工程师,博士,研究方向为信号处理、电子对抗技术、微波技术与天线。E-mail: yujiaao123@126.com

        敬东: (1969—),男,辽宁新民人,正高级工程师,博士,研究方向为新体制预警探测、多源信息融合等。E-mail: jingdong123@126.com

        通信作者:

        敬东 E-mail: jingdong123@126.com

      • 中图分类号: TN957.51

      A short-wave target time difference of arrival localization method corrected by calibration emitters

      • 摘要:

        为解决短波辐射源到达时间差 (time difference of arrival, TDOA) 定位(简称时差定位)方法受电离层影响导致的定位精度下降的问题,提出了一种利用参考源修正的短波辐射源目标时差定位方法。针对地球表面短波辐射源,基于电离层球面反射模型的电离层反射虚高近似方法,建立了利用参考修正的短波目标时差定位模型。考虑参考源与目标共用电离层反射区域对电离层虚高的影响,将各电离层反射点的距离相关性引入电离层虚高的协方差矩阵中,实现了目标定位精度的修正。通过推导和仿真所提模型的克拉美·罗下界,分析了参考源修正目标定位精度的可行性。进一步给出基于Armijo直线搜索Newton法的最大似然估计方法,通过仿真数据验证了所提算法的有效性,实现了良好的定位效果。

        Abstract:

        A time difference of arrival(TDOA) localization method is proposed for short-wave target corrected by calibration emitters. Based on the ionospheric reflection virtual height of spherical reflection model, a TDOA localization model is established for short-wave source targets and calibration emitters on the Earth’s surface. Considering the influence of ionospheric reflection region shared by calibration emitters and target, the distance correlation of each ionospheric reflection point is introduced into the covariance matrix of ionospheric virtual height, and the target localization accuracy is corrected. The feasibility of the target localization accuracy corrected by calibration emitters is analyzed by deriving and simulating the Cramer-Rao lower bound. Further, a maximum likelihood estimation method is proposed based on Armijo’s linear search Newton method. Effectiveness of the proposed method is verified by simulation data, and a good localization effect is achieved.

      • 短波辐射源目标侦察通过信号的电离层折射实现超视距探测,具有探测范围大、作用距离广的特点,广泛应用于民用电磁环境监测、国防等重要领域。到达时间差 (time difference of arrival, TDOA) 定位(简称时差定位)方法对直视条件下的辐射源目标实现了较好的定位效果,但短波信号在超视距传输条件下受电离层空间非均匀性和时变特性等影响严重,定位效果差,限制了时差定位技术对短波辐射源目标的应用效果[1-2]

        目前短波时差定位问题通过电离层球面反射模型、电离层准抛物模型和国际参考电离层(International Reference Ionosphere, IRI)模型模拟信号传输路径来建立短波时差定位模型,定位模型通过凸优化、约束总体最小二乘等迭代方法,或加权多维标度、两步加权最小二乘等闭式解方法实现目标定位参数的估计求解[3-6]。文献[7-8]提出了直视条件下利用参考源修正的时差定位闭式算法,验证了利用参考源可提高定位精度的有效性。经电离层传播的短波时差定位模型中包含电离层参数,导致模型非线性化,限制了闭式算法对短波目标的定位效果。文献[9]利用IRI模型获取更精确的电子密度信息,建立了以数值射线追踪为基础的短波时差定位模型,并通过克拉美·罗下界(Cramer-Rao lower bound, CRLB)分析了定位结果的有效性。文献[10]提出了利用随机空间谱估计的短波时差定位方法,提高了电离层信号多跳传输条件下的目标定位精度。文献[11-12]提出了利用电离层虚高参数的短波时差定位的方法,在假设基站间距较小、信号传输路径电离层虚高一致的条件下实现了较好的定位效果,但接收站基线扩大后其定位精度下降明显。

        针对接收站基线扩大、电离层非均匀性导致定位精度降低的问题,本文提出一种目标和参考源联合估计的短波时差定位方法。基于地球表面目标信号的电离层球面反射模型,利用目标、参考源到达接收站的时差,以及各传输路径电离层虚高作为观测数据,建立了利用参考源修正的目标时差定位模型。相比于假设目标到各接收站的电离层反射点虚高完全相关,本文将反射点共用电离层区域的距离相关性引入电离层虚高协方差矩阵中,通过联合估计实现了对目标定位精度的修正。通过理论推导和仿真分析所提模型的CRLB,验证了参考源修正目标定位精度的有效性。进一步给出了基于Armijo直线搜索Newton法的最大似然估计法进行目标定位参数估计,并通过仿真结果验证了所提算法的有效性。

        文献[11]验证了球面反射模型对解决电离层信号传播问题的有效性,为便于理论分析和模型设计,本文采用电离层球面反射模型建立包含目标、接收站和参考源的时差定位场景,如图1所示。第i个接收站位置{{\boldsymbol{s}}_i} = ({x_{{{\boldsymbol{{s}}}_i}}},{y_{{{\boldsymbol{s}}_i}}},{{\textit{z}}_{{{\boldsymbol{s}}_i}}}),目标位置{\boldsymbol{u}} = ({x_{\boldsymbol{u}}},{y_{\boldsymbol{u}}},{{\textit{z}}_{\boldsymbol{u}}}),第j个参考源位置{{\boldsymbol{c}}_j} = ({x_{{{{{\boldsymbol{c}}}}_{_j}}}},{y_{{{\boldsymbol{c}}_{_j}}}},{{\textit{z}}_{{{\boldsymbol{c}}_{_j}}}})。地球为等效半径R的球体,接收站、目标和参考源均位于地球表面,根据约束条件{{\boldsymbol{s}}_i^{\text{T}}}{{\boldsymbol{s}}_i} = {{\boldsymbol{u}}^{\text{T}}}{\boldsymbol{u}}{\text{ = }}{{\boldsymbol{c}}_j^{\text{T}}}{{\boldsymbol{c}}_j} = {R^2},位置矢量{\boldsymbol{u}}{{\boldsymbol{s}}_i}{{\boldsymbol{c}}_j}中仅有x,y两个分量为自变量,{\textit{z}}分量为因变量,为避免解模糊,设{\textit{z}} > 0

        图  1  含有参考源的时差定位场景
        Fig.  1  Scenario of TDOA localization with calibration emitter

        接收站{{\boldsymbol{s}}_i}收到的目标 {\boldsymbol{u}} 的信号经电离层球面反射点{{\boldsymbol{P}}_{\boldsymbol{u}}}反射,电离层虚高为{h_i},其传输距离 {r_i}

        {r_i} = 2{( {{R^2} + {{\left( {R + {h_i}} \right)}^2} - \left( {R + {h_i}} \right)\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_i}} \right\|} )^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} (1)

        以接收站{{\boldsymbol{s}}_1}为参考,则各接收站接收到目标信号并进行相干处理后得到的距离差为

        {r_{i1}} = {r_i} - {r_1} \text{,}i = 2,3, \cdots ,M (2)

        式(2)为目标距离差状态方程。将该距离差合并表示为 M - 1 列状态向量 {{\boldsymbol{r}}^{\mathrm{o}}} ,则目标距离差观测值可表示为服从正态分布的随机向量:

        {\boldsymbol{r}} \sim {\text{N}}\left( {{{\boldsymbol{r}}^{\mathrm{o}}},{{\boldsymbol{Q}}_{\text{α}}}} \right) (3)

        其均值为 {{\boldsymbol{r}}^{\mathrm{o}}} ,协方差矩阵为{{\boldsymbol{Q}}_{\text{α}} }{{\boldsymbol{Q}}_{\text{α}}}M - 1阶方阵,其对角元素为\sigma _{\mathrm{t}}^2{c^2},其余元素为0.5\sigma _{\mathrm{t}}^2{c^2}\sigma _{\mathrm{t}}^2为时差测量方差,c为光速。目标时差的观测值即可表示为{\boldsymbol{\tau }} = {\boldsymbol{r}}/c

        类似的,接收站{{\boldsymbol{s}}_i}收到的辐射参考点 {{\boldsymbol{c}}_j} 的信号经球面反射点{{\boldsymbol{P}}_{{{\boldsymbol{c}}}}}反射,电离层虚高为{h_{i,j}},其传输距离 {\tilde r_{i,j}}

        {\tilde r_{i,j}} = 2{( {{R^2} + {{\left( {R + {h_{i,j}}} \right)}^2} - \left( {R + {h_{i,j}}} \right)\| {{{\boldsymbol{c}}_j} + {{\boldsymbol{s}}_i}} \|} )^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} (4)

        接收站接收到N个参考源信号的距离差为

        {\tilde r_{i1,j}} = {\tilde r_{i,j}} - {\tilde r_{1,j}} \text{,}i = 2,3, \cdots ,M\text{,}j = 1,2, \cdots ,N (5)

        式(5)为参考源的距离差状态方程。将该距离差合并表示为N\left( {M - 1} \right)列状态向量 {{\boldsymbol{\tilde r}}^{\mathrm{o}}} ,则参考源距离差观测值可表示为服从正态分布的随机向量:

        {\boldsymbol{\tilde r}} \sim {\text{N}}\left( {{{{\boldsymbol{\tilde r}}}^{\mathrm{o}}},{{\boldsymbol{Q}}_{\mathrm{c}} }} \right) (6)

        其均值为 {{\boldsymbol{\tilde r}}^{\mathrm{o}}} ,协方差矩阵为{{\boldsymbol{Q}}_{\mathrm{c}}}{{\boldsymbol{Q}}_{\mathrm{c}}}N\left( {M - 1} \right)阶方阵,其N个对角线子方阵为{\left[ {{{\boldsymbol{Q}}_{\mathrm{c}}}} \right]_{ii}} = {{\boldsymbol{Q}}_ {\text{α}} }。参考源时差的观测值即可表示为{\boldsymbol{\tilde \tau }} = {\boldsymbol{\tilde r}}/c

        合并目标和参考源信号路径的电离层虚高表示为 \left( {N + 1} \right)M 列状态向量 {{\boldsymbol{h}}^{\mathrm{o}}}

        {{\boldsymbol{h}}^{\mathrm{o}}} = \left[ {{h_1}, \cdots ,{h_M},{h_{1,1}}, \cdots ,{h_{M,1}}, \cdots ,{h_{M,N}}} \right] (7)

        其观测值可表示为服从正态分布的随机向量:

        {\boldsymbol{h}} \sim {\text{N}}\left( {{{\boldsymbol{h}}^{\mathrm{o}}},{Q_{\mathrm{h}}}} \right) (8)

        其均值为 {{\boldsymbol{h}}^{\mathrm{o}}} ,协方差矩阵{{\boldsymbol{Q}}_{\mathrm{h}}}\left( {N + 1} \right)M阶方阵。

        短波时差定位模型中,对电离层虚高的测量和估计是影响辐射源目标定位精度的关键,参考电离层理论模型分析和电离层虚高测量经验数据分析[5-6, 10],在一定时间和频段范围内,当基站间距离小于某一门限(通常为30~100 km)时,同一地球表面短波辐射源到达各基站的反射路径对应的电离层虚高可认为具有较高的相关性,随着距离进一步增加,电离层虚高的相关性逐渐下降。同理,当目标辐射源与参考源距离靠近时,其对应的电离层虚高具有类似的距离相关性。当辐射源目标信号传输路径与参考源信号传输路径共用同一区域电离层时,利用信号路径电离层虚高的相关性可实现对目标定位结果的修正。为近似表征电离层虚高与反射点距离的相关性变化并不失一般性,本文假设该相关系数\rho 随反射点相对距离的函数为

        \rho = \exp \left( {\frac{{ - {{\left\| {{{\boldsymbol{P}}_{\boldsymbol{u}}} - {{\boldsymbol{P}}_{\boldsymbol{c}}}} \right\|}^2}}}{{{R_{{{\mathrm{cov}}} }^2}}}} \right) (9)

        设电离层相关性距离参数 {R_{{{\mathrm{cov}}} }} = 200 km,该函数表征的电离层虚高相关系数随反射点相对距离变化曲线如图2所示。

        图  2  电离层虚高相关系数随反射点相对距离变化
        Fig.  2  Ionospheric virtual height correlation coefficient vs. relative distance of reflection points

        利用该距离相关性,电离层虚高协方差矩阵中的元素 {\left[ {{{\boldsymbol{Q}}_{\mathrm{h}}}} \right]_{p,q}} 可表示为

        \begin{split} & {\left[ {{{\boldsymbol{Q}}_{\mathrm{h}}}} \right]_{p = mM + i,q = nM + j}} \\ & =\left\{ {\begin{array}{*{20}{l}} {\sigma _{\mathrm{h}}^2\exp \left[ {\dfrac{{ - {{\left\| {{{\boldsymbol{s}}_i} - {{\boldsymbol{s}}_j}} \right\|}^2}}}{{{{\left( {{R_{{{\mathrm{cov}}} }}/2} \right)}^2}}}} \right]}&{n = m} \\ {\sigma _{\mathrm{h}}^2\exp \left[ {\dfrac{{ - {{\left\| {{{\boldsymbol{c}}_m} - {\boldsymbol{u}} + {{\boldsymbol{s}}_i} - {{\boldsymbol{s}}_j}} \right\|}^2}}}{{{{\left( {{R_{{{\mathrm{cov}}} }}/2} \right)}^2}}}} \right]}&{m > 0,n = 0} \\ {\sigma _{\mathrm{h}}^2\exp \left[ {\dfrac{{ - {{\left\| {{\boldsymbol{u}} - {{\boldsymbol{c}}_n} + {{\boldsymbol{s}}_i} - {{\boldsymbol{s}}_j}} \right\|}^2}}}{{{{\left( {{R_{{{\mathrm{cov}}} }}/2} \right)}^2}}}} \right]}&{n > 0,m = 0} \\ {\sigma _{\mathrm{h}}^2\exp \left[ {\dfrac{{ - {{\left\| {{{\boldsymbol{c}}_m} - {{\boldsymbol{c}}_n} + {{\boldsymbol{s}}_i} - {{\boldsymbol{s}}_j}} \right\|}^2}}}{{{{\left( {{R_{{{\mathrm{cov}}} }}/2} \right)}^2}}}} \right]}&{m,n > 0,m \ne n} \end{array}} \right. \end{split} (10)

        式中, \sigma _{\mathrm{h}}^2 为电离层虚高方差。当 n = m 时,协方差矩阵反映了目标信号到各接收站传输路径电离层虚高的相关情况;当 m > 0,n=0 n > 0,m = 0 时,协方差矩阵反映了目标信号及参考源到各接收站传输路径的电离层虚高相关情况;当 m,n > 0,m \ne n 时,协方差矩阵反映了各参考源到各接收站传输路径的电离层虚高相关情况。即使协方差矩阵{{\boldsymbol{Q}}_{\mathrm{h}}}存在{\boldsymbol{h}}中的两分量完全相关,也不意味着两值相等。

        通过估计参数的CRLB来表征所提时差定位模型的定位精度。针对本文短波时差定位场景模型,设接收站位置和已知参考源位置为常量,其CRLB可表示为通过观测量 {\boldsymbol{r}} {\boldsymbol{\tilde r}} {\boldsymbol{h}}实现对参数 {{\boldsymbol{u}}^{\mathrm{o}}} {{\boldsymbol{h}}^{\mathrm{o}}} 无偏估计的最小方差值,该值等于Fisher信息矩阵的逆。Fisher信息矩阵通过计算观测量联合概率密度函数二阶偏导数的期望得到,通过计算Fisher信息矩阵的逆可得到所提模型的CRLB[4]。Fisher信息矩阵可表示为

        {\text{FIM}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{X}}&{\boldsymbol{Y}} \\ {{{\boldsymbol{Y}}^{\text{T}}}}&{\boldsymbol{Z}} \end{array}} \right] (11)

        {\boldsymbol{X}}子阵中包含的{{\boldsymbol{Q}}_{\mathrm{h}}} {{\boldsymbol{u}}^{\mathrm{o}}} 的函数,与传统方差矩阵为常数项矩阵不同,此处须对{{\boldsymbol{Q}}_{\mathrm{h}}}求偏导,故有

        {\boldsymbol{X}}= {\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}} \right) + {{\boldsymbol{X}_{{{\boldsymbol{Q}}}}'}} (12)
        {{\boldsymbol{X}}_{\boldsymbol{Q}}'} = \dfrac{1}{2}\left[ {\begin{array}{*{20}{c}} {{\text{tr}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {x_{\boldsymbol{u}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {x_{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)}&{{\text{tr}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {x_{\boldsymbol{u}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {y_{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)} \\ {{\text{tr}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {x_{\boldsymbol{u}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {y_{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)}&{{\text{tr}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {y_{\boldsymbol{u}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {y_{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)} \end{array}} \right] (13)
        {\boldsymbol{Y}} = {\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right) (14)
        {\boldsymbol{Z}} = {\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right) + {\left( {\frac{{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\mathrm{c}}^{ - 1}\left( {\frac{{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right) + {\boldsymbol{Q}}_{\mathrm{h}}^{ - 1} (15)

        式(12)~(15)中各项表达式详见附录A。

        为分析参考源对定位精度的影响,比较分析不包含参考源和包含参考源对目标CRLB的影响。不考虑参考源及其距离相关性对协方差矩阵{\boldsymbol{Q}}_{\mathrm{h}}^{}的影响时,FIM矩阵表示为

        {\text{FIM}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{X}}_0}}&{\boldsymbol{Y}} \\ {{{\boldsymbol{Y}}^{\text{T}}}}&{{{\boldsymbol{Z}}_0}} \end{array}} \right] (16)
        {\boldsymbol{X}}_0= {\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}} \right) (17)
        {{\boldsymbol{Z}}_0} = {\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {\frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{{h}}}^{\mathrm{o}}}}}} \right) + \sigma _{\mathrm{h}}^{ - 2}{\boldsymbol{I}} (18)
        {\text{CRLB}'_{\boldsymbol{u}}}= {{\boldsymbol{X}}_0^{ - 1}} + {{\boldsymbol{X}}_0^{ - 1}}{\boldsymbol{Y}}{\left( {{{\boldsymbol{Z}}_0} - {{\boldsymbol{Y}}^{\text{T}}}{{\boldsymbol{X}}_0^{ - 1}}{\boldsymbol{Y}}} \right)^{ - 1}}{{\boldsymbol{Y}}^{\text{T}}}{{\boldsymbol{X}}_0^{ - 1}} (19)

        包含参考源时差定位精度的CRLB可表示为

        \begin{split} {\text{CRL}}{{\text{B}}_{\boldsymbol{u}}} = & {\left( {{{\boldsymbol{X}}_0} + {{{\boldsymbol{X'}_Q}}}} \right)^{ - 1}} + {\left( {{{\boldsymbol{X}}_0} + {{{\boldsymbol{X'}_Q}}}} \right)^{ - 1}}{\boldsymbol{Y}} \\ & {\left[ {{{\boldsymbol{Z}}_0} + {\boldsymbol{\varGamma }} - {{\boldsymbol{Y}}^{\text{T}}}{{\left( {{\boldsymbol{X}} + {{{\boldsymbol{X'}_Q}}}} \right)}^{ - 1}}{\boldsymbol{Y}}} \right]^{ - 1}} {{\boldsymbol{Y}}^{\text{T}}}{\left( {{\boldsymbol{X}} + {{{\boldsymbol{X'}_Q}}}} \right)^{ - 1}} \end{split} (20)
        {\boldsymbol{\varGamma }} = {\left( {\frac{{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\mathrm{c}}^{ - 1}\left( {\frac{{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right) + \left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1} - \sigma _{\mathrm{h}}^{ - 2}{\boldsymbol{I}}} \right) (21)

        根据式(13)可得 {{\boldsymbol{X}'_{\boldsymbol{Q}}}} 正定,{\left( {{{\boldsymbol{X}}_0} + {{{\boldsymbol{X}}}'_{\boldsymbol{Q}}}} \right)^{ - 1}} < {\boldsymbol{X}}_0^{ - 1},则

        {\text{CRL}}{{\text{B}}_{\boldsymbol{u}}} \leqslant {\boldsymbol{X}}_0^{ - 1} + {\boldsymbol{X}}_0^{ - 1}{\boldsymbol{Y}}{\left[ {{{\boldsymbol{Z}}_0} - {{\boldsymbol{Y}}^{\text{T}}}{{\boldsymbol{X}}_0^{ - 1}}{\boldsymbol{Y}} + {\boldsymbol{\varGamma }}} \right]^{ - 1}}{{\boldsymbol{Y}}^{\text{T}}}{{\boldsymbol{X}}_0^{ - 1}} (22)

        根据式(14), {\boldsymbol{Y}} 的前m行子阵 {{\boldsymbol{Y}}_m} 为满秩,当 {\boldsymbol{Q}}_{\mathrm{h}}^{ - 1} \ne \sigma _{\mathrm{h}}^{ - 2}{\boldsymbol{I}} 时,参照文献[13],根据定义(10)总有 {{\boldsymbol{Y}}_m}{\boldsymbol{\varGamma }}{{\boldsymbol{Y}}_m^{\text{T}}} 正定,使得

        {\text{CRL}}{{\text{B}'_{\boldsymbol{u}}}} \geqslant {\boldsymbol{X}}_0^{ - 1} + {\boldsymbol{X}}_0^{ - 1}{\boldsymbol{Y}}{\left[ {{{\boldsymbol{Z}}_0} - {{\boldsymbol{Y}}^{\mathrm{T}}}{{\boldsymbol{X}}_0^{ - 1}}{\boldsymbol{Y}} + {\boldsymbol{\varGamma }}} \right]^{ - 1}}{{\boldsymbol{Y}}^{\mathrm{T}}}{{\boldsymbol{X}}_0^{ - 1}} (23)

        则有 {\text{CRL}}{{\text{B}'_{\boldsymbol{u}}}} \geqslant {\text{CRL}}{{\text{B}}_{\boldsymbol{u}}} ,即不含有参考源的目标定位误差大于含参考源的目标定位误差,证明了所提方法利用参考源修正目标定位精度的可行性。

        当目标逐渐靠近参考源,辐射源目标信号路径电离层虚高与参考源信号路径电离层虚高的相关性逐渐增加,并体现为正定矩阵 {\boldsymbol{Q}}_{\mathrm{h}}^{} 中非对角线元素值增加,使得式(22)有 {{\boldsymbol{Y}}_m}{\boldsymbol{\varGamma }}{{\boldsymbol{Y}}_m^{\text{T}}} 逐渐变大, {\text{CRL}}{{\text{B}}_{\boldsymbol{u}}} 值降低,即目标靠近参考源降低了目标定位误差。反之当目标远离参考源时,各路径电离层虚高相关性降低并趋于相互独立,并体现为 {\boldsymbol{Q}}_{\mathrm{h}}^{} \to \sigma _{\mathrm{h}}^2{\boldsymbol{I}} ,则有 {{\boldsymbol{Y}}_m}{\boldsymbol{\varGamma }}{{\boldsymbol{Y}}_m^{\text{T}}} \to {{\textit{0}}} {\text{CRL}}{{\text{B}}_{\boldsymbol{u}}} \to {\text{CRL}}{{\text{B}'_{\boldsymbol{u}}}} ,即目标定位误差趋近于不含参考源的目标定位方式,目标定位精度得不到来自参考源的修正。

        针对本文短波时差定位场景模型,基于参考源修正的目标定位状态方程由式(2)、(5)、(7)给出。根据式(3)、(6)、(8)的概率形式,并设接收站位置{{\boldsymbol{s}}_i}和参考源位置{{\boldsymbol{c}}_j}为常量,则观测值的联合概率密度函数及其对数化可表示为

        \begin{split} p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{r}}^{\mathrm{o}}},{{{\boldsymbol{\tilde r}}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right) = & p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right) \\ = & p\left( {{\boldsymbol{r}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right)p\left( {{\boldsymbol{\tilde r}};{{\boldsymbol{h}}^{\mathrm{o}}}} \right)p\left( {{\boldsymbol{h}};{{\boldsymbol{h}}^{\mathrm{o}}}} \right) \end{split} (24)
        \begin{split} & \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right) = - \frac{1}{2} \\ & [ {\left( {M - 1} \right)\ln 2{\text{π }} + \ln \left| {{{\boldsymbol{Q}}_{\text{α}} }} \right| + {{\left( {{\boldsymbol{r}} - {{\boldsymbol{r}}^{\mathrm{o}}}} \right)}^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {{\boldsymbol{r}} - {{\boldsymbol{r}}^{\mathrm{o}}}} \right)} \\ & +N\left( {M - 1} \right)\ln 2{\text{π }} + \ln \left| {{{\boldsymbol{Q}}_{\mathrm{c}}}} \right| + {\left( {{\boldsymbol{\tilde r}} - {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}} \right)^{\text{T}}}{\boldsymbol{Q}}_{\mathrm{c}}^{ - 1}\left( {{\boldsymbol{\tilde r}} - {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}} \right) \\ & + {N\left( {M - 1} \right)\ln 2{\text{π }} + \ln \left| {{{\boldsymbol{Q}}_{\mathrm{h}}}} \right| + {{\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)}^{\text{T}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)} ] \end{split} (25)

        根据最大似然估计方法,当估计的参数 {{\boldsymbol{u}}^{\mathrm{o}}} {{\boldsymbol{h}}^{\mathrm{o}}} 使得观测值 {\boldsymbol{r}} {\boldsymbol{\tilde r}} {\boldsymbol{h}}的概率最大时,该估计值 {{\boldsymbol{u}}^{\mathrm{o}}} 即为所提定位模型的定位结果:

        {{\boldsymbol{\hat u}}^{\mathrm{o}}} = \arg \mathop {\sup }\limits_{{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right) (26)

        在最大似然估计过程中,当目标远离参考源时,根据式(10)电离层虚高方差矩阵 {\boldsymbol{Q}}_{\mathrm{h}}^{} 的非对角线元素几乎为0,各信号路径的电离层虚高观测值相互独立,导致参考源的估计参数和目标的估计参数缺乏相关性,从而未能对最大似然概率产生影响,目标定位参数得不到修正。随着目标靠近参考源,矩阵 {\boldsymbol{Q}}_{\mathrm{h}}^{} 的非对角线元素增加,参考源信号路径电离层虚高与目标信号路径电离层虚高相关性增加,使得目标参数估计结果同时满足目标测量值和参考源测量值条件下的似然概率最大,从而实现了参考源对目标定位结果的修正。因式(9)定义的相关性是相对距离的连续性函数,该优化过程无须给定相对距离判断门限即可完成自动迭代优化过程。

        为实现对所提时差定位模型目标定位参数的最大似然估计,针对电离层虚高协方差矩阵{{\boldsymbol{Q}}_{\mathrm{h}}}是变量{\boldsymbol{u}}的函数并与参数{{\boldsymbol{s}}_i}{{\boldsymbol{c}}_j}有关,经典的两步加权最小二乘(two-step weighted least squares, TS-WLS)算法难以实现有效求解的问题,本文提出一种基于Armijo直线搜索Newton法的最大似然估计方法。

        合并目标位置和电离层虚高的参数估计向量,即 {\boldsymbol{\theta }} = {\left[ {{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right]^{\text{T}}} {{\boldsymbol{u}}^{\mathrm{o}}} 仅包含{\boldsymbol{u}}的2个x,y分量,残差量 {{\boldsymbol{b}}_k} 为观测量联合概率密度函数对估计量的导数。当估计量 {\boldsymbol{\theta }} 为最大似然估计时,观测量联合概率密度函数得到的观测量出现概率为最大。同时根据最大似然估计方法,当该残差量 {{\boldsymbol{b}}_k} 0矢量时估计量为最大似然估计。本文采用Newton法迭代求解该估计量,利用Fisher矩阵代替Newton法中的Hessian矩阵[14-15],则估计量由k步到k + 1步的最大似然估计迭代公式为

        {{\boldsymbol{\theta }}_{k + 1}} = {{\boldsymbol{\theta }}_k} + \mu {\text{FI}}{{\text{M}}^{ - 1}} \left( {{{\boldsymbol{\theta }}_k}} \right){{\boldsymbol{b}}_k} (27)

        {{\boldsymbol{b}}_k} 中各项表达式见附录B。

        为避免迭代过程中出现奇异值问题,采用小扰动正则化方法对协方差矩阵{{\boldsymbol{Q}}_{\mathrm{h}}}附加一个小扰动对角阵\delta {\boldsymbol{I}}。为避免迭代发散并保证估计值在合理的范围内,基于Armijo直线搜索法,对更新向量乘以步长\mu \in (0,1],使得更新向量对最大似然目标函数充分升高,即满足如下Armijo条件不等式:

        \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{\theta }}_{k + 1}}} \right) > \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{\theta }}_k}} \right) + \alpha \mu {\boldsymbol{b}}_k^{\mathrm{T}}{\text{FI}}{{\text{M}}^{ - 1}}\left( {{{\boldsymbol{\theta }}_k}} \right){{\boldsymbol{b}}_k} (28)

        迭代流程如图3所示,具体算法实现步骤如下:

        1)设置包含目标、参考源位置和电离层虚高的估计量初始值,并设迭代步数k = 1

        2)根据目标、参考源位置对电离层虚高进行观测或基于电离层模型仿真计算。

        3)根据式(10)计算电离层虚高协方差矩阵{{\boldsymbol{Q}}_{\mathrm{h}}},并附加小扰动对角阵\delta {\boldsymbol{I}}使其正则化。

        4)根据式(11)~(15)计算Fisher信息矩阵,根据式(27)计算残差向量 {{\boldsymbol{b}}_k}

        5)判断更新步估计量是否满足式(28)的Armijo约束条件,其中 \alpha \in (0,0.5] 。若不满足则回调步长\mu = \beta \mu ,其中回调系数\; \beta \in \left( {0,1} \right) ,并重新执行步骤4;若满足条件则得到合适的步长\mu

        6)根据迭代公式(27)计算并更新k + 1步估计量。

        7)判断迭代步数是否达到最大值或满足残差量小于阈值的收敛条件。若判断为真,则迭代结束并得到最终估计量;判断结果为否,则基于更新估计量返回执行步骤2继续迭代,并另迭代步数k = k + 1

        图  3  基于Armijo直线搜索Newton法的最大似然估计流程
        Fig.  3  Flowchart of maximum likelihood estimation based on Armijo’s linear search Newton method

        在实际应用中,参考源与接收站之间传输路径的电离层虚高观测值,可由根据IRI模型采用射线追踪技术进行反演解算,或通过位于反射点下方的电离层测高仪固定站测量。由于目标到接收站之间传输路径的电离层虚高未知,该电离层虚高观测值可初设为经验值,当所提迭代算法的目标位置结果收敛后,根据目标点位置对该电离层虚高观测值进行更新,并继续迭代至收敛。

        为验证参考源修正目标定位方法在球面反射电离层模型条件下的定位效果,基于CRLB仿真分析不同参数条件对目标定位精度的影响,并仿真验证所提基于Armijo直线搜索Newton法的最大似然估计方法的定位效果。

        仿真基于图1所示的以地心为原点三维坐标系,单位为km。设地球为圆球模型,接收站、参考源、目标位置及序号设置见表1,因各点位置在地球半径为R的球面上,表1中仅给出位置的xy坐标分量,且分量z>0。

        表  1  接收站和参考源以及目标位置坐标
        Tab.  1  Coordinate position of receiving stations, calibration emitters and their target location
        序号 接收站/km 参考源/km 目标位置/km
        1 [0, 0] [500 \sqrt2 ,500 \sqrt2] [900, 900]
        2 [1000, 0] [750 \sqrt2, 750 \sqrt2] [750, 750]
        3 [0, 1000] [500 \sqrt3, 500] [710, 710]
        4 [500, 500] [500, 500 \sqrt3] [1200, 1200]
        5 [200, 500] [ 750 \sqrt3, 750] [1550, 800]
        6 [500, 200] [750, 750 \sqrt3] [950, 800]
        下载: 导出CSV 
        | 显示表格

        在电离层传输条件下,设接收站位置和已知参考源位置具有精度较高的测量值,因此主要对电离层虚高方差 \sigma _{\mathrm{h}}^2 、时间测量方差\sigma _{\mathrm{t}}^2,以及接收站和参考源位置分布进行分析。

        1)电离层虚高方差 \sigma _{\mathrm{h}}^2 影响分析

        选取接收站1,2,3,4和参考源1,2,对不同位置的目标1,2,3进行仿真,考虑时间测量精度约为15 ns,可设 \sigma _{\mathrm{t}}^2{c^2} = 2 \times {10^{ - 5}}\;{\text{k}}{{\text{m}}^2} c为光速。分析电离层虚高方差 \sigma _{\mathrm{h}}^2 对定位精度的影响,定位精度通过计算 {\text{CRL}}{{\text{B}}_{\boldsymbol{u}}} 迹的对数值 {\text{10lg}}\left( {{\text{tr}}\left( {{\text{CRL}}{{\text{B}}_{\boldsymbol{u}}}} \right)} \right) 来表征。仿真结果如图4所示,随着电离层虚高方差的增加,含参考源和无参考源的目标定位精度逐渐变差。由于目标1距离参考源较远,其目标定位精度变化与无参考源基本吻合,随着目标2,3逐渐靠近参考源1其定位精度明显提高。在本文所设仿真条件下,在 \sigma _{\mathrm{h}}^2 = {10^2} 时,靠近参考源1的目标3理论定位精度可提高约20 dB。

        图  4  电离层虚高方差 \sigma _{\mathrm{h}}^2 对定位精度的影响
        Fig.  4  Influence of ionospheric virtual height variance \sigma _{\mathrm{h}}^2 on localization accuracy

        2)时间测量方差\sigma _{\mathrm{t}}^2影响分析

        选取接收站1,2,3,4和参考源1,2,对不同位置的目标1,2,3进行仿真,电离层虚高方差为 \sigma _{\mathrm{h}}^2 = {10^2}\;{\text{k}}{{\text{m}}^2} ,分析时间测量方差\sigma _{\mathrm{t}}^2对CRLB的影响。仿真结果如图5所示,对于目标1由于其距离参考源较远,电离层虚高方差对定位精度影响较大,\sigma _{\mathrm{t}}^2的变化对定位精度影响较小,随着目标2,3逐渐靠近参考源1,电离层虚高的影响被参考源修正,\sigma _{\mathrm{t}}^2的降低明显改善了定位精度。

        图  5  时间测量方差\sigma _{\mathrm{t}}^2对定位精度的影响
        Fig.  5  Influence of time measurement variance \sigma _{\mathrm{t}}^2 on localization accuracy

        3)接收站和参考源分布影响分析

        仿真分析不同接收站和参考源分布条件下2 000 km×2 000 km区域内目标定位精度,定位精度的等高线分布结果如图6所示,该图为地球曲面的局部展开面。当无参考源时,选取4个接收站(1~4号)和6个接收站(1~6号),该区域内目标定位精度基本不变,如图6(a)、(b)所示;当设置2个参考源(1,2号)时,接收站数量由4个增加至6个,目标定位精度 {\text{10lg}}\left( {{\text{tr}}\left( {{\text{CRL}}{{\text{B}}_{\boldsymbol{u}}}} \right)} \right) < 30\;{\text{dB}} 的区域明显扩大至1 200 km×1 200 km,如图6(c)、(d)所示;当设置4个参考源(3~6号)时,接收站数量由4个增加至6个,目标定位精度 {\text{10lg}}\left( {{\text{tr}}\left( {{\text{CRL}}{{\text{B}}_{\boldsymbol{u}}}} \right)} \right) <30\;{\text{dB}} 的区域扩大至1 400 km×1 400 km,区域内定位精度提高,由于参考源与目标共用的电离层反射区域增加,使得无参考源邻近的部分区域(如(900, 900)附近)同样得到较高的定位精度,如图6(e)、(f)所示。

        图  6  接收站和参考源对区域目标定位精度影响
        Fig.  6  Influence of the distribution of receiving stations and calibration emitters on localization accuracy

        算法场景1选取4个接收站(1~4号)和2个参考源(1,2号),采用无约束和Armijo约束的Newton法的最大似然估计方法分别对目标1,3的位置进行计算。设迭代目标初始值为{{\boldsymbol{u}}^{\mathrm{o}}} = (1 \; 000,1 \; 000),向量 {{\boldsymbol{h}}^{\mathrm{o}}} 中各元素初始值均为150,算法迭代步数为100,小扰动对角阵系数\delta = {10^{ - 8}}。Armijo约束条件中 \alpha = 0.5 ,回调初始步长 \mu = 1 ,回调系数\;\beta = 0.2。两种算法的迭代过程对比如图7所示,随着迭代步数的增加目标定位精度的CRLBu和相对真实目标距离(绝对误差)逐渐收敛。在迭代过程中式(27)中的FIM−1项受 {{\boldsymbol{Q}}_{\mathrm{h}}} 的不断更新影响,随着目标3估计值靠近参考源1, {{\boldsymbol{Q}}_{\mathrm{h}}} 在迭代过程中变化明显,无约束的Newton法易产生迭代抖动和收敛速度下降的现象。采用Armijo条件不等式后,约束更新向量指向最大似然目标函数的升高方向,实现了更好的收敛效果。比较目标1,3进行迭代收敛后的 {\text{10lg}}\left( {{\text{tr}}\left( {{\text{CRL}}{{\text{B}}_{\boldsymbol{u}}}} \right)} \right) 值分别为27.39,15.44 dB,单次仿真的相对真实目标距离分别为23.90,11.11 km。通过1 000次蒙特卡洛随机观测值仿真解算定位结果均方误差(mean squared error, MSE)的对数值分别为28.14,16.53 dB,并接近CRLB,靠近参考源的目标实现了更高的定位精度,验证了所提方法利用参考源改善目标定位精度的有效性。

        图  7  算法场景1条件下所提估计方法单次仿真迭代过程
        Fig.  7  Iterative process of the proposed estimation method for algorithm scenario 1 in a single simulation

        算法场景2选取6个接收站(1~6号)和4个参考源(3~6号),采用所提估计方法对目标4,5,6的位置进行计算。通过所提算法对场景2目标4,5,6进行迭代求解如图8所示,迭代收敛后的 {\text{10lg}}\left( {{\text{tr}}\left( {{\text{CRL}}{{\text{B}}_{\boldsymbol{u}}}} \right)} \right) 值分别为30.81,25.48,19.50 dB,单次仿真的相对真实目标距离分别为35.26,28.53,13.08 km。通过1 000次蒙特卡洛仿真解算定位结果MSE的对数值分别为32.34,27.78,22.13 dB。虽然场景2中目标未靠近参考源,但参考源与目标共用电离层反射区域的影响,同样实现了定位精度的提高。

        图  8  算法场景2条件下所提估计方法单次仿真迭代过程
        Fig.  8  Iterative process of the proposed estimation method for algorithm scenario 2 in a single simulation

        为进一步验证所提算法的适用性,采用基于国际电离层模型IRI2016的PHaRLAP工具箱进行电离层射线追踪仿真,并设置地球模型为WGS84模型,各站点位置采用经纬度坐标表示。设置辐射源目标位置为(31.5°N,118.5°E),工作频率为15.01 MHz。参考源1位置为(31.0°N,121.0°E),参考源2位置为(31.0°N,120.0°E),工作频率分别为15.12,15.22 MHz。接收站位置及辐射源目标、参考源到达各接收站信号传输时间如表2所示,模型建立时间为2023-01-15T02:00UT。

        表  2  接收站位置及各信号到达时间
        Tab.  2  Position and arrival time of the receiving stations
        接收站
        序号
        接收站
        经纬度
        目标信号
        传输时间
        /ms
        参考源1
        信号
        传输时间/ms
        参考源2
        信号
        传输时间/ms
        1 (30.5°N, 104.0°E) 5.085 5.842 5.565
        2 (39.5°N, 116.0°E) 3.777 4.070 4.258
        3 (46.0°N, 129.0°E) 6.504 6.361 6.818
        4 (22.5°N, 113.5°E) 4.279 4.533 4.122
        下载: 导出CSV 
        | 显示表格

        在迭代过程给出目标位置估计的同时,信号路径电离层虚高也由IRI模型的射线追踪算法给出。因此即使未能得到目标路径电离层虚高的真实观测结果,本算法仍可在迭代中根据IRI模型推断目标信号路径的电离层虚高,并在迭代优化过程中实现较好的定位效果。采用2个参考源时目标参数估计迭代过程如图9所示,其中40步前目标位置参数估计结果远离参考源,参数估计误差较大;在第40和60步出现阶跃,此时目标与参考源靠近,迭代算法通过电离层虚高相关性方差矩阵 {\boldsymbol{Q}}_{\mathrm{h}}^{} 的影响,自动判断靠近参考源的目标参数估计结果得到更大的似然概率,从而改善了目标电离层虚高估计结果,同时目标位置得到进一步优化。

        图  9  采用2个参考源时目标参数迭代过程
        Fig.  9  Iterative process of target parameters estimation corrected by 2 calibration emitters

        利用所提模型的定位结果如图10所示,采用2个参考源时,对目标的定位结果为(30.996°N,118.346°E),相对真实目标距离57.77 km;仅采用参考源1时,对目标的定位结果为(31.015°N, 117.463°E),相对真实目标距离112.41 km。验证了所提算法使得目标定位精度通过参考源修正得到有效提升,其在更真实的WGS84地球模型和IRI电离层模型条件下仍具有较好的适用性。

        图  10  结合IRI2016模型的射线追踪定位结果
        Fig.  10  Ray tracing localization results combined with IRI2016 model

        本文针对地海表面短波辐射源目标,提出一种目标和参考源联合估计的短波时差定位方法,建立了目标和参考源联合估的时差定位模型,将反射点共用电离层区域的距离相关性引入电离层虚高协方差矩阵中,通过联合估计实现了对目标定位精度的修正。通过推导和仿真所提模型的CRLB,分析了参考源修正定位结果的有效性。进一步给出基于最大似然估计的目标定位Newton迭代解算方法,通过仿真数据验证了所提算法的有效性,实现了良好的定位效果。研究结果对推动短波辐射源定位相关的技术应用具有重要意义。

        1.3节中,式(12)~(15)中多个偏导项的具体表达式如(A1)~(A5)所示。{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}} \mathord{\left/ {\vphantom {{\partial {{\boldsymbol{r}}^o}} {\partial {{\boldsymbol{u}}^o}}}} \right. } {\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}共有{M - 1}行,第i - 1\left(i \in \left[ {2,M} \right]\right){{\partial r_{i1}^{\mathrm{o}}} \mathord{\left/ {\vphantom {{\partial r_{i1}^{\mathrm{o}}} {\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}} \right. } {\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}

        \frac{{\partial {r_{i1}^{\mathrm{o}}}}}{{\partial {\boldsymbol{u}^{\mathrm{o}}}}} = {\left[ \begin{gathered} \frac{{{x_{\boldsymbol{u}}}{{\textit{z}}_{{{\boldsymbol{s}}_i}}} - {{\textit{z}}_{\boldsymbol{u}}}{x_{{{\boldsymbol{s}}_i}}}}}{{{{\textit{z}}_{\boldsymbol{u}}}}} \cdot \frac{{2\left( {{{R}} + {h_i}} \right)}}{{\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_i}} \right\|{r_i}}} - \frac{{{x_{\boldsymbol{u}}}{{\textit{z}}_{{{\boldsymbol{s}}_1}}} - {{\textit{z}}_{\boldsymbol{u}}}{x_{{{\boldsymbol{s}}_1}}}}}{{{{\textit{z}}_{\boldsymbol{u}}}}} \cdot \frac{{2\left( {{{R}} + {h_1}} \right)}}{{\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_1}} \right\|{r_1}}} \\ \frac{{{y_{\boldsymbol{u}}}{{\textit{z}}_{{{\boldsymbol{s}}_i}}} - {{\textit{z}}_{\boldsymbol{u}}}{y_{{{\boldsymbol{s}}_i}}}}}{{{{\textit{z}}_{\boldsymbol{u}}}}} \cdot \frac{{2\left( {{{R}} + {h_i}} \right)}}{{\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_i}} \right\|{r_i}}} - \frac{{{y_{\boldsymbol{u}}}{{\textit{z}}_{{{\boldsymbol{s}}_1}}} - {{\textit{z}}_{\boldsymbol{u}}}{y_{{{\boldsymbol{s}}_1}}}}}{{{{\textit{z}}_{\boldsymbol{u}}}}} \cdot \frac{{2\left( {{{R}} + {h_1}} \right)}}{{\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_1}} \right\|{r_1}}} \\ \end{gathered} \right]^{\text{T}}} (A1)

        {{\partial {{\boldsymbol{r}}^{\mathrm{o}}}} \mathord{\left/ {\vphantom {{\partial {{\boldsymbol{r}}^{\mathrm{o}}}} {\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right. } {\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}共有 {M - 1} 行,第i - 1(i \in \left[ {2,M} \right]){{\partial r_{i1}^{\mathrm{o}}} \mathord{\left/ {\vphantom {{\partial r_{i1}^o} {\partial {{\boldsymbol{h}}^o}}}} \right. } {\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}

        \dfrac{{\partial r_{i1}^{\mathrm{o}}}}{{\partial {\boldsymbol{h}^{\mathrm{o}}}}} = \left[ {\begin{array}{*{20}{c}} { - \dfrac{{4\left( {{{R}} + {h_1}} \right) - 2\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_1}} \right\|}}{{{r_1}}}}&{{{\textit{0}}_{\left( {i - 2} \right) \times 1}}}&{\dfrac{{4\left( {{{R}} + {h_i}} \right) - 2\left\| {{\boldsymbol{u}} + {{\boldsymbol{s}}_i}} \right\|}}{{{r_i}}}}&{{{{\textit{0}}}_{\left( {M - i + MN} \right) \times 1}}} \end{array}} \right] (A2)

        {{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}} \mathord{\left/ {\vphantom {{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}} {\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}} \right. } {\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}共有N\left( {M - 1} \right)行,第\left( {j - 1} \right)\left( {M - 1} \right) + \left( {i - 1} \right)(j \in \left[ {1,N} \right]i \in \left[ {2,M} \right])行{{\partial \tilde r_{i1,j}^{\mathrm{o}}} \mathord{\left/ {\vphantom {{\partial \tilde r_{i1,j}^o} {\partial {{\boldsymbol{h}}^o}}}} \right. } {\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}

        \dfrac{{\partial \tilde r_{i1,j}^{\mathrm{o}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}} = \left[ {\begin{array}{*{20}{c}} {{{{\textit{0}}}_{jM \times 1}}}&{ - \dfrac{{4\left( {R + {h_{1,j}}} \right) - 2\left\| {{{\boldsymbol{c}}_j} + {{\boldsymbol{s}}_1}} \right\|}}{{{{\tilde r}_{1,j}}}}}&{{{{\textit{0}}}_{\left( {i - 2} \right) \times 1}}}&{\dfrac{{4\left( {R + {h_{i,j}}} \right) - 2\left\| {{{\boldsymbol{c}}_j} + {{\boldsymbol{s}}_i}} \right\|}}{{{{\tilde r}_{i,j}}}}}&{{{{\textit{0}}}_{\left( {M - i + M\left( {N - j} \right)} \right) \times 1}}} \end{array}} \right] (A3)
        \dfrac{\partial {\left[{{\boldsymbol{Q}}}_{{\mathrm{h}}}\right]}_{p=mM+i,q=nM+j}}{\partial {x}_{\boldsymbol{u}}}=\left\{\begin{array}{cc}{\left[{{\boldsymbol{Q}}}_{{\mathrm{h}}}\right]}_{pq}\cdot \dfrac{2\left({x}_{{\boldsymbol{c}}_{m}}+{x}_{{\boldsymbol{s}}_{i}}-{x}_{{\boldsymbol{s}}_{j}}\right){\textit{z}}_{\boldsymbol{u}}-2\left({\textit{z}}_{{\boldsymbol{c}}_{m}}+{\textit{z}}_{{\boldsymbol{s}}_{i}}-{\textit{z}}_{{\boldsymbol{s}}_{j}}\right){x}_{\boldsymbol{u}}}{{\textit{z}}_{\boldsymbol{u}}{R}_{\mathrm{cov}}^{2}{}}& m > 0,n=0\\ {\left[{{\boldsymbol{Q}}}_{{\mathrm{h}}}\right]}_{pq}\cdot \dfrac{-2\left(-{x}_{{\boldsymbol{c}}_{n}}+{x}_{{\boldsymbol{s}}_{i}}-{x}_{{\boldsymbol{s}}_{j}}\right){\textit{z}}_{\boldsymbol{u}}+2\left(-{\textit{z}}_{{\boldsymbol{c}}_{n}}+{\textit{z}}_{{\boldsymbol{s}}_{i}}-{\textit{z}}_{{\boldsymbol{s}}_{j}}\right){x}_{\boldsymbol{u}}}{{\textit{z}}_{\boldsymbol{u}}{R}_{\mathrm{cov}}^{2}{}}& n > 0,m=0\\ 0& 其它\end{array}\right. (A4)
        \dfrac{\partial {\left[{{\boldsymbol{Q}}}_{{\mathrm{h}}}\right]}_{p=mM+i,q=nM+j}}{\partial {y}_{\boldsymbol{u}}}=\left\{\begin{array}{cc}{\left[{{\boldsymbol{Q}}}_{{\mathrm{h}}}\right]}_{pq}\cdot \dfrac{2\left({y}_{{\boldsymbol{c}}_{m}}+{y}_{{\boldsymbol{s}}_{i}}-{y}_{{\boldsymbol{s}}_{j}}\right){\textit{z}}_{\boldsymbol{u}}-2\left({\textit{z}}_{{\boldsymbol{c}}_{m}}+{\textit{z}}_{{\boldsymbol{s}}_{i}}-{\textit{z}}_{{\boldsymbol{s}}_{j}}\right){y}_{\boldsymbol{u}}}{{\textit{z}}_{\boldsymbol{u}}{R}_{\mathrm{cov}}^{2}{}}& m > 0,n=0\\ {\left[{{\boldsymbol{Q}}}_{{\mathrm{h}}}\right]}_{pq}\cdot \dfrac{-2\left(-{y}_{{\boldsymbol{c}}_{n}}+{y}_{{\boldsymbol{s}}_{i}}-{y}_{{\boldsymbol{s}}_{j}}\right){\textit{z}}_{\boldsymbol{u}}+2\left(-{\textit{z}}_{{\boldsymbol{c}}_{n}}+{\textit{z}}_{{\boldsymbol{s}}_{i}}-{\textit{z}}_{{\boldsymbol{s}}_{j}}\right){y}_{\boldsymbol{u}}}{{\textit{z}}_{\boldsymbol{u}}{R}_{\mathrm{cov}}^{2}{}}& n > 0,m=0\\ 0& 其它\end{array} \right. (A5)

        1.6节中,{\boldsymbol{b}}表达式为

        {\boldsymbol{b}} = \left[ \begin{gathered} \frac{{\partial \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right)}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}} \\ \frac{{\partial \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right)}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}} \\ \end{gathered} \right] (B1)
        \begin{split} \dfrac{{\partial \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right)}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}} = & \dfrac{\partial }{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}\left[ { - \dfrac{1}{2}{{\left( {{\boldsymbol{r}} - {{\boldsymbol{r}}^{\mathrm{o}}}} \right)}^{\text{T}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {{\boldsymbol{r}} - {{\boldsymbol{r}}^{\mathrm{o}}}} \right)} \right] + \dfrac{\partial }{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}\left[ { - \dfrac{1}{2}\ln \left| {{{\boldsymbol{Q}}_{\mathrm{h}}}} \right| - \dfrac{1}{2}{{\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)}^{\text{T}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)} \right] \\ = & \dfrac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{u}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {{\boldsymbol{r}} - {{\boldsymbol{r}}^{\mathrm{o}}}} \right) - \dfrac{1}{2}\left[ {\begin{array}{*{20}{c}} {{\text{tr}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {x_{\boldsymbol{u}}}}}} \right)} \\ {{\text{tr}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {y_{\boldsymbol{u}}}}}} \right)} \end{array}} \right] + \dfrac{1}{2}\left[ {\begin{array}{*{20}{c}} {{{\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)}^{\text{T}}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {x_{\boldsymbol{u}}}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}} \right)\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)} \\ {{{\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)}^{\text{T}}}\left( {{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\dfrac{{\partial {{\boldsymbol{Q}}_{\mathrm{h}}}}}{{\partial {y_{\boldsymbol{u}}}}}{\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}} \right)\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right)} \end{array}} \right] \end{split} (B2)
        \frac{{\partial \ln p\left( {{\boldsymbol{r}},{\boldsymbol{\tilde r}},{\boldsymbol{h}};{{\boldsymbol{u}}^{\mathrm{o}}},{{\boldsymbol{h}}^{\mathrm{o}}}} \right)}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}} = \frac{{\partial {{\boldsymbol{r}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\text{α}} ^{ - 1}\left( {{\boldsymbol{r}} - {{\boldsymbol{r}}^{\mathrm{o}}}} \right) + \frac{{\partial {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}}}{{\partial {{\boldsymbol{h}}^{\mathrm{o}}}}}{\boldsymbol{Q}}_{\boldsymbol{c}}^{ - 1}\left( {{\boldsymbol{\tilde r}} - {{{\boldsymbol{\tilde r}}}^{\mathrm{o}}}} \right) + {\boldsymbol{Q}}_{\mathrm{h}}^{ - 1}\left( {{\boldsymbol{h}} - {{\boldsymbol{h}}^{\mathrm{o}}}} \right) (B3)
      • 图  1   含有参考源的时差定位场景

        Fig.  1   Scenario of TDOA localization with calibration emitter

        图  2   电离层虚高相关系数随反射点相对距离变化

        Fig.  2   Ionospheric virtual height correlation coefficient vs. relative distance of reflection points

        图  3   基于Armijo直线搜索Newton法的最大似然估计流程

        Fig.  3   Flowchart of maximum likelihood estimation based on Armijo’s linear search Newton method

        图  4   电离层虚高方差 \sigma _{\mathrm{h}}^2 对定位精度的影响

        Fig.  4   Influence of ionospheric virtual height variance \sigma _{\mathrm{h}}^2 on localization accuracy

        图  5   时间测量方差\sigma _{\mathrm{t}}^2对定位精度的影响

        Fig.  5   Influence of time measurement variance \sigma _{\mathrm{t}}^2 on localization accuracy

        图  6   接收站和参考源对区域目标定位精度影响

        Fig.  6   Influence of the distribution of receiving stations and calibration emitters on localization accuracy

        图  7   算法场景1条件下所提估计方法单次仿真迭代过程

        Fig.  7   Iterative process of the proposed estimation method for algorithm scenario 1 in a single simulation

        图  8   算法场景2条件下所提估计方法单次仿真迭代过程

        Fig.  8   Iterative process of the proposed estimation method for algorithm scenario 2 in a single simulation

        图  9   采用2个参考源时目标参数迭代过程

        Fig.  9   Iterative process of target parameters estimation corrected by 2 calibration emitters

        图  10   结合IRI2016模型的射线追踪定位结果

        Fig.  10   Ray tracing localization results combined with IRI2016 model

        表  1   接收站和参考源以及目标位置坐标

        Tab.  1   Coordinate position of receiving stations, calibration emitters and their target location

        序号 接收站/km 参考源/km 目标位置/km
        1 [0, 0] [500 \sqrt2 ,500 \sqrt2] [900, 900]
        2 [1000, 0] [750 \sqrt2, 750 \sqrt2] [750, 750]
        3 [0, 1000] [500 \sqrt3, 500] [710, 710]
        4 [500, 500] [500, 500 \sqrt3] [1200, 1200]
        5 [200, 500] [ 750 \sqrt3, 750] [1550, 800]
        6 [500, 200] [750, 750 \sqrt3] [950, 800]
        下载: 导出CSV

        表  2   接收站位置及各信号到达时间

        Tab.  2   Position and arrival time of the receiving stations

        接收站
        序号
        接收站
        经纬度
        目标信号
        传输时间
        /ms
        参考源1
        信号
        传输时间/ms
        参考源2
        信号
        传输时间/ms
        1 (30.5°N, 104.0°E) 5.085 5.842 5.565
        2 (39.5°N, 116.0°E) 3.777 4.070 4.258
        3 (46.0°N, 129.0°E) 6.504 6.361 6.818
        4 (22.5°N, 113.5°E) 4.279 4.533 4.122
        下载: 导出CSV
      • [1] 张旭辉,姜春华,刘桐辛,等. 电离层虚高对超视距雷达多站联合定位精度的影响[J]. 电波科学学报,2022,37(5):761-767.

        ZHANG X H,JIANG C H,LIU T X,et al. Effect of the ionospheric virtual height on the joint positioning accuracy of multi-station over-the-horizon radar system[J]. Chinese journal of radio science,2022,37(5):761-767. ( in Chinese

        [2]

        YANG L,GAO H,LING Y,et al. Localization method of wide-area distribution multistatic sky-wave over-the-horizon radar[J]. IEEE geoscience and remote sensing letters,2022(19). DOI: 10.1109/LGRS.2020.3018322

        [3]

        HUANG S,PUN Y M,SO M C,et al. A provably convergent projected gradient-type algorithm for TDOA-based geolocation under the quasi-parabolic ionosphere model[J]. IEEE signal processing letters,2020(27):1335-1339. doi: 10.1109/LSP.2020.3010676

        [4]

        QI H,WU X,JIA L. Semidefinite programming for unified TDOA-based localization under unknown propagation speed[J]. IEEE communications letters,2020,24(9):1971-1975. doi: 10.1109/LCOMM.2020.2996970

        [5] 刘桐辛,周晨,杨国斌,等. 基于多站观测的短波电离层信道一致性与相参性研究[J]. 电波科学学报,2022,37(3):364-371. doi: 10.12265/j.cjors.2021127

        LIU T X,ZHOU C,YANG G B,et al. Research on the consistency and coherence of the ionospheric short-wave transmission channel based on multi-station sounding[J]. Chinese journal of radio science,2022,37(3):364-371. (in Chinese) doi: 10.12265/j.cjors.2021127

        [6]

        BOURGEOIS D,MORISSEAU C,FLECHEUX M. Over-the-horizon radar target tracking using multi-quasi-parabolic ionospheric modelling[J]. IET proceedings radar sonar navigation,2006,153(5):409-416. doi: 10.1048/IP-RSN.20050100

        [7]

        HO K C,YANG L. On the use of a calibration emitter for source localization in the presence of sensor position uncertainty[J]. IEEE transactions on signal processing,2008,56(12):5758-5772. doi: 10.1109/TSP.2008.929870

        [8] 秦兆涛,王俊,陶磊岩,等. 基于标校源辅助的不相交多目标到达时差定位[J]. 北京航空航天大学学报,2015,44(5):1026-1036. doi: 10.13700/j.bh.1001-5965.2017.0370

        QIN Z T,WANG J,TAO L Y,et al. TDCA localization of multiple disjoint sources based on a calibration emitter[J]. Journal of Beijing University of Aeronautics and Astronautics,2015,44(5):1026-1036. (in Chinese) doi: 10.13700/j.bh.1001-5965.2017.0370

        [9] 李深,周晨,王君明,等. 基于经验电离层模型的短波时差定位理论分析[J]. 系统工程与电子技术,2023,45(7):1911-1919.

        LI S,ZHOU C,WANG M J,et al. Theoretical analysis of shortwave TDOA geolocation based on empirical ionospheric model[J]. Systems engineering arid electronics,2023,45(7):1911-1919. (in Chinese)

        [10]

        XIA N,XING B H. A direct localization method for HF source geolocation and experimental results[J]. IEEE antennas and wireless propagation letters,2021,20(5):725-732. doi: 10.1109/LAWP.2021.3061507

        [11] 张铁男. 基于电离层反射信号的多站短波时差定位技术研究[D]. 哈尔滨:哈尔滨工业大学,2019.

        ZHANG T N. Study of multi-station short-wave time-difference-of-arrival localization technique based on ionosphere-layer reflected signal[D]. Harbin:Harbin Institute of Technology,2019. (in Chinese)

        [12] 王鼎,尹洁听,朱中梁. 针对超视距短波辐射源的测角与测时差协同定位方法[J]. 中国科学:信息科学,2022,52:1942-1973. doi: 10.1360/SSI-2021-0331

        WANG D,YIN J T,ZHU Z L. Novel cooperative localization method of over-the-horizon shortwave emitters based on direction-of-arrival and time-difference-of-arrival measurements[J]. Scientia sinica information,2022,52:1942-1973. (in Chinese) doi: 10.1360/SSI-2021-0331

        [13]

        LE Y,HO K C. Alleviating sensor position error in source localization using calibration emitters at inaccurate locations[J]. IEEE transactions on signal processing,2010,58(1):67-83. doi: 10.1109/TSP.2009.2028947

        [14] KAY S M. 统计信号处理基础估计与检测理论[M]. 北京:电子工业出版社,2014:55-61.
        [15] 张显达. 矩阵分析与应用 [M]. 2版. 北京:清华大学出版社,2004:267-274.
      • 期刊类型引用(0)

        其他类型引用(1)

      图(10)  /  表(2)
      计量
      • 文章访问数:  187
      • HTML全文浏览量:  24
      • PDF下载量:  55
      • 被引次数: 1
      出版历程
      • 收稿日期:  2023-06-28
      • 录用日期:  2023-11-27
      • 网络出版日期:  2023-11-27
      • 刊出日期:  2024-06-29

      目录

      /

      返回文章
      返回