收藏 分销(赏)

基于传播算子的卫星导航系统干扰源直接定位方法.pdf

上传人:自信****多点 文档编号:633518 上传时间:2024-01-19 格式:PDF 页数:7 大小:1.90MB
下载 相关 举报
基于传播算子的卫星导航系统干扰源直接定位方法.pdf_第1页
第1页 / 共7页
基于传播算子的卫星导航系统干扰源直接定位方法.pdf_第2页
第2页 / 共7页
亲,该文档总共7页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述

1、第 21 卷 第 8 期2023 年 8 月Vol.21,No.8Aug.,2023太赫兹科学与电子信息学报Journal of Terahertz Science and Electronic Information Technology基于传播算子的卫星导航系统干扰源直接定位方法唐元春1,陈端云2,夏炳森1(1.国网福建省电力有限公司经济技术研究院,福建 福州 350000;2.国网福建省电力有限公司,福建 福州 350001)摘要:针对现有单星定位技术精确度不足问题,提出一种基于传播算子的卫星导航系统干扰源直接定位方法,解决现有技术无法快速有效地定位中高轨卫星干扰源的问题。对卫星接收信号

2、建模,计算每个时刻的接收信号协方差矩阵,通过对协方差矩阵分块,计算得到传播算子,进而得到导向矢量正交投影子空间的估计;然后建立损失函数,通过搜索得到干扰源位置估计结果。理论分析和仿真结果证明,该算法在复杂度大大降低的同时,定位精确度与传统基于多重信号分类(MUSIC)的直接定位算法相当,远优于现有的测向定位技术,具有重要的工程应用价值。关键词:导航卫星;传播算子;干扰源;直接定位中图分类号:TN96 文献标志码:Adoi:10.11805/TKYDA2021118Direct positioning method of interference source of satellite navi

3、gation system Direct positioning method of interference source of satellite navigation system based on propagation operatorbased on propagation operatorTANG Yuanchun1,CHEN Duanyun2,XIA Bingsen1(1.State Grid Fujian Economic Research Institue,Fuzhou Fujian 350000,China;2.State Grid Fujian Electric Pow

4、er Co.,Ltd,Fuzhou Fujian 350001,China)AbstractAbstract:To address the problem of insufficient precision of existing single-satellite positioning techniques,a direct positioning method of interference source for satellite navigation system based on propagation operator is proposed.Firstly,the receive

5、d signal of the satellite is modeled and the covariance matrix at each moment is calculated.After blocking the covariance matrix,the propagation operator is obtained to estimate the orthogonal projection subspace of the steering vector.Then the loss function can be established and the position of th

6、e interference source can be estimated by searching the loss function.Theoretical analysis and simulation results prove that the localization accuracy of the proposed method is comparable to that of the traditional Multiple Signal Classification(MUSIC)-based direct positioning algorithm,while the ca

7、lculation complexity is greatly reduced.Meanwhile,the localization accuracy is much higher than that of the existing bearing-only localization techniques,which is valuable for practical engineering applications.KeywordsKeywords:navigation satellite;propagation operator;interference source;direct pos

8、itioning卫星导航定位系统具有覆盖范围大,通信距离远,不受地理环境影响等优势。我国自主研制的北斗导航系统发展迅速,除了导航定位和授时外,还兼具短报文通信功能,已在交通、测绘、电信、气象等方面发挥着至关重要的作用,是我国综合国力的象征之一。由于卫星导航系统工作在一个开放的环境中,其导航定位等功能极易受到各种人为或无意的干扰。随着云计算、大数据、物联网、移动互联网、人工智能等新兴技术带动的新一轮科技革命和产业变革,无线电网络日益增多,卫星所处的电磁频谱环境日益恶化,复杂干扰信号的存在导致卫星导航质量下降,甚至会给国家带来巨大的经济损失。因此,快速有效地通过受扰卫星对干扰源进行定位是当务之急。

9、根据定位所需的卫星数量分类,卫星对地面干扰源的定位可分为单星定位体制和多星定位体制1-5。传统的双星和三星定位技术虽然定位精确度高,但需要多颗卫星协同工作,邻星的选择是一个难题。单星定位只需一文章编号:2095-4980(2023)08-0985-07收稿日期:2021-03-23;修回日期:2023-06-27太赫兹科学与电子信息学报第 21 卷颗卫星就能完成定位,不仅节省了有限的轨道资源,同时避免了邻星选择的困难。本文研究的主要对象为单星定位体制。现有的单星定位体制大都是测向交叉定位,整个定位过程分成参数估计和干扰源定位 2 个阶段。通过参数估计阶段对干扰源方向的估计,结合卫星与干扰源的空

10、间几何关系,进行干扰源位置解算6-7。但这种传统两步定位方法并没有充分考虑定位过程的整体性,每次参数估计步骤都是一个独立的过程,干扰源定位精确度完全依赖于参数估计的准确性,损失了整个定位过程中干扰源位置不变这一有效信息。由于低轨卫星的轨道高度相对较低,卫星可以截获功率相对较强的干扰源信号,因而可以获得较为准确的干扰源方位估计,卫星对干扰源的定位精确度容易达到需求。但对于卫星导航系统的中高轨卫星,由于卫星距离干扰源较远,一方面使测向带来的位置误差巨大;另一方面接收信号功率衰减严重,卫星接收到的干扰源信号较弱8-9,导致方位信息测量精确度不够,中间参数估计错误还会直接导致定位失败,不利于实际应用,

11、实现的技术难度非常大10。传统多重信号分类(MUSIC)等直接定位方法不需要估计定位参数过程,直接对原始信号进行采样处理,利用信号中的位置信息构建目标函数,通过穷尽搜索等算法实现定位11。本文针对现有的中高轨单星测向定位体制存在的问题,在干扰源位置不变的假设下,提出一种基于传播算子的卫星导航系统直接定位方法,系统性地考虑整个定位过程,利用传播算子及正交投影子空间的性质构造定位损失函数,通过直接搜索损失函数得到干扰源位置估计。算法避免了中间参数的估计步骤,且由于直接从原始数据层抽取位置信息,避免了信息的二次损耗,相比于传统测向定位算法,定位精确度得到了大幅提升。本文所提方法避免了传统基于 MUS

12、IC 等直接定位算法的特征值分解步骤,大大降低了计算复杂度。文中也给出了算法的计算复杂度分析,仿真结果表明该算法具有较高的估计精确度,计算效率高,是一种适用于中高轨卫星的快速有效的干扰源定位方法,同时也可拓展到无线电监测、雷达、勘探等领域。1系统模型系统模型如图 1 所示,建立地固坐标系OeX4Y4Z4,以地心为原点,OeX4由地心向外指向格林尼治子午线,OeZ4垂直于地球赤道平面与地球自转轴重合指向北极,OeY4与OeX4和OeZ4构成右手直角坐标系。在地理坐标系中,采用经度、纬度和高程来描述空间的位置,表示为H。地球半径为R,根据地理坐标系与地固坐标系之间的转换关系,可以计算出地固坐标系下

13、的空间位置为:x=(R+H)coscosy=(R+H)sincosz=(R+H)sin(1)干扰源地理位置矢量为p k=kkHkT,根据卫星星历得到卫星地固位置矢量为ul=xlylzlT,卫星接收机上集成了一阵元数为M的 L 阵,该 L 阵由 2 个阵元数相同的均匀线阵组成,所有阵元均被视为无差别的点,且对于各个方向的接收增益相等。传统的测向定位算法通常需要先从卫星截获的信号中估计出干扰源信号的方位角和俯仰角,利用卫星与干扰源之间的相对运动,通过空间几何关系才能实现定位。由于利用了方向信息,传统测向定位对卫星的姿态测量和控制精确度要求较高。本文采用的定位原理不需要通过中间参数的估计,通过构造损

14、失函数就能直接从原始数据层实现干扰源位置估计,大大避免了中间参数估计带来的误差。2定位算法2.1 算法描述为了简化模型,假设卫星接收的干扰源信号为窄带信号,且一段时间内对卫星存在连续干扰,则卫星在第l个位置的快拍t的接收信号为:xl(t)=Alsl(t)+nl(t)(2)式中:t=12J,J为快拍数;sl(t)=bl1s1(t)bl2s2(t)blKsK(t)T为阵列天线的参考阵元接收到的信号矢量,blk OeX4Z4Y4 u1ulpkFig.1 System model图1 系统模型986第 8 期唐元春等:基于传播算子的卫星导航系统干扰源直接定位方法为卫星接收信号的传输衰减系数;sk(t)

15、为干扰源k的发射信号,k=12K,K为干扰源个数;nl(t)(l=12L)为相互独立的零均值加性高斯白噪声,方差为2,且信号与噪声之间互不相关;Al=al(p1)al(p2)al(pK)为阵列流型矩阵,al(pk)=ejkTl(pk)d1ejkTl(pk)d2ejkTl(pk)dMT,其中,pk为干扰源地固位置矢量,dm表示阵列的第m个阵元相对于参考阵元的位置矢量,kl(pk)为卫星在位置ul处的波数向量,表示为:kl(pk)=2pk-ulpk-ul(3)式中:为波长;表示向量的 2 范数。由于实际信号采样长度有限,接收信号协方差矩阵由采样协方差矩阵代替,可计算为:Rl=1Jt=1Txl(t)

16、xHl(t)=AlRssAHl+2IM(4)式中:Rss=Esl(t)sHl(t)为信源信号协方差矩阵;IM为M M维单位矩阵。考虑到阵列流型矩阵Al为M K维列满秩矩阵,将Al分块为:Al=A1;A2(5)式中:A1 CK K为满秩矩阵;A2 C(M-K)K。显然,矩阵A2可以通过矩阵A1的线性变换得到:PA1=A2(6)式中 P 为(M-K)K维传播算子。将式(5)代入式(4),可得Rl=A1RssAHlA2RssAHl+2IM(7)考虑无噪声环境下的协方差矩阵Rl=Rl-2IM,令Rl=Ga;Gb,其中Ga为矩阵Rl的前K行,Gb为矩阵Rl的后(M-K)行,利用式(6)的关系,则有PGa

17、=Gb(8)由于实际环境中不可避免地存在噪声,因此,根据最小二乘理论,传播算子矩阵的估计为:P=minGHaPH-GHb2(9)式中:Ga和Gb分别由矩阵Rl的前K行和后(M-K)行构成,由式(9)可解得传播算子矩阵的估计值为:P=(GaGHa)-1GaGHbH(10)构造矩阵Pc=IK;P(11)式中IK为K K维单位矩阵,显然有PcA1=Al(12)根据式(12)可得矩阵Pc和矩阵Al由相同的列空间组成,即spanPc=spanAl。因此,可得导向矢量的正交投影子空间的估计为:Ql=Pc(PHcPc)-1PHc(13)利用导向矢量正交投影子空间的性质,可以构造出损失函数:f-1(p)=l=

18、1LaHl(p)(IM-Ql)al(p)(14)将搜索区域内的位置坐标 p 分别代入式(14),对上述损失函数进行搜索,找到式(14)的最大值,该最大值对应的位置坐标即为干扰源的位置估计结果。987太赫兹科学与电子信息学报第 21 卷2.2 搜索区域的选取由于地表范围很大,如果对整个地球表面进行一次全搜索,则定位效率很低。因此考虑利用卫星的对地覆盖区域来缩小定位搜索区域,从而加快卫星对干扰源的定位速度,这也是成功定位干扰源的保障。对于位于地面的静止干扰源(即干扰源在地理坐标系下的高度为 0),其可能存在的地固空间是无限的,而全球经纬度范围是有限的,因而考虑将干扰源可能存在的区域根据经纬度范围网

19、格化,如图 2所示。根据图 2 所示的卫星覆盖区域示意图,卫星位置点为S,干扰源可能存在的位置点为E,显然卫星覆盖区域内的位置点到卫星与地心连线的投影大于卫星与地球的切点到卫星与地心连线的投影,即向量OeE到向量OeS的投影应大于dmin。从图 2 所示的几何关系中,可以计算出dmin=R2(R+H)(15)同时通过前面介绍的系统模型可知地心Oe为坐标系原点,显然可得OeE=piOeS=ul(16)式中pi为干扰源可能存在的位置矢量,由离散化的经纬度网格点p i转换的地固坐标系下的坐标,则卫星在时刻l的对地覆盖区域可表示为:areal=pi|pieiR2(R+H)(17)式中:piei表示向量

20、pi与向量ei的点积,ei=ulul为卫星与地心连线方向上的单位向量。取L个时刻卫星覆盖区域的交集,可得卫星覆盖区域的重叠区域为:Area=pi|piarea1area2areaL2.3 算法步骤综上所述,本文所提算法的流程可归纳如下:步骤 1:选取卫星在L个不同时刻接收的干扰源信号,分别计算每个时刻的接收信号协方差矩阵;步骤 2:对协方差矩阵分块,计算得到传播算子,进而得到导向矢量正交投影子空间的估计;步骤 3:根据卫星的轨道信息计算出被干扰时卫星覆盖区域的重叠区域;步骤 4:利用导向矢量正交投影子空间的性质建立损失函数,将L个时刻的导向矢量正交投影子空间融合,通过搜索重叠区域得到干扰源位置

21、估计结果。2.4 复杂度分析测向定位(Angle of Arrival,AOA)方法把参数估计和干扰源定位分成 2 个独立的阶段,算法的计算复杂度主要取决于所选的测向算法的复杂度,而干扰源定位精确度完全依赖于参数估计的精确度。对于中高轨卫星,对测向精确度要求很高,复杂度较低的测向算法无法满足定位需求,因此以 MUSIC 测向算法进行分析。本文只考虑复数乘法的次数,MUSIC 算法的计算复杂度主要包括:协方差矩阵的求解,复杂度为O(JM2);协方差矩阵的特征值分解,复杂度为O(M3);谱函数搜索,复杂度为O()1M(M-K),1为角度搜索的网格点数。基于 MUSIC 的直接定位算法(DPD-MU

22、SIC)计算复杂度主要包括:协方差矩阵的求解,复杂度为O(JM2);协方差矩阵的特征值分解,复杂度为O(M3);谱函数搜索,复杂度为O()2M(M-K),2为位置搜索的网格点数。本文所提的算法无需特征值分解步骤,计算复杂度主要包括:协方差矩阵的求解,复杂度为O(JM2);传播算子矩阵Hdminsatellite u1ROeSElongitude latitude interference source piFig.2 Sketch map of satellite coverage area图2 卫星覆盖区域示意图988第 8 期唐元春等:基于传播算子的卫星导航系统干扰源直接定位方法的求解,复

23、杂度为O(K3+KM2+K2M);谱函数搜索,复杂度为O()2()K3+(K2+K)(M+1)。表 1 为 3 种干扰源定位算法复杂度对比。图 3 比较了特定参数下几种算法的复杂度随搜索间隔的变化。搜索间隔,在 AOA 方法中为方位角和俯仰角搜索间隔,方位角搜索范围为 0180,俯仰角搜索范围为090;在 DPD-MUSIC 和所提算法中为经度和纬度搜索间隔,经度搜索范围为 0180,纬度搜索范围为 090,搜索网格数目1和2根据设置的搜索范围和搜索间隔可以计算出。其他仿真参数设置为:L 阵阵元数 M=9,辐射源个数 K=1,每个位置的采样快拍数J=100。由图 3 可以看出,所提算法具有最低

24、的计算复杂度,且这种优势随着搜索间隔的减小而更加明显。3算法仿真为评估本文所提方法的有效性,使用蒙特卡洛仿真实验比较本文算法与传统测向定位算法 AOA 和 DPD-MUSIC的定位性能。仿真采用 MEO11 卫星轨道数据,卫星 4.4 s采集一次辐射源信号数据,卫星第一次采集时间与轨道数据零点重合。仿真实验中,采用 Matlab 作为仿真平台,采用均方根误差(Root Mean Square Error,RMSE)衡量算法的定位性能。假设pk为干扰源的真实位置,p ki为第i次蒙特卡洛仿真实验的位置估计结果,则 RMSE 为:RMSEk=1Ni=1Npk-p ki2(18)图 4 和图 5 分

25、别为本文所提方法在信噪比为 20 dB 和 10 dB 时的干扰源位置估计散点图。仿真场景设置为:选取卫星导航系统在轨 MEO 卫星在飞行轨道上的 3 个位置进行干扰源定位,将卫星的经度、纬度和高度转换成地固坐标系中的位置,卫星在地固坐标系下的 3 个位置坐标如表 2 所示。卫星接收机上集成了一个阵元数M=9的 L阵,每个位置的采样快拍数J=100;地面干扰源为连续窄带信号,其形式可表示为式(2);地面干扰源在地理坐标系下的位置为(82.4725 km,35.7645 km,0 km);经纬度搜索间隔均为 0.05,经度搜索范围为 0180,纬度搜索范围为 090;地球半径为R=6 400 k

26、m。根据式(1)中给出的地理坐标和地固坐标的转换方法,干扰源对应的地固坐标可以表示为(680 km,5148 km,3740 km)。从图中可以看出,在信噪比为 20 dB 时,x 方向上的位置估计误差均在50 km 以内,y 方向上的位置估计误差在 20 km 以内;在信噪比为 10 dB 时,x 方向上的位置估计误差均在 100 km以内,y 方向上的位置估计误差在 60 km 以内。由此可得,本文所提方法可有效实现 MEO 卫星对干扰源的定位。Fig.3 Comparison of the complexity versus search step for azimuth/elevati

27、on/longitude/latitude图3 特定参数下的算法复杂度随方位角/俯仰角/经度/纬度搜索间隔变化表1 算法复杂度对比Table1 Comparison of computational complexityalgorithm typeAOADPD-MUSICproposed algorithmcomputational complexityO(M3+JM2+1M(M-K)O(M3+JM2+2M(M-K)O(J+K)M2+K2M+K3+2(K3+(K2+K)(M+1)Fig.5 Scatter plot of interference source location under R

28、SN=10 dB and snapshot=100图5 RSN=10 dB,快拍数为100条件下的干扰源定位散点图Fig.4 Scatter plot of interference source location under RSN=20 dB and snapshot=100图4 RSN=20 dB,快拍数为100条件下的干扰源定位散点图989太赫兹科学与电子信息学报第 21 卷L 阵 和 面 阵 示 意 图 分 别 如 图 6 和 图 7 所 示,图 8 给出了在卫星配置 L 阵或面阵的情况下,本文所提方法与基于 MUSIC 的直接定位方法、AOA 方法在不同信噪比下的定位性能比较。由图

29、可以看出,随着信噪比或快拍数的增加,干扰源位置估计误差越来越小。本文所提方法在降低复杂度的情况下,定位误差与高复杂度的 DPD-MUSIC 相当,且二者定位误差远小于 AOA 方法。图 9 为本文所提方法在快拍数为 100 和 200 的情况下的定位误差收敛情况。仿真参数设置为:地面干扰源在地理坐标系下的位置为(82.4725 km,35.7645 km,0 km),选取卫星导航系统在轨 MEO 卫星在一段连续时间 500 s 内的位置,信噪比为 10 dB,L 阵阵元数为 11。由图可以看出,本文所提方法定位误差可以稳定收敛到 20 km 以内,具有很高的定位精确度,完全满足卫星干扰源定位精

30、确度要求。4结论本文针对现有的单星测向定位体制无法解决中高轨卫星干扰源定位问题,给出了一种基于传播算子的卫星导航系统干扰源直接定位方法。该方法有效避免了中间参数估计步骤,直接从原始数据层抽取干扰源位置信息,定位误差小,无需传统直接定位技术中复杂的特征值分解步骤,计算复杂度低,仅需一颗卫星即可完成干扰源定位。理论分析和仿真结果证明,该算法在复杂度大大降低的同时,定位精确度与传统基于 MUSIC 的直接定位算法相当,远优于现有的测向定位技术,具有重要的工程应用价值。参考文献:1 郝才勇,刘恒.采用空间稀疏性的单星无源定位方法J.电信科学,2016,32(4):65-69.(HAO Caiyong,

31、LIU Heng.Method of single satellite passive geolocation using spatial sparsityJ.Telecommunications Science,2016,32(4):65-69.)YN.XZMkkFig.6 Schematic diagram of L-shaped array图6 L阵示意图 YN.XZMsubarray 2subarray 1subarray N.kkinformation sourceFig.7 Schematic diagram of uniform planar array图7 均匀面阵示意图 01

32、015205RSN/dB106105104RMSE/m0 0.10.30.21051.401.351.301.25AOA(area array)DPD-MUSIC(area array)proposed(area array)AOA(L-shape)DPD-MUSIC(L-shape)proposed(L-shape)Fig.8 Comparison of the RMSE versus SNR under snapshot=100图8 快拍数为100条件下几种算法的定位RMSE随信噪比变化snapshot-100snapshot-2004.034.134.234.434.33t/(103 s

33、)RMSE/km9080706050403020100Fig.9 Convergence of localization error图9 干扰源定位误差收敛情况表2 仿真中MEO卫星在地固坐标系下的位置参数Table2 Position parameters of MEO satellite in the ground-fixed coordinate systemsatellite orbital altitudesatellite position 1satellite position 2satellite position 320 195 km(13 294,20 276,10 846

34、)km(13 295,20 290,10 819)km(13 296,20 304,10 791)km990第 8 期唐元春等:基于传播算子的卫星导航系统干扰源直接定位方法 2 DEMPSTER A G,CETIN E.Interference localization for satellite navigation systemsJ.Proceedings of the IEEE,2016,104(6):1318-1326.3 姜勤波,刘壮华,郑健,等.基于最大后验概率的单星多目标无源定位算法J.系统工程与电子技术,2014,36(10):1906-1912.(JIANG Qinbo,LI

35、U Zhuanghua,ZHENG Jian,et al.Algorithm for multi-targets passive localization with single satellite based on posterior probability maximizationJ.Systems Engineering and Electronics,2014,36(10):1906-1912.)4 HAO Benjian,AN Di,WANG Linlin,et al.A new passive localization method of the interference sour

36、ce for satellite communicationsC/2017 the 9th International Conference on Wireless Communications and Signal Processing(WCSP).Nanjing,China:IEEE,2017:1-6.5 GUO X S,LU H R,ZHANG Y,et al.Study on joint passive localization using time difference of arrival and frequency difference of arrival to improve

37、 the accuracy of four-satellite localizationJ.Advances in Mechanical Engineering,2017,9(12):1687814017737451.6 LUO J A,SHI Y F,HAN S T,et al.Bearing-only multi-target localization for wireless array networks:a spatial sparse representation approachC/2018 21st International Conference on Information

38、Fusion(FUSION).Cambridge,UK:IEEE,2018:1726-1731.7 LUO J A,SHAO X H,PENG D L,et al.A novel subspace approach for bearing-only target localizationJ.IEEE Sensors Journal,2019,19(18):8174-8182.8 王晓君,甄云双.宽带干扰背景下的多阵元卫星导航信号模拟J.太赫兹科学与电子信息学报,2020,18(2):208-214.(WANG Xiaojun,ZHEN Yunshuang.Simulation of satel

39、lite navigation signal on array antenna under the context of broadband interferenceJ.Journal of Terahertz Science and Electronic Information Technology,2020,18(2):208-214.)9 闻长远,岳富占,仇跃华.高轨GPS信号可用性分析J.电子设计工程,2014,22(2):29-33.(WEN Changyuan,YUE Fuzhan,QIU Yuehua.Analysis of high altitude GPS signal av

40、ailabilityJ.Electronic Design Engineering,2014,22(2):29-33.)10 王枭轩,陆兆峰,左小清,等.北斗系统融合单点性能分析J.城市勘测,2019(6):75-80.(WANG Xiaoxuan,LU Zhaofeng,ZUO Xiaoqing,et al.BDS system fusion single point performance analysisJ.Urban Geotechnical Investigation&Surveying,2019(6):75-80.)11 吴癸周,郭福成,张敏.信号直接定位技术综述J.雷达学报,2020,9(6):998-1013.(WU Guizhou,GUO Fucheng,ZHANG Min.Direct position determination:an overviewJ.Journal of Radars,2020,9(6):998-1013.)作者简介:唐元春(1973-),男,学士,高级工程师,主要研究方向为通信工程.email:.陈端云(1977-),男,硕士,高级工程师,主要研究方向为通信与电子工程.夏炳森(1981-),男,硕士,高级工程师,主要研究方向为通信与信息系统.991

展开阅读全文
相似文档                                   自信AI助手自信AI助手
猜你喜欢                                   自信AI导航自信AI导航
搜索标签

当前位置:首页 > 学术论文 > 论文指导/设计

移动网页_全站_页脚广告1

关于我们      便捷服务       自信AI       AI导航        获赠5币

©2010-2024 宁波自信网络信息技术有限公司  版权所有

客服电话:4008-655-100  投诉/维权电话:4009-655-100

gongan.png浙公网安备33021202000488号   

icp.png浙ICP备2021020529号-1  |  浙B2-20240490  

关注我们 :gzh.png    weibo.png    LOFTER.png 

客服