• 中国科学引文数据库(CSCD)
  • 中文科技期刊数据库
  • 中国核心期刊(遴选)数据库
  • 日本科学技术振兴机构数据库(JST)
  • 中国学术期刊(网络版)(CNKI)
  • 中国学术期刊综合评价数据库(CAJCED)
  • 中国超星期刊域出版平台

卫星导航干扰源无人机定位观测路径规划方法

姚胤博, 许睿, 徐悦, 刘睿, 陈奇东

姚胤博, 许睿, 徐悦, 刘睿, 陈奇东. 卫星导航干扰源无人机定位观测路径规划方法[J]. 全球定位系统, 2024, 49(6): 72-83. DOI: 10.12265/j.gnss.2024097
引用本文: 姚胤博, 许睿, 徐悦, 刘睿, 陈奇东. 卫星导航干扰源无人机定位观测路径规划方法[J]. 全球定位系统, 2024, 49(6): 72-83. DOI: 10.12265/j.gnss.2024097
YAO Yinbo, XU Rui, XU Yue, LIU Rui, CHEN Qidong. Method for configuring observation point paths for UAV satellite navigation interference source localization[J]. GNSS World of China, 2024, 49(6): 72-83. DOI: 10.12265/j.gnss.2024097
Citation: YAO Yinbo, XU Rui, XU Yue, LIU Rui, CHEN Qidong. Method for configuring observation point paths for UAV satellite navigation interference source localization[J]. GNSS World of China, 2024, 49(6): 72-83. DOI: 10.12265/j.gnss.2024097

卫星导航干扰源无人机定位观测路径规划方法

详细信息
    作者简介:

    姚胤博: (2000—),男,硕士,研究方向为卫星导航干扰源定位. E-mail:sz2203108@nuaa.edu.cn

    许睿: (1983—),女,副教授,研究方向为卫星导航抗干扰/抗欺骗研究. E-mail:ruixu@nuaa.edu.cn

    徐悦: (1994—),女,硕士,工程师,研究方向为卫星导航对抗等. E-mail:951061712@qq.com

    刘睿: (1989—),男,博士,高级工程师,研究方向为卫星导航对抗等. E-mail:15375516858@163.com

    陈奇东: (1980—),男,博士,正高级工程师,研究方向为卫星导航对抗、电波环境应用等. E-mail:chenqd0929@126.com

    通信作者:

    许睿 E-mail:ruixu@nuaa.edu.cn

  • 中图分类号: P228.4

Method for configuring observation point paths for UAV satellite navigation interference source localization

  • 摘要:

    本文面向单无人机干扰源搜寻系统中观测点配置需求,针对干扰源位置未知、测角设备存在盲区以及精度因子(dilution of precision,DOP)计算复杂问题,比较了不同DOP与误差分布的吻合性,提出合适的观测点配置方法辅助完成干扰源搜寻的工作. 交叉定位中观测点与干扰源的相对分布对定位结果存在影响,因此引入DOP来描述分布的优劣. 将合适的DOP作为规划指标,提出定步长方法、不定步方法、均匀方法提供监测观测点配置,对不同初始点分布时各方法的配置效果做出比较. 并针对上述3种方法的缺陷提出改进的配置算法,均通过仿真实验进行验证. 所提出的改进配置方法符合单无人机监测过程中顺序增加观测点的实际情况,即使在初始角度分布极差时也可达到与理想精度较为接近的精度,可满足监测定位工作的需求.

    Abstract:

    This paper addresses the optimization of observation point configuration in a single UAV interference source search system, focusing on unknown interference source positions, goniometric equipment blind spots, and complex calculations of dilution of precision. It compares various dilution of precision with error distribution and proposes a method to enhance the search process. By using dilution of precision to evaluate distribution effectiveness, the paper presents three methods: fixed-step, variable-step, and uniform. These methods are compared under different initial distributions. Based on the shortcomings of the three methods above, an improved algorithm is proposed and validated through simulations. This algorithm supports the sequential addition of observation points, achieving near-ideal accuracy even with poor initial angle distributions, thereby meeting monitoring and localization requirements.

  • GNSS作为现今一种重要导航技术,在各领域中发挥着重要的作用. 但GNSS容易受到电磁波干扰[1-2],为了能够有效应对干扰攻击,快速处理干扰源,相关研究已提出了针对抗干扰源的定位系统.

    对干扰源的定位本质上是基于各种观测信息对干扰源的位置进行估计[3-4],包括到达时间定位法(location by difference of times, TOA)、到达角定位法(angle of arrival, AOA)、接收信号强度(received signal strength indicator, RSS)等方法[5]. 由于实际进行观测时必然存在误差,观测值的误差会传递到定位的结果,该过程受观测点与干扰源之间的几何分布影响[6-7],合理的观测点设置可以改善定位精度[8],例如文献[9-10]依据AOA定位的精度因子推算最优观测点分布. 但目前的观测点配置方法基于全局配置[11-12],需要提前预知干扰源位置以及可达观测点数量,临时增加观测点会涉及到已有观测点的修改,难以满足实际搜寻的需要.

    本文以一种采用AOA定位的无人机干扰监测系统为背景,对整个干扰源监测过程中观测点位置对最终定位精度产生的影响进行分析,比较了几种用于衡量观测点位置分布优劣的指标. 将最合适的指标作为依据,提出符合实际限制条件的观测点规划方法,不改变已有点的位置,仅配置新的观测点,最终达到需要的定位精度水平.

    AOA定位法是通过阵列天线、空间谱等测向技术[13]获得干扰源所发出信号到达观测点的入射角度,结合观测点自身定位模块提供的位置信息进行干扰源的定位方法[14]. 无人机干扰源监测系统进行干扰源搜寻的工作中,无人机搭载测向载荷获得干扰信号方向信息,定位模块提供无人机所在的位置信息. 方向信息和位置信息共同组成观测数据,无人机抵达不同观测点进行测向后,获得多组观测数据组,利用AOA算法进行干扰源定位. 由于测向模块的测向天线存在限制,监测系统的测向结果并没有给出俯仰角,因此本研究主要在二维平面内开展,监测系统测向定位原理如图1所示.

    图  1  测向定位原理图

    图1所示的$xOy$平面为研究场景,采用的定位算法原理如下:假设定位目标的干扰源位置为$({x_r},{y_r})$,各观测点的坐标为$\left( {{x_i},{y_i}} \right)$,$i = 1,2 \cdots ,N$,将观测点$i$与干扰源连线与$x$轴正半轴之间的夹角记作${\theta _i}$,并将观测点$i$和观测点$j$对应的角${\theta _i}$${\theta _j}$的差值记作${\theta _{ij}}$,将观测点$i$到干扰源之间的距离记作${d_i}$.

    根据观测点与干扰源之间的几何分布关系,可以得到下式:

    $$ \begin{cases} {{x_r} - {x_i} = {d_i}\cos\; {\theta _i}} \\ {{y_r} - {y_i} = {d_i}\sin\; {\theta _i}} \end{cases} $$ (1)

    式中,$ d_i $为两者之间的距离.

    $$ \frac{{{y_r} - {y_i}}}{{{x_r} - {x_i}}} = \frac{{\sin\; {\theta _i}}}{{\cos\; {\theta _i}}} $$ (2)

    将所有观测点的整理后得到如下的方程组

    $$ \begin{array}{*{20}{l}}\qquad\qquad\qquad\qquad\quad\boldsymbol{AX}=\boldsymbol{B} \\ \boldsymbol{A}=\begin{bmatrix}\mathrm{cos}\; \theta_1 & -\mathrm{sin}\; \theta_1 \\ \mathrm{cos}\; \theta_2 & -\mathrm{sin}\; \theta_2 \\ \vdots & \vdots \\ \theta_i & -\mathrm{sin}\; \theta_i\end{bmatrix},\; \boldsymbol{B}=\begin{bmatrix}\mathrm{cos}\; \theta_1y_1-\mathrm{sin}\; \theta_1x_1 \\ \mathrm{cos}\; \theta_2y_2-\mathrm{sin}\; \theta_2x_2 \\ \vdots \\ \mathrm{cos}\; \theta_iy_i-\mathrm{sin}\; \theta_ix_i\end{bmatrix} \\ \boldsymbol{X}=\begin{bmatrix}y_r \\ x_r\end{bmatrix}\end{array} $$ (3)

    可得到干扰源在$ x、y $平面上坐标的最小二乘解为

    $$ {\boldsymbol{X}} = {({{\boldsymbol{A}}^{\text{T}}}{\boldsymbol{A}})^{ - 1}}{{\boldsymbol{A}}^{\text{T}}}{\boldsymbol{B}}. $$ (4)

    在测量值没有误差的前提下,上式的结果可以较为准确的反映干扰源的位置. 但在实际监测过程中,测向误差不可避免,无人机自身的位置测量也存在误差,导致干扰源定位误差同时受测向误差与位置误差影响. 同时,观测点与干扰源之间的相对分布也会影响干扰源的定位精度,极端情况如2个观测点与干扰源接近共线,即夹角差${\theta _{ij}}$约为180°时,较小的测角误差会带来极大的定位误差.

    为了衡量观测点与干扰源的相对几何构型对定位精度的影响,以便展开后续的观测点配置工作,本文将精度因子(dilution of precision, DOP)的概念引入干扰源定位的过程中. DOP是用于描述测量误差放大系数的一个量,常用于GNSS定位中,用于描述不同可见星分布情况下定位精度误差的放大量. DOP越小,误差放大系数越低,定位精度越高. 结合AOA定位原理,DOP可以采用3种不同的计算方法.

    位/距无关DOP,记为$ {\text{DO}}{{\text{P}}_{\text{1}}} $,与卫星导航中的DOP类似,来源于定位方程中的几何矩阵. AOA定位过程中所用的传递矩阵${\boldsymbol{A}}$如式(3)所示.

    $ \boldsymbol{H}=\boldsymbol{A}^{\text{T}}\boldsymbol{A} $,则有

    $$ {\boldsymbol{H}} = {\begin{bmatrix} { \displaystyle\sum\limits_{i = 1}^N {{\cos }^2}{\theta _i}}&{ - \displaystyle\sum\limits_{i = 1}^N \sin \;{\theta _i}\cos\; {\theta _i}} \\ { - \displaystyle\sum\limits_{i = 1}^N \sin\; {\theta _i}\cos\; {\theta _i}}& { \displaystyle\sum\limits_{i = 1}^N {{\sin }^2}{\theta _i}} \end{bmatrix}} $$ (5)

    ${{\boldsymbol{H}}^{ - 1}}$记作$ {\text{DO}}{{\text{P}}_{\text{1}}} $的特征矩阵${{\boldsymbol{P}}_1}$,整理得到$ {\text{DO}}{{\text{P}}_{\text{1}}} $表达式

    $$ {\text{DO}}{{\text{P}}_{\text{1}}} = \sqrt {{\mathrm{tr}}\left( {{{\boldsymbol{P}}_1}} \right)} = \sqrt {\frac{N}{{\displaystyle\sum\limits_{i = 1}^{N - 1} \displaystyle \sum\limits_{j = i + 1}^N {{\sin }^2}{\theta _{ij}}}}} $$ (6)

    式(6)中,$N$表示观测点的个数. 该式展示出$ {\text{DO}}{{\text{P}}_{\text{1}}} $仅与观测点数量和顺次观测点夹角${\theta _{ij}}$,与观测点位置以及观测点与干扰源距离无关.

    $ {\text{DO}}{{\text{P}}_{\text{1}}} $不能反映出观测点位置误差或观测点到干扰源的距离误差对定位结果的影响. 将多个干扰源的式(1)组合后,进一步整理为

    $$ \left\{ \begin{gathered} (x - {x_1})\tan \;{\theta _1} = y - {y_1} \\ (x - {x_2})\tan \;{\theta _2} = y - {y_2} \\ \cdots \\ (x - {x_i})\tan \;{\theta _i} = y - {y_i} \\ \end{gathered} \right. $$ (7)

    对(7)求导,整理可得

    $$ \left\{ \begin{gathered} - \tan \;{\theta _1}{\text{d}}x + {\text{d}}y = (x - {x_1}){\sec ^2}{\theta _1}{\text{d}}{\theta _1} \\ - \tan \;{\theta _2}{\text{d}}x + {\text{d}}y = (x - {x_2}){\sec ^2}{\theta _2}{\text{d}}{\theta _2} \\ \cdots \\ - \tan\; {\theta _i}{\text{d}}x + {\text{d}}y = (x - {x_i}){\sec ^2}{\theta _i}{\text{d}}{\theta _i} \\ \end{gathered} \right. $$ (8)

    经整理可得以下方程

    $$ \begin{array}{*{20}{l}}\qquad\qquad\qquad\quad\bar{\boldsymbol{A}}\text{d}\boldsymbol{X}=\bar{\boldsymbol{B}} \\ \bar{\boldsymbol{A}}=\begin{bmatrix}-\tan\; \theta_1 & 1 \\ -\tan\; \theta_2 & 1 \\ \vdots & \vdots \\ -\tan\; \theta_i & 1\end{bmatrix},\ \bar{\boldsymbol{B}}=\begin{bmatrix}\left(x_r-x_1\right)\sec^2\theta_1\text{d}\theta_1 \\ \left(x_r-x_2\right)\sec^2\theta_2\text{d}\theta_2 \\ \vdots \\ \left(x_r-x_i\right)\sec^2\theta_i\text{d}\theta_i\end{bmatrix} \\ \text{d}\boldsymbol{X}=\begin{bmatrix}{{\mathrm{d}}}x_r \\ {{\mathrm{d}}}y_r\end{bmatrix} \\ \end{array} $$ (9)

    距离相关$ {\text{DO}}{{\text{P}}_{\text{2}}} $的特征矩阵$ {{\boldsymbol{P}}_2} $如下:

    $$ \begin{gathered}\text{d}\boldsymbol{X}=(\bar{\boldsymbol{A}}^{\text{T}}\bar{\boldsymbol{A}})^{-1}\bar{\boldsymbol{A}}^{\text{T}}\boldsymbol{B}\text{,} \\ \text{DO}\text{P}_{\text{2}}=\sqrt{\mathrm{tr}\left(\boldsymbol{P}_2\right)}=\sqrt{\mathrm{tr}\left(\boldsymbol{E}\left(\text{d}\boldsymbol{X}\text{d}\boldsymbol{X}^{\text{T}}\right)\right)} \\ \end{gathered} $$ (10)

    假设各观测点所观测到的角度噪声为方差一致的高斯白噪声,则$ {\text{DO}}{{\text{P}}_{\text{2}}} $的表达式如下:

    $$ {\text{DO}}{{\text{P}}_2} = \sqrt {\frac{N}{{\displaystyle\sum\limits_{i = 1}^{N - 1} \displaystyle\sum\limits_{j = i + 1}^N {{\sin }^2}{\theta _{ij}}}}} \sqrt {\displaystyle\sum\limits_{i = 1}^N d_i^2} $$ (11)

    式中:${\sigma ^2}$是角度测量值的噪声方差.

    观察上式可以发现,$ {\text{DO}}{{\text{P}}_{\text{2}}} $不仅与观测点和干扰源间的角度差有关,也与观测点与干扰源间的距离有关,距离越小$ {\text{DO}}{{\text{P}}_{\text{2}}} $越小. 比较式(6)与式(12),发现$ {\text{DO}}{{\text{P}}_{\text{2}}} $涉及夹角${\theta _{ij}}$部分与$ {\text{DO}}{{\text{P}}_{\text{1}}} $相同,是在$ {\text{DO}}{{\text{P}}_{\text{1}}} $基础上乘以距离相关项,即

    $$ {\text{DO}}{{\text{P}}_2} = {\text{DO}}{{\text{P}}_1}\sqrt {\mathop \sum \limits_{i = 1}^N d_i^2} $$ (12)

    当所有观测点与干扰源距离相等,即在以干扰源为圆心的圆周上,则距离变成常数,此时$ {\text{DO}}{{\text{P}}_2} $的变化趋势则与$ {\text{DO}}{{\text{P}}_{\text{1}}} $一致.

    $ {\text{DO}}{{\text{P}}_{\text{1}}} $中,已经求过最小二乘算法的特征矩阵,考虑到实际测量存在的测量误差,设测量误差矩阵为${\boldsymbol{e}}$,式(3)可变为以下的形式:

    $$ {\boldsymbol{A}}{{\boldsymbol{X}}_T} = {\boldsymbol{B}} + {\boldsymbol{e}} $$ (13)

    $\theta = {\left( {{\theta _1},{\theta _2}, \ldots ,{\theta _n}} \right)^{\mathrm{T}}}$,在上式中对$\theta $求偏导

    $$ \frac{{{\text{d}}{\boldsymbol{B}}}}{{{\text{d}}{\boldsymbol{\theta}} }} = {\text{diag}}\left( {\frac{{{\text{d}}{b_1}}}{{{\text{d}}{\theta _1}}},\frac{{{\text{d}}{b_2}}}{{{\text{d}}{\theta _2}}}, \cdots ,\frac{{{\text{d}}{b_n}}}{{{\text{d}}{\theta _n}}}} \right) $$ (14)

    其中$\dfrac{{{\text{d}}{b_i}}}{{{\text{d}}{\theta _i}}} = {x_i}{\text{cos}}\;{\theta _i} + {y_i}{\text{sin}}\;{\theta _i}$.

    ${\boldsymbol{e}}$的协方差矩阵${{\boldsymbol{P}}_{\boldsymbol{e}}}$表达式为

    $$ {{\boldsymbol{P}}_{\boldsymbol{e}}} = {\boldsymbol{E}}\left( {{\boldsymbol{e}}{{\boldsymbol{e}}^{\text{T}}}} \right) = {\text{diag}}\left( {{g_1},{g_2}, \cdots ,{g_n}} \right){\sigma ^2} $$ (15)

    其中$ g_i=\left(x_i\text{cos}\ \theta_i+y_i\text{sin}\ \theta_i\right)^2 $.

    可以得到$ {\text{DO}}{{\text{P}}_{\text{3}}} $的特征矩阵${{\boldsymbol{P}}_3}$

    $$ \begin{array}{l}\text{Δ}{\boldsymbol{X}}={\boldsymbol{X}}-{{\boldsymbol{X}}}_{{\boldsymbol{T}}}=({{\boldsymbol{A}}}^{{\mathrm{T}}} {\boldsymbol{A}})^{-1}{{\boldsymbol{A}}}^{{\mathrm{T}}}\left(-{\boldsymbol{e}}\right)\\ {{\boldsymbol{P}}}_{3}={\boldsymbol{E}}\left(\text{Δ}{\boldsymbol{X}}\text{Δ}{{\boldsymbol{X}}}^{\text{T}}\right)={{\boldsymbol{H}}}^{-1}{{\boldsymbol{A}}}^{\text{T}}{{\boldsymbol{P}}}_{{\boldsymbol{e}}}({{{\boldsymbol{H}}}^{-1}{\boldsymbol{A}}^{\rm {T}}})^{\rm {T}}\end{array} $$ (16)

    $ {\text{DO}}{{\text{P}}_{\text{3}}} $

    $$ {\text{DO}}{{\text{P}}_{\text{3}}} = \sqrt {{\mathrm{tr}}\left( {{{\boldsymbol{P}}_3}} \right)} = \frac{{\sqrt {{l_{11}} + {l_{22}}} }}{{\displaystyle\sum\limits_{i = 1}^{N - 1} \displaystyle\sum\limits_{j = i + 1}^N {{\sin }^2}{\theta _{ij}}}} $$ (17)

    其中

    $$ \begin{array}{*{20}{l}} {{l_{11}} = \displaystyle\sum \limits_{i = 1}^N {{\left( {{f_{11}}\sin \;{\theta _i} - {f_{12}}\cos \;{\theta _i}} \right)}^2}{g_i}} \\ {{l_{22}} = \displaystyle\sum \limits_{i = 1}^N {{\left( {{f_{21}}\sin \;{\theta _i} - {f_{22}}\cos \;{\theta _i}} \right)}^2}{g_i}} \end{array} $$ (18)
    $$ \begin{gathered}f_{11}=\sum\limits_{i=1}^N\cos^2\theta_i \\ f_{22}=\sum\limits_{i=1}^N\sin^2\theta_i \\ f_{12}=f_{21}=\sum\limits_{i=1}^N\sin\ \theta_i\cos\ \theta_i \\ \end{gathered} $$ (19)

    $ {\text{DO}}{{\text{P}}_{\text{3}}} $的表达式最为复杂,可以反映出观测点到干扰源距离对定位精度的影响,但计算时无需干扰源到观测点间的距离,仅需要角度和位置信息(${g_i}$).

    假设在以干扰源为原点的坐标轴上存在2个观测点,并且两者到干扰源的距离相等(均为$r$),此时DOP2与DOP1的关系如下:

    $$ {\text{DO}}{{\text{P}}_2} = {\text{DO}}{{\text{P}}_1}\sqrt 2 r $$ (20)

    对于DOP3,将上述情况带入式(16)可得到

    $$ {\text{DO}}{{\text{P}}_{\text{3}}} = \frac{{\sqrt {{l_{11}} + {l_{22}}} }}{{{{\sin }^2}{\theta _{12}}}} $$ (21)

    其中分母部分相关项代入距离$r$后可简化为

    $$ g_i=\left(x_i\mathrm{cos\ }\theta_i+y_i\mathrm{sin}\ \theta_i\right)^2=\left(\frac{x_i}{\mathrm{cos}^2\theta_i}+\frac{y_i}{\mathrm{sin}^2\theta_i}\right)^2=r^2 $$ (22)
    $$ \begin{split} {l_{11}} =& \mathop \sum \limits_{i = 1}^N {\left( {{f_{11}}\sin\; {\theta _i} - {f_{12}}\cos\; {\theta _i}} \right)^2}{g_i} \\ =& ( ( {{{\cos }^2}{\theta _1} + {{\cos }^2}{\theta _2}} )\sin\; {\theta _1}\\ &-( {\sin\; {\theta _1}\cos\; {\theta _2} + \sin\; {\theta _2}\cos\; {\theta _1}} )\cos\; {\theta _2} )^{\text{2}}{r^2} \\ {l_{22}} =& \mathop \sum \limits_{i = 1}^N {\left( {{f_{21}}\sin\; {\theta _i} - {f_{22}}\cos\; {\theta _i}} \right)^2}{g_i} \\ =& ( ( {\sin\; {\theta _1}\cos\; {\theta _2} + \sin\; {\theta _2}\cos\; {\theta _1}} )\sin\; {\theta _1} \\&- ( {{{\sin }^2}{\theta _1} + {{\sin }^2}{\theta _2}} )\cos\; {\theta _2} )^2{r^2} \end{split} $$ (23)

    最终,式(21)分母部分可简化为

    $$ {l_{11}} + {l_{22}} = 2{\sin ^2}{\theta _{12}}{r^2}, $$ (24)

    将式(24)代入式(21),可得$ {\text{DO}}{{\text{P}}_{\text{3}}} $表达式

    $$ {\text{DO}}{{\text{P}}_{\text{3}}} = \frac{{\sqrt {2{{\sin }^2}{\theta _{12}}{r^2}} }}{{{{\sin }^2}{\theta _{12}}}} = \sqrt {\frac{2}{{\sin\; {\theta _{12}}}}} r = {\text{DO}}{{\text{P}}_1}r $$ (25)

    式(25)显示在2个观测点距干扰源等距时,${\text{DO}}{{\text{P}}_{\text{3}}}$同样正比于$ {\text{DO}}{{\text{P}}_1} $,是DOP1$r$倍,与DOP2相比,是DOP2$\sqrt 2 $倍.

    $ {\text{DO}}{{\text{P}}_{\text{1}}} $的计算不受距离影响,当$N$=2时,$ {\text{DO}}{{\text{P}}_{\text{1}}} $形式如下:

    $$ {\text{DO}}{{\text{P}}_{\text{1}}} = \sqrt {\frac{2}{{{{\sin }^2}{\theta _{12}}}}} $$ (26)

    显然,当sin($ {\theta _{12}} = \pm 90\text{°} $)=±1时,$ {\text{DO}}{{\text{P}}_{\text{1}}} $取得最小值$\sqrt 2 $,此时夹角应满足$ {\theta _{12}} = $±90°.

    根据式(20)、(24)与(25)可知,在2点到观测点距离相同的情况下,3种DOP成等比例关系,而当观测点与干扰源位置形成的夹角为±90°时,3种DOP取到最小值.

    本部分通过仿真不同的观测点分布场景下DOP值和定位精度的变化情况,分析3种DOP是否能反映定位精度.

    首先进行2点定位的仿真,仿真场景条件如图2所示,在图2所示的$xOy$场景中,设定干扰源位置位于坐标轴的中心,观测点$ {O}_{1}、{O}_{2} $到干扰源的距离$ {d_1} $${d_2}$相等,均为500 m.

    图  2  测向交叉定位图(θ12=60°)

    ${O_1}$始终位于x轴正半轴上,${O_2}$则以干扰源为圆心转动,此时$ {\theta _1} $=0,${\theta _{12}}$=${\theta _2}$. 在仿真实验中改变${\theta _{12}}$,比较不同角度下各DOP和最终定位精度的变化情况.

    在实际定位过程中,定位误差的产生是因为测量值存在噪声,观测噪声通过定位方程传递给最后的定位结果.

    为了分析DOP值与误差大小的关系,需要在观测数据存在噪声的前提下进行实验. 给模拟实验的观测值增加位置噪声和角度噪声,此处以及后续的所有仿真中的噪声类型均为高斯白噪声,其中位置噪声的正态分布形式为$ \sigma _{{x}} = \sigma _{{y}} \sim {\mathrm{N}}\left( {0,10} \right) $,单位为m;角度噪声的正态分布形式为$ \sigma _{\theta}^{} \sim {\mathrm{N}}\left( {0,5} \right) $,单位为°. 进行多次定位,取均方根误差(root mean squared error,RMSE)作为定位精度的结果,最后得到${\theta _{12}}$与RMSE和各DOP的图像如图3所示.

    图  3  双观测点RMSE与DOP结果

    图3可知,在夹角${\theta _{ij}}$介于40°~140°之间时,RMSE与DOP的整体变化趋势接近,但数值上差距较大:$ {\text{DO}}{{\text{P}}_1} $过小,即使与测向误差综合也远远小于相应的定位误差;而$ {\text{DO}}{{\text{P}}_2} $$ {\text{DO}}{{\text{P}}_{\text{3}}} $的数值大小均以千计,显然大于定位误差的数值.

    但是在该角度区间内,3种DOP的变化趋势与干扰源定位误差的变化趋势基本一致. 夹角太小或太大时,前后观测点与干扰源近似共线,DOP与误差均急剧增加. 在夹角为80°~100°之间时,定位误差较小,最小值接近65 m,各DOP在此刻也相对较小.

    上述的仿真分析说明DOP可以反映角度对定位精度的影响,接下来需要研究其是否能反映距离的影响. 单从DOP值的表达式来看,$ {\text{DO}}{{\text{P}}_1} $无法反映距离对定位精度的影响,而$ {\text{DO}}{{\text{P}}_2} $$ {\text{DO}}{{\text{P}}_{\text{3}}} $均能反映.

    为了比较距离对DOP和定位精度的影响,改变$ {d_1} $${d_2}$为300 m和100 m,重新进行仿真,分别观察3种距离角度变化时定位精度和各DOP值的变化情况. 实验结果如图4所示.

    图  4  不同距离双观测点的RMSE和DOP结果

    图4可知,当观测点到干扰源的距离变化时,其定位结果的RMSE随着距离变小而变小,距离为500 m时最小定位RMSE约为60 m;距离为300 m时最小定位RMSE约为40 m;距离为100 m时,最小定位RMSE约为19 m. 此时3种距离下$ {\text{DO}}{{\text{P}}_1} $均相同,其无法反映距离变化带来的定位误差变化. $ {\text{DO}}{{\text{P}}_2} $$ {\text{DO}}{{\text{P}}_{\text{3}}} $随着距离变小而变小,可反映定位误差随距离的变化趋势.

    结合上述的仿真实验发现,在2点定位时,3种DOP值都可以反映角度对定位误差的影响,2点夹角在90°附近区域时定位精度较好,与1.4中的推导结果一致,而当夹角接近0°或180°时,即2个观测点与干扰源近似共线时,定位误差增加,而观测点距干扰越远,误差增加情况越严重.

    上一部分的模拟对2点定位时分布的2个基本因素的影响进行了分析,本部分绘制DOP值和定位精度分布图,对已有2点时继续增加新的观测点产生的影响进行比较.

    分布图的绘制过程为:选定随机2个观测点为固定点,把以干扰源为中心、附近范围2 000 m×2 000 m的区域划分为101×101的点阵. 将点阵中各点分别与固定点组合,进行多次交叉定位解算.

    为了避免特殊性,绘制分布图时选择2个到干扰源距离不同的观测点作为2个固定点,具体分布形式图5所示.

    图  5  固定点位置情况

    将3种DOP值和定位精度等数值以热力分布图的形式呈现出来,可较为直观的体现3种DOP值是否能够衡量定位结果的误差,继而确认最适合作为优化定位误差的指标.

    此时分布图结果如图6所示.

    图  6  定位误差与DOP值热力分布图

    图6可得出以下2个方面的结论:

    1)从DOP分布形状看,$ {\text{DO}}{{\text{P}}_{\text{3}}} $与RMSE最近似. 距离干扰源500 m以内处,$ {\text{DO}}{{\text{P}}_{\text{3}}} $与RMSE分布均为沙漏形状,在角度约为100°与280°处最小,在10°与190°处最大. 图6 (c)、(d)中红色虚线所围成的区域即对应$ {\text{DO}}{{\text{P}}_{\text{3}}} $和RMSE最小的区域.

    $ {\text{DO}}{{\text{P}}_1} $$ {\text{DO}}{{\text{P}}_2} $的分布均与RMSE分布相差较大. $ {\text{DO}}{{\text{P}}_{\text{2}}} $虽然能反映距离的影响,但是不同角度间的区别较小,同时$ {\text{DO}}{{\text{P}}_{\text{2}}} $最小的区域与RMSE最小的区域无法对应;$ {\text{DO}}{{\text{P}}_{\text{1}}} $则无法反映距离的影响. 以上说明当观测点到干扰源的距离不同时,$ {\text{DO}}{{\text{P}}_1} $$ {\text{DO}}{{\text{P}}_2} $均会受到较大的影响.

    2)从DOP的具体数值看,$ {\text{DO}}{{\text{P}}_2} $$ {\text{DO}}{{\text{P}}_{\text{3}}} $较大,数量级可达${10^3}$. 虽然$ {\text{DO}}{{\text{P}}_{\text{3}}} $大小分布和RMSE结果对应,但要通过其数值和测量误差的方差计算得到定位的RMSE比较困难. 相比之下,$ {\text{DO}}{{\text{P}}_1} $的大小虽然较为合理,但是由于其分布本身就与RMSE有较大差距,因此也无法进行合适的计算.

    根据分布图的变化,总结几何分布对定位精度的影响如下:观测噪声不变时,观测点与干扰源间的角度、距离均会对最终定位精度产生影响. 角度存在一个最优角度轴,在该轴线附近一定范围内的角度均会获得较好的定位精度;在角度确定的前提下,观测点距离干扰源越近,定位结果精度越高.

    综上所述,$ {\text{DO}}{{\text{P}}_1} $可以反映角度对定位误差的影响,但是无法反映距离的影响;$ {\text{DO}}{{\text{P}}_2} $可以反映角度和距离对定位误差的影响,但是某些情况下无法准确反映;$ {\text{DO}}{{\text{P}}_{\text{3}}} $的分布图与RMSE结果最为符合,对实际定位误差的描述最为准确.

    已有研究显示,若所有观测点与干扰源保持相同距离,并且分布为以干扰源为圆心的内接多边形时定位效果最好[15]. 然而,由于干扰源位置未知,并且难以提前确定需要多少个观测点,平均分布观测点的规划方法较难实现.

    此外,在实际运行过程中,观测点过于接近干扰源或将因远近效应导致测向模块失效. 因此,在实际规划时需考虑以下限制条件:

    1)无干扰源位置信息时,初始2个观测点位置无法指定,近似为随机位置.

    2)监测行动实际由1台移动站完成,后续新增观测点不能改变已有观测点位置,观测点的配置是按顺序增加的过程.

    3)规划监测载体路径时,设置观测点到干扰源的最小距离,观测点与干扰源间的距离不得小于该限制距离.

    在使用DOP值描述定位精度的过程中,$ {\text{DO}}{{\text{P}}_1} $仅能对角度的影响做出反映,$ {\text{DO}}{{\text{P}}_2} $$ {\text{DO}}{{\text{P}}_3} $可以同时反映角度和距离的影响,但$ {\text{DO}}{{\text{P}}_2} $的计算过程中需要重新计算距离,且描述精度不及$ {\text{DO}}{{\text{P}}_3} $. 理论上说,使用$ {\text{DO}}{{\text{P}}_3} $指导观测点规划的效果最好.

    下面将对增加观测点时$ {\text{DO}}{{\text{P}}_{\text{3}}} $的变化特性进行分析. 根据式(16),$ {\text{DO}}{{\text{P}}_{\text{3}}} $的分子部分记作${{{D}}_{{\mathrm{num}}}}$,分母记作${D}_{{\mathrm{den}}}$,即

    $$ \begin{array}{l}D_{\text{num}}=\sqrt{\begin{array}{l}g_1(f_{11}\mathrm{sin}\ \theta_1-f_{12}\mathrm{cos\ }\theta_1)^2+\dots+ \\ g_1(f_{21}\mathrm{sin}\ \theta_1-f_{22}\mathrm{cos}\ \theta_1)^2+\dots\end{array}} \\ D_{\text{den}}=\displaystyle\sum\limits_{i=1}^{N-1}\displaystyle\sum\limits_{j=i+1}^N\mathrm{sin}^2\theta_{ij} \\ \text{DOP}_{\text{3}}=\dfrac{D_{\text{num}}}{D_{\text{den}}}\end{array} $$ (27)

    在本文提出的定位环境中,每次进行规划时,先前已有观测点的位置和角度都被视为已有的常值,因此对$ {\text{DO}}{{\text{P}}_{\text{3}}} $进行求导时,可仅把最新的规划点${\theta _x}$作为未知数进行求导.

    $$ {{D}}_{3}{{'}}=({{D}}_{\text{num}}{{'}}\times {{D}}_{\text{den}}-{{D}}_{\text{num}}\times {{D}}_{\text{den}}{{'}})/({{D}}_{\text{den}})^{2} $$ (28)

    对式(22),由于分母始终为正数,此时仅考虑分子部分.

    $$ {{{D}}_{{d}}} = {{{D}}_{{\text{num}}}}{{'}} \times {{{D}}_{{\text{den}}}} - {{{D}}_{{\text{num}}}} \times {{{D}}_{{\text{den}}}}{'}$$ (29)

    假设N为1个较大值时,存在观测点$ {m_1} $,其角度${\theta _{{m_1}}}$使${\text{DO}}{{\text{P}}_{\text{3}}}$取到最小值,则有

    $$ \begin{gathered} {{{D}}_{{d}}}\left( {{\theta _{{m_1}}}} \right) = 0 \\ \mathop {\lim }\limits_{x \to \theta _{_{{m_1}}}^ - } {{{D}}_{{d}}}\left( x \right) = 0 \\ \mathop {\lim }\limits_{x \to \theta _{_{{m_1}}}^ + } {{{D}}_{{d}}}\left( x \right) = 0 \\ \end{gathered} $$ (30)

    上式中,第1个极限式数值小于0,第2个极限式数值大于0,两式符号相反.

    此时再增加1个新观测点$ {m_2} $,其对应角度为${\theta _{{m_2}}}$. 由于N为较大值,因此${{{D}}_{{d}}}$中的求和项目近似相等,仅影响与${{{D}}_{{\text{num}}}}'$中与$2{\theta _x}$有关的项. 当${\theta _{{m_2}}} = {\theta _{{m_1}}} + 90\text{°} $ 时,有

    $$ \begin{split} \mathrm{cos\ }2\theta_{m_2} & =-\mathrm{cos}\ 2\theta_{m_1} \\ \mathrm{sin}\ 2\theta_{m_2} & =-\mathrm{sin}\ 2\theta_{m_1} \\ \displaystyle\sum\limits_{i=1}^{N-1}\mathrm{sin}\ 2\theta_{im_2} & =\displaystyle\sum\limits_{i=1}^{N-2}\left(\mathrm{sin}\ 2\theta_{im_2}+\mathrm{sin}\ 2\theta_{m_1m_2}\right) \\ &=-\displaystyle\sum\limits_{i=1}^{N-2}\mathrm{sin}\ 2\theta_{im_1} \end{split} $$ (31)

    上式中${\theta _{i{m_2}}}$${\theta _{i{m_1}}}$${\theta _{{m_1}{m_2}}}$与1.1中规定一致,分别表示下标对应角度的差,通过观察可得

    $$ {{{\dot D}}_{{d}}}\left( {{\theta _{{m_2}}}} \right)\sim - {{{D}}_{{d}}}\left( {{\theta _{{m_1}}}} \right) = 0 $$ (32)

    对于$2{\theta _x}$,其角度增加180°之后变化趋势相反,此时${\theta _{{m_2}}}$为DOP3的1个极小值点,在${\theta _{{m_2}}}$可近似取得DOP3的最小值,此时有

    $$ \begin{gathered} \mathop {\lim }\limits_{x \to \theta _{_{{m_2}}}^ - } {{{D}}_{{d}}}\left( x \right) = 0 \\ \mathop {\lim }\limits_{x \to \theta _{_{{m_2}}}^ + } {{{D}}_{{d}}}\left( x \right) = 0 \\ \end{gathered} $$ (33)

    同样的,第1个极限式数值大于0,第2个极限式数值小于0,两式符号相反.

    此时,$ {m_2} $对应的$ {\text{DO}}{{\text{P}}_{\text{3}}} $最小,可以认为是新增的最优观测点,并且$ {m_1} $$ {m_2} $之间的角度差为90°. 在观测点$ {m_2} $的基础上继续增加观测点时,后续最优观测点也会出现类似情况,其角度仍与上一个最优观测点相差90°. 按顺序配置观测点时,当观测点的总数增加到一定大小,新配置的观测点角度与最后一个观测点角度差会始终为90°.

    结合推导结果和限制条件,对无人机监测系统在实际观测过程中的观测点提出3种规划方法:

    1)定步长法:规定观测点数量单次规划移动步长,以最新观测点为圆心,计算向所有方向前进1个步长后得到的新观测点的$ {\text{DO}}{{\text{P}}_{\text{3}}} $的大小,选择$ {\text{DO}}{{\text{P}}_{\text{3}}} $最小的点作为新的观测点. 重复上述步骤直至观测点数量达到限制.

    2)不定步长法:根据初始观测点获得的粗定位结果,在以粗定位结果为中心,定位盲区距离为半径的圆上选择$ {\text{DO}}{{\text{P}}_{\text{3}}} $最小的角度轴上的点,该点作为新的观测点,重复上述步骤直至观测点数量达到限制.

    3)均匀分布法:根据初始观测点获得粗定位结果,直接以该粗定位结果为中心进行目标数量的均匀分布配置进行定位.

    图7展示了3种方法在没有噪声时的规划结果. 第1种方法以螺旋形式向着干扰源的位置逐渐靠近;第2种方法则是在快速接近干扰源之后形成以干扰源为中心的正方形分布. 第3种方法则表现为获得粗定位结果后,舍弃初始点,直接进行均匀分布配置.

    图  7  无噪声规划效果

    规划方法2得到的结果与前文推导得到的结论符合,即当观测点数量足够多时,新增观测点与上一个点的夹角呈90°. 结合1.4与3.1中的推导,在顺序配置观测点的前提下,可将$ {\text{DO}}{{\text{P}}_{\text{3}}} $简化为$ {\text{DO}}{{\text{P}}_{\text{1}}} $,或是直接简化DOP值的计算,采用基于最新点增加90°的方式进行后续配置.

    图7的结果是在没有噪声的环境下产生的,但实际测量数据一定存在噪声. 当存在噪声时,这3种方法在初始2点的角度较好时,可获得定位精度理想的观测点配置结果. 但是由于每次配置都需遍历各角度的DOP值,计算量较大.

    图8所示,在遇到初始2点分布较差(即${\theta _{{{ij}}}} < 20\text{°} $)时,配置效果较差,导致最终定结果准确性受到影响.

    图  8  初始点分布较差时的规划效果

    基于上述规划方法计算量大以及对分布不好的初始观测点效果较差的2项缺陷,提出优化方法思路如下:

    1)根据初始观测点得到粗定位结果.

    2)不计算DOP值,直接以到粗定位结果距离减半、角度增加90°的原则增加新观测点.

    3)以新观测点重新计算定位结果,与原先的定位结果比较. 若定位结果的距离差距较大,则认为原先的定位结果与观测点不可靠,剔除该定位结果和未参与构建新点的观测点.

    4)以新得到的定位结果和观测点为基础继续增加观测点,重复2、3步骤,直到观测点数量达到要求.

    算法流程图如图9所示.

    相较于原来的3种规划方法,该算法在增加新点时简化了计算DOP值的过程,使得规划的计算量降低,并且通过比较前后的定位结果差距,将对定位精度贡献较差的观测点舍去,使得定位精度和规划速度都得到提升.

    图  9  算法流程图

    为了比较观测点配置后定位精度的效果,现对上述的方法进行误差仿真实验. 实验中,有干扰源设置为坐标系中心,初始2点的分布包括以下4种情况:

    情况1:(500,0),(536,449)

    情况2:(300,0),(230,193)

    情况3:(1500,0),(643,766)

    情况4:(300,0),(295,52)

    这4种情况的研究主要是为了比较4种方法在不同角度和距离时的精度情况. 各情况中初始点相对于干扰源的分布情况如图10所示.

    图  10  仿真实验点的分布图

    设置定步长法步长为100 m,定位盲区为200 m. 所有观测数据增加相同的角度噪声与位置噪声,进行500次仿真实验. 作为比较,将真实干扰源位置为中心的均匀分布的定位结果作为理想精度.

    实验结果如表1图11所示.

    表1可以看出,对于定步长方法,当其初始观测点距离较远时,不论角度如何,得到的RMSE均较差;距离较近时,即使初始角度较差,仍然能得到良好的RMSE. 不定步长方法在角度良好时效果良好,但距离过远时依然会受到影响;角度较差时,不管距离如何,不定步长法的RMSE都急剧增加.

    表  1  不同情况不同规划方法效果比较
    方法 初始点情况 最终RMSE/m 单次运行时间/s
    定步长 1 39.64 0.014307
    2 13.34
    3 91.71
    4 15.65
    不定步长 1 16.19 0.017772
    2 11.29
    3 26.83
    4 326.52
    均匀 1 13.17 0.000499
    2 10.09
    3 31.59
    4 5949.46
    优化 1 16.49 0.000351
    2 12.63
    3 16.29
    4 12.72
    理想 - 9.64 -
    下载: 导出CSV 
    | 显示表格
    图  11  各情况规划方法的定位精度

    均匀方法的问题相似,在初始角度较差时定位精度下降的更明显,RMSE接近6 000 m,此时可以认为已经无法完成定位工作. 相比之下,优化方法受到距离和角度的影响均较小,能够适合各种初始点分布情况.

    除了定位精度适应性方面的差距,从单次运行时间看来,定步长和不定步长的方法由于每次增加新的观测点都需要重新计算DOP,因此运行速度较慢. 优化方法省略了计算过程,直接通过增加90°来确定新的观测点,明显提高了运行速度.

    图11可知,角度合适的前提下,初始点距离干扰源较远时,定步长方法的RMSE与其它的方法差距较大,均匀规划的RMSE结果出现了一定的波动;距离干扰源较近时,所有方法均能达到较好的精度,但由于前2种方法不会舍弃初始观测点,在新增观测点数量较少时效果较差.

    角度不合适的前提下,不定步长规划和均匀规划的RMSE结果会显著增加,尤其是均匀规划;定步长规划在只有3~4个点时的RMSE显著增加,观测点数量增加之后逐渐恢复到正常水平. 只有优化规划始终在该情况下保持正常效果,因此认为优化规划在距离较远、角度较差的初始情况中也有良好的规划效果.

    本文对2点以及多点交叉定位时观测点与干扰源的相对分布情况对干扰源定位精度产生的影响展开仿真分析,给出3种DOP值来描述几何构型的好坏. 在对2点分布情况的分析中发现,2点角度差为80°~100°时,定位误差最小;2点距离干扰源的距离越近,定位误差越小. 在多点定位的分布图绘制中,DOP值和RMSE分布情况也符合上述规律,并且$ {\text{DO}}{{\text{P}}_{\text{3}}} $最能描述定位精度的变化.

    由于实际的干扰源监测中存在已有观测点无法更换、测向盲区等限制,因此基于实际的限制条件先提出3种观测点规划方法:定步长法、不定步长法、均匀分布法,进行测试后发现存在计算量大以及初始点角度差时定位效果不好的问题,因此在3种方法的基础上提出优化方法,能够较好地解决上述问题. 经过仿真比较发现,本文提出的优化规划可以在无法改变原有观测点的基础上获得较好的结果,定位精度足以满足干扰源监测系统的精度要求.

  • 图  1   测向定位原理图

    图  2   测向交叉定位图(θ12=60°)

    图  3   双观测点RMSE与DOP结果

    图  4   不同距离双观测点的RMSE和DOP结果

    图  5   固定点位置情况

    图  6   定位误差与DOP值热力分布图

    图  7   无噪声规划效果

    图  8   初始点分布较差时的规划效果

    图  9   算法流程图

    图  10   仿真实验点的分布图

    图  11   各情况规划方法的定位精度

    表  1   不同情况不同规划方法效果比较

    方法 初始点情况 最终RMSE/m 单次运行时间/s
    定步长 1 39.64 0.014307
    2 13.34
    3 91.71
    4 15.65
    不定步长 1 16.19 0.017772
    2 11.29
    3 26.83
    4 326.52
    均匀 1 13.17 0.000499
    2 10.09
    3 31.59
    4 5949.46
    优化 1 16.49 0.000351
    2 12.63
    3 16.29
    4 12.72
    理想 - 9.64 -
    下载: 导出CSV
  • [1] 靳睿敏, 郭艺, 杨会贇, 等. 基于导航接收机的GNSS弱干扰检测识别技术[J]. 全球定位系统, 2022, 47(6): 91-95.
    [2] 陈立军, 潘正军, 陈孝如. 无人机GPS欺骗干扰攻击对策[J]. 全球定位系统, 2023, 48(4): 91-98.
    [3] 吴涛, 胡艳霞, 田甜, 等. GNSS干扰定位技术分析[J]. 全球定位系统, 2023, 48(5): 103-111.
    [4] 张治, 卢鸿谦, 班晓军, 等. 基于MEMS/GNSS的多机协同无源定位[J]. 光学精密工程, 2021, 29(12): 2844-2854.
    [5] 聂大惟, 朱海, 吴飞, 等. 基于RSSI概率分布与贝叶斯估计的加权定位方法[J]. 全球定位系统, 2022, 47(2): 52-59. DOI: 10.12265/j/gnss.2021080902
    [6] 严宏基, 舒红, 孙红星. 几何精度衰减因子对两种到达时间差算法定位精度的影响分析[J]. 测绘地理信息2020, 45(3): 61-65.
    [7]

    FANG X P, LI J B, ZHANG S X, et al. Optimal AOA sensor-source geometry with deployment region constraints[J]. IEEE communications letters, 2022, 26(4): 793-797. DOI: 10.1109/LCOMM.2022.3144152

    [8]

    HAMDOLLAHZADEH M, AMIRI R, BEHNIA F. Optimalsensor placement for multi-source AOA localisation withdistance-dependent noise model[J]. IET radar, sonar &navigation, 2019(6): 13. DOI: 10.1049/iet-rsn.2018.5426

    [9] 王国刚. 两站测向交叉定位相对误差几何稀释度研究[J]. 舰船电子工程, 2015, 35(12): 58-61.
    [10] 罗双喜. 多站交叉定位相对GDOP及其测向站分布问题研究[J]. 指挥控制与仿真, 2020, 42(2): 7-11. DOI: 10.3969/j.issn.1673-3819.2020.02.002
    [11]

    PANWAR K, FATIMA G, BABU P. Optimal sensorplacement for hybrid source localization using fused TOA-RSS-AOA measurements[J]. IEEE transactions on aerospaceand electronic systems, 2023, 59(2): 1643-1657. DOI: 10.1109/TAES.2022.3202879

    [12] 聂文梅, 宋晓霞. 基于自适应粒子群优化算法的无线传感器网络覆盖控制[J]. 沈阳工业大学学报, 2023, 45(4): 459-464.
    [13] 姚志成, 许佳诺, 杨剑, 等. 一种改进的干涉仪-MUSIC联合测向算法[J]. 兵器装备工程学报, 2023, 44(11): 240-247.
    [14] 赵辰乾, 王松波, 刘益辰. 地球表面测向交叉定位算法[J]. 舰船电子对抗, 2021, 44(5): 86-89. DOI: 10.16426/j.cnki.jcdzdk.2021.05.019
    [15] 王本才, 何友, 王国宏, 等. 多站无源定位最佳配置分析[J]. 中国科学: 信息科学, 2011, 41(10): 1251-1267.
图(11)  /  表(1)
计量
  • 文章访问数:  96
  • HTML全文浏览量:  12
  • PDF下载量:  29
  • 被引次数: 0
出版历程
  • 收稿日期:  2024-05-26
  • 网络出版日期:  2024-12-24
  • 刊出日期:  2024-12-14

目录

/

返回文章
返回
x 关闭 永久关闭