非参数联合估计信号及位置的被动定位方法与流程

文档序号:11861993阅读:465来源:国知局
非参数联合估计信号及位置的被动定位方法与流程

本发明属于雷达信号处理技术领域,尤其涉及一种非参数联合估计信号及位置的被动定位方法。



背景技术:

近年来,由于日益激烈的电子战,使得在复杂环境下,对微弱目标进行高精度定位受到越来越多的关注。同时,芯片技术的飞速发展,使处理器的数据处理能力得到了极大改善,大数据的处理变得越来越现实。

被动定位技术由其诸多优点如隐蔽性、低功耗等,受到了越来越多的关注。被动定位通常采用几个分布较开的雷达基站截获目标发射机的发射信号,相比主动雷达,由于雷达自身不需要发射源发射电磁波,因此被动定位也被称为无源定位。

通常的被动定位技术主要针对目标为发射机的情况,也可称之为发射机定位,接收雷达通过处理截获的发射信号对目标实现定位。较为常用的定位方法是时差定位TDOA、到达角定位AOA。它们的基本思路是,首先在各基站接收机中提取出有关目标发射机的到达时间或到达角度,再将这些参数传输给处理中心,通过求解方程式的方法估计出目标位置。这种方法虽然传输的数据量小,计算简单,但是其定位精度低,无法满足高精度定位的需要。

文献“Direct Position Determination of Narrowband Radio Frequency Transmitters,IEEE Signal Process.Lett.,vol.11,no.5,pp.513-516,May 2004”提出了一种联合处理各基站接收机观测数据的直接定位方法DPD,该方法没有传统定位方法的参数提取过程,尽可能地保留了目标信息,仿真显示有明显的优势。这种方法考虑了两种情况,一种是目标发射机发射的信号已知,如训练信号或同步信号,我们称这种情况的DPD定位算法为DPD-known算法,由于知道了发射信号的信号波形,因此该方法是理论上最优的定位算法。另一种情况是目标发射的信号完全未知,主要针对非合作的发射机,在电子战中最为常见。由于不知道目标信号的形式,DPD定位方法通过使代价函数特征值最大的方法实现了发射机定位,仿真显示其定位性能仍然优于传统的定位方法,称这种方法为DPD-unknown算法。但是,该方法忽略了发射信号的波形信息,其定位精度受限,不能适应微弱目标的更高精度定位。



技术实现要素:

本发明的发明目的是:为了解决微弱非合作目标发射机精度受限的问题,本发明提出了一种高精度的非参数联合估计信号及位置的被动定位方法。

本发明的技术方案是:一种非参数联合估计信号及位置的被动定位方法,包括以下步骤:

A、初始化雷达系统参数,其中系统参数包括接收机个数L、各接收机位置坐标(xl,yl,zl)、l=1,2,...,L、采样间隔Ts、各接收机通道上的噪声协方差矩阵Rl

B、从各接收机中读取量测并对量测进行采样,得到L个离散的量测向量rl

C、分别划分目标位置的网格搜索区间(xgrid,ygrid,zgrid),信号发射时间的网格搜索区间tgrid,以及发射信号的长度的网格搜索区间Tgrid

D、选取一个网格点(xgrid,ygrid,zgrid,tgrid,Tgrid)的坐标参数作为发射信号的位置坐标(xgrid,ygrid,zgrid),发射时间tgrid,发射信号长度Tgrid,计算发射信号s0落在量测信号中对应的起始时间tl,s和终止时间tl,e

E、根据步骤D中得到的起始时间和终止时间计算发射信号s0在每个量测向量rl中对应的区间[nl,s,nl,e];

F、从每个接收机中的量测向量rl中提取第nl,s个元素到第nl,e个元素之间的离散样本点构成发射信号的样本向量s0,l,并计算提取的L个样本向量的功率系数E;

G、对步骤F中提取的每个样本向量s0,l进行归一化处理,将归一化处理后的L个样本向量组成矩阵S0′,并计算矩阵S0′的协方差矩阵Σ0

H、计算步骤G中协方差矩阵Σ0的特征值并将L个特征值按大小进行排列,再计算最大特征值对应特征向量ω1及最大特征值的成分系数ε,得到估计的发射信号并计算代价函数值

I、判断数据平面上的所有网格点是否已经被遍历;若数据平面上的所有网格点没有被遍历,则返回步骤D;若数据平面上的所有网格点被遍历,则选取最大代价函数值对应的网格点作为估计值,完成被动定位。

进一步地,所述步骤D中发射信号s0落在量测信号中对应的起始时间tl,s和终止时间tl,e的计算公式具体为:

<mrow> <msub> <mi>t</mi> <mrow> <mi>l</mi> <mo>,</mo> <mi>s</mi> </mrow> </msub> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>t</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> </mrow>

<mrow> <msub> <mi>t</mi> <mrow> <mi>l</mi> <mo>,</mo> <mi>e</mi> </mrow> </msub> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>t</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>T</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> </mrow>

其中,c为光速。

进一步地,所述步骤E中发射信号s0在每个量测向量rl中对应的区间[nl,s,nl,e]的计算公式具体为:

其中,表示向下取整。

进一步地,所述步骤F中计算提取的L个样本向量的功率系数E的计算公式具体为:

<mrow> <mi>E</mi> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <mo>|</mo> <msub> <mi>s</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>l</mi> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>/</mo> <msub> <mi>T</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> </mrow>

其中,|·|为求向量的模。

进一步地,所述步骤G中对步骤F中提取的每个样本向量s0,l进行归一化处理具体为:对提取的每个样本向量s0,l求平均得到并将样本向量s0,l中的每个元素减去均值得到归一化后的样本向量s′0,l

进一步地,所述步骤H中最大特征值的成分系数ε的计算公式具体为:

<mrow> <mi>&epsiv;</mi> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>/</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>&lambda;</mi> <mi>l</mi> </msub> <mo>;</mo> </mrow>

其中,λ1、λ2...λL为协方差矩阵Σ0的L个特征值,λ1≥λ2≥...≥λL

进一步地,所述步骤H中估计的发射信号具体表示为:

<mrow> <msub> <mover> <mi>s</mi> <mo>^</mo> </mover> <mn>0</mn> </msub> <mo>=</mo> <msub> <mi>S</mi> <mn>0</mn> </msub> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> </mrow>

其中,S0=[s0,1,s0,2,...,s0,L]。

进一步地,所述步骤H中计算代价函数值包括以下分步骤:

S1、根据估计的发射信号构造匹配向量sl,表示为:

<mrow> <msub> <mi>s</mi> <mi>l</mi> </msub> <mo>=</mo> <mo>&lsqb;</mo> <msub> <mn>0</mn> <mrow> <mi>l</mi> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <msub> <mover> <mi>s</mi> <mo>^</mo> </mover> <mn>0</mn> </msub> <mo>;</mo> <msub> <mn>0</mn> <mrow> <mi>l</mi> <mo>,</mo> <mn>2</mn> </mrow> </msub> <mo>&rsqb;</mo> </mrow>

其中,0l,1为由0组成的(nl,s-1)×1的列向量,0l,2为由0组成的(N-nl,e)×1的列向量;

S2、计算网格点(xgrid,ygrid,zgrid,tgrid,Tgrid)上的代价函数值,计算公式具体为:

本发明的有益效果是:本发明采用网格搜索的方法从接受机量测中提取出发射信号的样本向量,再对样本向量进行非参数化的处理估计出发射信号,然后利用估计出的信号计算代价函数,从而估计出发射机的位置,有效实现了对微弱发射机的高精度定位,且可适用的发射信号广泛,可以直接应用到现有被动雷达的高精度定位中。

附图说明

图1是本发明的非参数联合估计信号及位置的被动定位方法流程示意图。

图2是本发明实施例中发射机发射信号1时的定位性能对比示意图。

图3是本发明实施例中发射机发射信号2时的定位性能对比示意图。

图4是本发明实施例中发射机发射信号3时的定位性能对比示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅用以解释本发明,并不用于限定本发明。

如图1所示,为本发明的非参数联合估计信号及位置的被动定位方法流程示意图。一种非参数联合估计信号及位置的被动定位方法,包括以下步骤:

A、初始化雷达系统参数,其中系统参数包括接收机个数L、各接收机位置坐标(xl,yl,zl)、l=1,2,...,L、采样间隔Ts、各接收机通道上的噪声协方差矩阵Rl

B、从各接收机中读取量测并对量测进行采样,得到L个离散的量测向量rl

C、分别划分目标位置的网格搜索区间(xgrid,ygrid,zgrid),信号发射时间的网格搜索区间tgrid,以及发射信号的长度的网格搜索区间Tgrid

D、选取一个网格点(xgrid,ygrid,zgrid,tgrid,Tgrid)的坐标参数作为发射信号的位置坐标(xgrid,ygrid,zgrid),发射时间tgrid,发射信号长度Tgrid,计算发射信号s0落在量测信号中对应的起始时间tl,s和终止时间tl,e

E、根据步骤D中得到的起始时间和终止时间计算发射信号s0在每个量测向量rl中对应的区间[nl,s,nl,e];

F、从每个接收机中的量测向量rl中提取第nl,s个元素到第nl,e个元素之间的离散样本点构成发射信号的样本向量s0,l,并计算提取的L个样本向量的功率系数E;

G、对步骤F中提取的每个样本向量s0,l进行归一化处理,将归一化处理后的L个样本向量组成矩阵S0′,并计算矩阵S0′的协方差矩阵Σ0

H、计算步骤G中协方差矩阵Σ0的特征值并将L个特征值按大小进行排列,再计算最大特征值对应特征向量ω1及最大特征值的成分系数ε,得到估计的发射信号并计算代价函数值

I、判断数据平面上的所有网格点是否已经被遍历;若数据平面上的所有网格点没有被遍历,则返回步骤D;若数据平面上的所有网格点被遍历,则选取最大代价函数值对应的网格点作为估计值,完成被动定位。

在步骤A中,本发明对接收机个数L、各接收机位置坐标(xl,yl,zl)、采样间隔Ts、各接收机通道上的噪声协方差矩阵Rl等雷达系统参数进行初始化,具体为设置目标发射机的位置,雷达基站接收机的个数L=4,雷达基站接收机位置,初始化采样周期Ts=10-7s,噪声协方差矩阵Rl

在步骤B中,从各接收机中读取量测并对量测进行采样,产生用于仿真的雷达接收机量测向量rl,这里的量测向量为N×1的列向量,N表示量测样本个数,rl=[rl[0],rl[1],...,rl[N-1]]T,l=1,2,...,L。其中N表示量测样本个数,T表示转置操作。

在步骤D中,本发明选取一个网格点(xgrid,ygrid,zgrid,tgrid,Tgrid),并将其坐标参数作为发射信号的位置坐标(xgrid,ygrid,zgrid),发射时间tgrid和发射信号长度Tgrid,从而计算发射信号s0落在量测信号中对应的起始时间tl,s和终止时间tl,e,计算公式具体为:

<mrow> <msub> <mi>t</mi> <mrow> <mi>l</mi> <mo>,</mo> <mi>s</mi> </mrow> </msub> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>t</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> </mrow>

<mrow> <msub> <mi>t</mi> <mrow> <mi>l</mi> <mo>,</mo> <mi>e</mi> </mrow> </msub> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>t</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>T</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> </mrow>

其中,c为光速。

在步骤E中,本发明根据步骤D中得到的起始时间和终止时间计算发射信号s0在每个量测向量rl中对应的区间[nl,s,nl,e],计算公式具体为:

其中,表示向下取整。

在步骤F中,本发明采用从每个接收机中的量测向量rl中提取第nl,s个元素到第nl,e个元素之间的离散样本点的方式,构成发射信号的样本向量s0,l,再计算提取的L个观测向量的功率系数E,计算公式具体为:

<mrow> <mi>E</mi> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <mo>|</mo> <msub> <mi>s</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>l</mi> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>/</mo> <msub> <mi>T</mi> <mrow> <mi>g</mi> <mi>r</mi> <mi>i</mi> <mi>d</mi> </mrow> </msub> </mrow>

其中,|·|为求向量的模。

在步骤G中,本发明对步骤F中提取的每个样本向量s0,l进行归一化处理,具体为:对提取的每个观测向量s0,l求平均得到并将观测向量s0,l中的每个元素减去均值得到归一化后的观测向量s′0,l;再将归一化处理后的L个观测向量组成矩阵S0′,表示为S0′=[s′0,1,s′0,2,...,s′0,L],计算矩阵S0′的协方差矩阵Σ0,计算公式具体为:

<mrow> <msub> <mi>&Sigma;</mi> <mn>0</mn> </msub> <mo>=</mo> <msubsup> <mi>S</mi> <mn>0</mn> <mrow> <mo>&prime;</mo> <mi>H</mi> </mrow> </msubsup> <msubsup> <mi>S</mi> <mn>0</mn> <mo>&prime;</mo> </msubsup> </mrow>

其中,H表示共轭转置操作。

在步骤H中,本发明计算步骤G中协方差矩阵Σ0的特征值,并将计算得到的L个特征值按照从大到小的顺序进行排列,表示为λ1≥λ2≥...≥λL,再计算最大特征值λ1在所有特征值中所占比例ε,计算公式具体为:

<mrow> <mi>&epsiv;</mi> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>/</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>&lambda;</mi> <mi>l</mi> </msub> <mo>,</mo> </mrow>

计算最大特征值λ1对应的特征向量ω1,估计得到发射信号表示为其中S0=[s0,1,s0,2,...,s0,L],并计算代价函数值包括以下分步骤:

S1、根据估计的发射信号构造匹配向量sl,表示为:

<mrow> <msub> <mi>s</mi> <mi>l</mi> </msub> <mo>=</mo> <mo>&lsqb;</mo> <msub> <mn>0</mn> <mrow> <mi>l</mi> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <msub> <mover> <mi>s</mi> <mo>^</mo> </mover> <mn>0</mn> </msub> <mo>;</mo> <msub> <mn>0</mn> <mrow> <mi>l</mi> <mo>,</mo> <mn>2</mn> </mrow> </msub> <mo>&rsqb;</mo> </mrow>

其中,0l,1为由0组成的(nl,s-1)×1的列向量,0l,2为由0组成的(N-nl,e)×1的列向量;匹配向量sl的大小量测向量rl一致,为N×1的列向量;

S2、根据步骤S1中的匹配向量sl,及功率系数E,成分系数ε、量测向量rl计算网格点(xgrid,ygrid,zgrid,tgrid,Tgrid)上的代价函数值,计算公式具体为:

在步骤I中,本发明需要判断数据平面上的所有网格点是否已经被遍历;若数据平面上的所有网格点没有被遍历完,则返回步骤D,重新选择一个没有被遍历的网格点;若数据平面上的所有网格点已经被被遍历,则选取最大代价函数值对应的网格点作为估计值其中发射机的位置被估计为发射机发射信号的时间被估计为发射信号的长度被估计为根据即可完成被动定位。

本发明针对一个发射未知信号的发射机,采用多基地雷达接收机截获信号并将数据传输给处理中心进行定位;本发明采用了非参数化的方法估计出发射信号,能够应用于任意脉冲信号,其适用范围广泛,有效解决了现有的定位算法在低信噪比下定位性能不高的问题,从而实现对微弱目标的高精度定位,联合地估计了发射机发射的信号和位置,使定位精度得到了极大的提高。

在给定仿真参数下,对每一个信噪比(SNR)进行500次蒙特卡洛仿真实验,以均方根误差(RMSE)为定位性能依据,对比DPD-unknown算法,发现其定位性能有明显的优势。同时,以DPD-known算法作为定位误差的下界对比本发明方法在估计信号上的性能损失。为了方便表示,本发明的方法在图中被标注为DPD-enhanced。

仿真中,选取了三种不同发射信号分别验证本发明的定位性能。信号1为单载频脉冲信号,表达式为exp(-j2πf0t),信号频率f0为3MHz,信号长度为10us;信号2为高斯脉冲信号,其表达式为其中信号长度为Tp=20us,f0为3MHz,这里将信号能量归一化为1。信号3为线性调频信号,表达式为exp[-j2π(f0+0.5*kt)t],初始频率f0为3MHz,调频斜率k为1011Hz/s,信号长度为20us。三种不同信号的定位性能对比如图2、3、4所示。

如图2、3、4所示,为本发明实施例中发射机发射信号1、2、3时的定位性能对比示意图。,在高信噪比情况下,三个定位算法的定位误差都较低且相同;而随着信噪比的降低,三种定位方法性能出现差异。DPD-known由于假设知道发射信号的波形,因此始终获得了最好的定位性能。然而,在发射信号波形不知道的情况下,本发明的方法明显优于原有的DPD-unknown方法。在发射信号为信号1和信号2的情况下,本发明的方法较原有的DPD-unknown算法定位性能提高了约4dB;在发射信号为信号3的情况下,本发明的方法较原有的DPD-unknown算法定位性能提高了约3dB。由此可以看出,对大多数发射信号,本发明的方法均能明显优于原有的DPD-unknown算法。

本领域的普通技术人员将会意识到,这里所述的实施例是为了帮助读者理解本发明的原理,应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合,这些变形和组合仍然在本发明的保护范围内。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1