高级搜索

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

一种Alpha稳定分布噪声下目标螺旋桨特征提取方法

王彬 侯越圣

王彬, 侯越圣. 一种Alpha稳定分布噪声下目标螺旋桨特征提取方法[J]. 电子与信息学报, 2020, 42(10): 2478-2484. doi: 10.11999/JEIT190916
引用本文: 王彬, 侯越圣. 一种Alpha稳定分布噪声下目标螺旋桨特征提取方法[J]. 电子与信息学报, 2020, 42(10): 2478-2484. doi: 10.11999/JEIT190916
Bin WANG, Yuesheng HOU. Extraction of Target Propeller Features in Alpha Distribution Noise[J]. Journal of Electronics & Information Technology, 2020, 42(10): 2478-2484. doi: 10.11999/JEIT190916
Citation: Bin WANG, Yuesheng HOU. Extraction of Target Propeller Features in Alpha Distribution Noise[J]. Journal of Electronics & Information Technology, 2020, 42(10): 2478-2484. doi: 10.11999/JEIT190916

一种Alpha稳定分布噪声下目标螺旋桨特征提取方法

doi: 10.11999/JEIT190916
基金项目: 国家自然科学基金 (61572518)
详细信息
    作者简介:

    王彬:女,1971年生,副教授,研究方向为信号分析与处理

    侯越圣:男,1995年生,硕士生,研究方向为信号分析与处理

    通讯作者:

    王彬 commutech@163.com

  • 中图分类号: TN911.7

Extraction of Target Propeller Features in Alpha Distribution Noise

Funds: The National Natural Science Foundation of China(61572518)
  • 摘要: 为了解决Alpha稳定分布噪声下目标螺旋桨特征提取问题,该文提出一种基于分数低阶循环谱的特征提取方法。首先,从理论上推导出脉冲噪声条件下舰船辐射噪声分数低阶循环谱,指出分数低阶循环谱中出现峰值与螺旋桨特征的关系。然后根据该关系,提出基于分数低阶循环谱的螺旋桨特征估计方法。最后,通过仿真实验验证该方法的性能,并通过实测数据进一步验证了算法有效性。
  • 螺旋桨是舰船在海上航行的动力来源,不同类型目标的螺旋桨差别较大,水面舰艇的螺旋桨通常有3~4片桨叶,潜艇螺旋桨叶片数较多通常有7~8叶。如果能对目标螺旋桨特征进行准确的识别,将为判断目标类型提供重要的参考依据[1]

    辐射噪声幅度上呈现的规律变化会使其听觉上具有鲜明的节奏感,陶笃纯[2,3]对节奏感的产生机理进行研究,将这种现象建模为多个具有随机幅度、相同形状和相同重复周期的随机高斯脉冲,并以该模型为基础证明对包络进行功率谱分析可估计目标轴频、叶频和螺旋桨叶片数等螺旋桨特征。Lourens等人[4]在假设空化噪声和背景噪声服从高斯分布且相互独立情况下,推导出了轴频最大似然估计器。童峰等人[5]针对DEMON过程中低信噪比和本舰干扰的问题,提出一种基于模糊判决的线谱识别方式,并利用多普勒效应区分DEMON谱中线谱的来源,抑制本舰干扰提取目标轴频。实际中,目标产生的空化噪声以高斯色噪声刻画更为贴近,Clark等人[6]提出了高斯色噪声下的特征估计方法。Wisdom等人[7]通过增加相关时间解决DEMON过程中调制频率变化的问题。3/2维谱具有抑制加性高斯噪声的优点,被用于代替传统DEMON过程中的功率谱估计步骤,可获得更好的谱分析结果用于螺旋桨特征提取[8-10]。不同频段获得的DEMON谱效果差别较大,杨日杰等人[11]将不同频段获得的DEMON谱进行融合处理,增强特征线谱,并提出一种序列匹配方法实现轴频自动提取。

    总地来说,现在有的螺旋桨特征提取方法大都假设环境噪声是高斯分布的。但在浅海环境下,由于人类活动、海洋生物和自然因素的影响,使得接收信号存在大量的脉冲干扰。这种噪声不再服从高斯分布,通常用Alpha稳定分布来刻画[12,13]。Alpha稳定分布噪声的2阶及其以上各阶统计量是不存在的,这使得适用于高斯噪声环境下的传统方法在Alpha稳定分布噪声下性能退化甚至失效。

    本文提出基于分数低阶循环谱的螺旋桨特征提取方法,解决Alpha稳定分布噪声下目标螺旋桨特征提取问题。针对舰船辐射噪声表现出的非平稳性,本文从分数低阶循环平稳角度出发,从理论上对分数低阶循环谱提取包络调制频率的可行性进行论证。然后,基于理论分析结果给出了Alpha稳定分布噪声下目标螺旋桨特征提取方法,在不增加过多计算量的情况下可实现目标螺旋桨特征提取。理论推导和实验结果均证明了所提方法在Alpha稳定分布噪声下的有效性。

    在浅海脉冲环境下,接收到的辐射噪声r(t)可表示为

    r(t)=s(t)+nenv(t)
    (1)

    其中s(t)表示舰船辐射噪声,nenv(t)为环境噪声,s(t)可表示为

    s(t)=(1+m(t))ncav(t)
    (2)

    ncav(t)为螺旋桨空化噪声,服从均值为0、方差为σcav2的高斯分布。m(t)为舰船辐射噪声的包络调制信号,可将其建模为具有相同重复周期、相同形状、幅度随机的高斯脉冲随机过程。螺旋桨叶片数为B,每B个脉冲为1组,则2K+1组高斯脉冲序列可表示为[2]

    m(t)=Kk=K[B1i=0ξ(k)nu(tiTBkT)]
    (3)

    其中T为螺旋桨旋转周期,单个脉冲表达式为

    u(t)=12πet22σ2
    (4)

    每个高斯脉冲的持续时间T0=1/Bf0, f0为轴频。ξ(k)n是每个脉冲的幅值,其概率密度函数可认为在(ξn/2,3ξn/2)内均匀分布,ξn(n=1,2,···,B)是每个桨叶产生高斯脉冲幅度均值。

    式(1)中的nenv(t)表示与s(t)相互独立的环境噪声,服从Alpha稳定分布。对于Alpha稳定分布噪声而言,α是噪声特征指数。当α<2时可用于刻画脉冲噪声,α越小表示噪声的脉冲性越强,当α=2时表示高斯白噪声。因此,用Alpha稳定分布噪声来表征浅海环境噪声更具普适性。由于Alpha稳定分布噪声方差的理论值是无穷大的,因此通常用混合信噪比(Mixed Signal-to-Noise Ratio, MSNR)

    MSNR=10lg(σ2s/γ)
    (5)

    来描述信号功率和噪声功率的关系,σ2s为信号平均功率,γ(γ0)为噪声分散系数。

    DEMON谱是提取螺旋桨特征的重要工具,其中线谱出现的频率与螺旋桨轴频、叶频和螺旋桨叶片数等特征有关。DEMON谱获取流程为:首先利用带通滤波器滤除低频线谱,然后送入解调器获得信号的包络,而后对解调获得的包络波形做谱估计即可获得DEMON谱,如图1所示[14]

    图 1  DEMON谱估计流程

    目前,常用的解调方法有平方解调和绝对值解调,但这两种方法在Alpha稳定分布噪声环境中性能下降。主要原因是服从Alpha稳定分布的随机变量,当0<α<2

    E[|X|p]=,pα
    (6)
    E[|X|p]<,0p<α
    (7)

    其中,E[]表示期望。特别地,分数低阶Alpha稳定分布随机变量不存在有限的2阶矩。显然,信号的2阶统计量不收敛,所以在Alpha稳定分布噪声环境下,基于平方/绝对值解调的DEMON谱分析性能退化甚至失效。

    对于特征指数为α的Alpha稳定分布噪声,只有阶数小于α的分数低阶统计量是有限的。分数低阶循环谱是其中的一类重要概念,在Alpha稳定分布噪声下处理非平稳信号中发挥了重要作用[15,16]

    目标螺旋桨特征提取是通过分析辐射噪声信号的包络调制频率实现,为了表示方便将辐射噪声模型简化为单频包络调制信号的形式

    x(t)=(1+mcos(2πf0t))n(t)
    (8)

    其中m为调制度,满足0<m<1, f0为调制频率,n(t)为高斯白噪声。其分数低阶循环谱表达式如式(9)所示(推导过程略)

    Sεxb(f)={(1+m2b(b1)4)δ(f)+b24m2δ(ff0) +b2(b1)264m2δ(f2f0),ε=0(1+m2b(b1)4)b2ejθδ(ff0/2)+b2(b1)16m2ejθδ(f5f0/2), ε=±f0b24m2ej2θδ(f)+(1+m2b(b1)4)b(b1)8mej2θδ(ff0),  ε=±2f0b2(b1)16m2ej3θδ(f3f0/2), ε=±3f0b2(b1)264m2ej4θδ(f),  ε=±4f00,其他
    (9)

    式中f为频率、ε为循环频率、b为分数低阶循环谱的阶数。从式(9)中可以看出ε=0的分数低阶循环谱中会出现与包络调制频率有关的谱线。对于式(3)所表示的高斯脉冲包络而言,其ε=0的分数低阶循环谱表达式为(推导过程略)

    Sεmb(f)=|gb(f)|2T{(Ub0Ub5)+2πT(Ub5+2Ub1cos2πfT5+2Ub2cos4πfT5)n=δ(2πf2πnT)}
    (10)

    其中

    U0b=(ξb0)2+(ξb1)2+(ξb2)2+(ξb3)2+(ξb4)2U1b=ξb0ξb1+ξb1ξb2+ξb2ξb3+ξb3ξb4+ξb4ξb0U2b=ξb0ξb2+ξb1ξb3+ξb2ξb4+ξb3ξb0+ξb4ξb1U5b=(ξb0)2+(ξb1)2+(ξb2)2+(ξb3)2+(ξb4)2}
    (11)

    当螺旋桨为非均匀调制时,即每片桨叶的脉冲幅值不完全相等的情况,ε=0的分数低阶循环谱的第1根线谱对应螺旋桨轴频,能量最高线谱对应螺旋桨叶频,两者频率倍数差为螺旋桨叶片数。如果为均匀调制,即每片桨叶的脉冲幅值相等时,只有叶频谱线,没有轴频谱线。

    通过上一小节对信号分数低阶循环谱表达式的推导,发现其频率轴和循环频率轴都会出现与包络调制频率有关的线谱信息,可用于分析辐射噪声的包络调制频率。但计算分数低阶循环谱的过程复杂,计算量太大。ε=0的分数低阶循环谱也出现了与调制频率有关的峰值。从计算简便角度出发,本文提出如下螺旋桨特征估计方法,该方法流程如图2所示,具体步骤如下:

    图 2  螺旋桨特征提取方法流程

    步骤 1 对信号进行带通滤波,去除低频线谱的影响,选择幅度调制较强的频段进行分析;

    步骤 2 计算滤波后信号ε=0的分数低阶循环谱;

    步骤 3 从获得的谱中提取目标螺旋桨特征,出现的第1根线谱对应螺旋桨轴频,能量最高的线谱对应螺旋桨叶频,两者的倍数差为螺旋桨叶片数。

    以五叶非均匀调制螺旋桨为例进行仿真,设螺旋桨轴频3.6 Hz,各桨片的脉冲幅值的均值分别为1.0, 0.5, 0.7, 0.4, 0.6,包络调制后的仿真信号时域波形如图3所示。图4图6分别是在MSNR=10 dB, α=2.0, α=1.5α=1.2情况下的分析结果,图(a),图(b),图(c)依次为平方解调、绝对值解调和本文方法。

    图 3  舰船辐射噪声波形
    图 4  α=2.0时的谱分析结果
    图 5  α=1.5时的谱分析结果
    图 6  α=1.2时的谱分析结果

    α=1.5时,平方解调DEMON已经失效,绝对值解调DEMON谱的线谱分量也几乎被淹没难以辨认。在上述条件下,本文方法性能稳定,可实现目标螺旋桨特征提取。从图6(c)中可以看出,第1根线谱频为3.611 Hz,能量最高的线谱频率为18.03 Hz为螺旋桨叶频,两者相差约5倍。估计结果与设定轴频3.6 Hz的非均匀调制五叶螺旋桨相一致。仿真条件不变,改变噪声特征参数从α=1.1开始以0.1为间隔到α=2.0,每个间隔做100次实验。通过搜寻谱中能量最大值出现的频率,对螺旋桨叶频进行估计。以正确率衡量不同方法的叶频估计性能。若目标叶频的真实值为Y,估计值为ˆY,当两者满足|(YˆY)/Y|101认为估计正确,定义正确率

    正确率(%)=正确估计次数总实验次数×100%
    (12)

    对不同方法性能进行衡量,不同方法在不同噪声特征指数下的正确率如图7所示。

    图 7  噪声特征指数对叶频估计性能影响

    图7中本文方法曲线是设定阶数b=0.2获得的结果,可以看出本文方法在1.1α2.0均可稳健估计螺旋桨的叶频。

    仿真条件不变,在α=1.5情况下测试不同混合信噪比下叶频估计性能。混合信噪比从–10~26 dB以2 dB为间隔,每组参数下进行100次测试,正确率结果如图8所示。脉冲噪声对特征提取的影响在于其功率较大,会将有用线谱淹没。噪声特征指数不变的情况下,随着混合信噪比的升高,噪声功率相对减小,因此基于传统解调的叶频提取方法可在高混合信噪比条件下实现叶频估计。在α=1.5情况下,绝对值检波在混合信噪比高于10 dB时性能达到稳定,平方检波则需要更高的混合信噪比条件才能取得较好的性能。本文方法,通过分数低阶处理抑制了脉冲噪声的功率,因此在混合信噪比高于–2 dB时已经具有较好性能。

    图 8  混合信噪比对叶频估计性能影响

    在人工湖中进行实验多次实验,采用icListen HF水听器采集天和防务TH-B050R型AUV的辐射噪声。AUV动力装置为一个三叶螺旋桨,螺旋桨的转速可通过岸段控制台在300/400/500(r/min)中选定。实验区域水深1.5 m,测试时水听器位于距离岸边0.5 m、水深0.5 m处。控制AUV做匀速直线运动从水听器附近通过,图9所示的辐射噪声信号波形是在AUV与水听器的横正距离1.5 m,螺旋桨转速为400 r/min(6.67 r/s)情况下采集获得。

    图 9  AUV辐射噪声时域波形

    螺旋桨旋转引起的辐射噪声信号幅度调制是全频带的,但由于环境噪声分布不均、不同频段调制指数不同,在对包络中包含的螺旋桨特性进行提取时,不同频段的结果会出现差异。采用本文方法对不同频段进行分析获得的ε=0的分数低阶循环谱如图10所示。

    图 10  AUV辐射噪声不同频段循环谱分析结果

    图10(b)中出现第1根谱线频率为6.256 Hz对应目标轴频,与设定转速400 r/min(6.667 r/s)相差0.411 Hz。从图10(b)可以看出有3根较为明显的线谱,第3根线谱频率为18.73 Hz对应目标叶频,约为轴频的3倍,可认为该目标是由三叶螺旋桨推进,估计结果与目标螺旋桨叶片数相同。对比图10不同频段的解调结果,可以看出对实际信号进行分析时,解调频带的选择会对最后结果有较大影响,在某区间不是很明显的线谱在其他频带可能较为明显。例如,在图10(a)可以看到基频和倍频线谱,在图10(b)可以看到较为明显的3根线谱,图10(c)16 kHz以上的频段获得的则看不出包络调制现象。但是如何选取到最优的频段并没有理论上的指导。

    通过岸段控制台改变AUV螺旋桨转速,获得不同工况下的AUV辐射噪声。采用本文方法对8~16 kHz频段进行分析获得的轴频估计结果如表1所示。多次测试结果获得的轴频估计值与螺旋桨转速设定值较为接近,与转速设定值的相对误差均在6%以内。

    表 1  AUV轴频估计结果
    距离(m)11222.5
    螺旋桨转速
    设定值(r/min)
    300400300500400
    轴频估计值(Hz)4.716.344.737.936.28
    估计误差(Hz)0.290.330.270.400.39
    下载: 导出CSV 
    | 显示表格

    本文提出一种Alpha稳定分布噪声下目标螺旋桨特征提取方法。首先,对分数低阶循环谱估计信号包络调制频率的可行性进行理论分析。然后探究了ε=0的分数低阶循环谱中出现的峰值与螺旋桨特征的关系。此基础上,给出了一种计算简便的螺旋桨特征估计方法。仿真结果表明:本文方法对噪声特征指数1α2的Alpha稳定分布噪声具有稳健性,在混合信噪比高于–2 dB能够有效估计目标轴频。采用本文方法对实测辐射噪声进行分析,螺旋桨特征提取结果与AUV螺旋桨特征一致,验证了本文方法有效性。

  • 图  1  DEMON谱估计流程

    图  2  螺旋桨特征提取方法流程

    图  3  舰船辐射噪声波形

    图  4  α=2.0时的谱分析结果

    图  5  α=1.5时的谱分析结果

    图  6  α=1.2时的谱分析结果

    图  7  噪声特征指数对叶频估计性能影响

    图  8  混合信噪比对叶频估计性能影响

    图  9  AUV辐射噪声时域波形

    图  10  AUV辐射噪声不同频段循环谱分析结果

    表  1  AUV轴频估计结果

    距离(m)11222.5
    螺旋桨转速
    设定值(r/min)
    300400300500400
    轴频估计值(Hz)4.716.344.737.936.28
    估计误差(Hz)0.290.330.270.400.39
    下载: 导出CSV
  • 白敬贤, 高天德, 夏润鹏. 基于DEMON谱信息提取算法的目标识别方法研究[J]. 声学技术, 2017, 36(1): 88–92. doi: 10.16300/j.cnki.1000-3630.2017.01.016

    BAI Jingxian, GAO Tiande, and XIA Runpeng. Target recognition based on the information extraction algorithm of DEMON spectrum[J]. Technical Acoustics, 2017, 36(1): 88–92. doi: 10.16300/j.cnki.1000-3630.2017.01.016
    陶笃纯. 舰船噪声节奏的研究(I)——数学模型及功率谱密度[J]. 声学学报, 1983, 8(2): 65–76. doi: 10.15949/j.cnki.0371-0025.1983.02.001

    TAO Dunchun. A study on ship radiated noise rhythms (I) — mathematical model and power spectrum density[J]. Acta Acustica, 1983, 8(2): 65–76. doi: 10.15949/j.cnki.0371-0025.1983.02.001
    陶笃纯. 舰船噪声节奏的研究(Ⅱ)——自相关函数及节奏信息的提取[J]. 声学学报, 1983, 8(5): 280–289. doi: 10.15949/j.cnki.0371-0025.1983.05.004

    TAO Duchun. A study on ship-radiated noise rhythms (Ⅱ) — atuocorrelation function and extraction of rhythm information[J]. Acta Acustica, 1983, 8(5): 280–289. doi: 10.15949/j.cnki.0371-0025.1983.05.004
    LOURENS J G and DU PRCEZ J A. Passive sonar ML estimator for ship propeller speed[J]. IEEE Journal of Oceanic Engineering, 1998, 23(4): 448–453. doi: 10.1109/48.725238
    童峰, 陆佶人, 方世良. 水声目标识别中一种轴频提取方法[J]. 声学学报, 2004, 29(5): 398–402. doi: 10.15949/j.cnki.0371-0025.2004.05.003

    TONG Feng, LU Jiren, and FANG Shiliang. A method to estimate the propeller shaft rate in underwater acoustic vessel classification[J]. Acta Acustica, 2004, 29(5): 398–402. doi: 10.15949/j.cnki.0371-0025.2004.05.003
    CLARK P, KIRSTEINS I, and ATLAS L. Multiband analysis for colored amplitude-modulated ship noise[C]. 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, Dallas, USA, 2010: 3970–3973.
    WISDOM S, ATLAS L, and PITTON J. Extending coherence time for analysis of modulated random processes[C]. 2014 IEEE International Conference on Acoustics, Speech and Signal Processing, Florence, Italy, 2014: 340–344.
    LI Sichun and YANG Desen. DEMON feature extraction of acoustic vector signal based on 3/2-D spectrum[C]. The 2007 2nd IEEE Conference on Industrial Electronics and Applications, Harbin, China, 2007: 2239–2243.
    温旋旋. 基于高阶统计量的舰船辐射噪声包络特征提取[D].[硕士论文], 哈尔滨工程大学, 2017.

    WEN Xuanxuan. Envelope signature extraction of ship-radiated noise based on HOS[D].[Master dissertation], Harbin Engineering University, 2017.
    许劲峰, 郑威. 基于EMD-1(1/2)维谱的舰船辐射噪声调制特征提取[J]. 舰船电子工程, 2018, 38(10): 197–203. doi: 10.3969/j.issn.1672-9730.2018.10.048

    XU Jinfeng and ZHENG Wei. Modulation feature extraction of ship radiated nose based on EMD 1(1/2)-dimensional spectrum[J]. Ship Electronic Engineering, 2018, 38(10): 197–203. doi: 10.3969/j.issn.1672-9730.2018.10.048
    杨日杰, 郑晓庆, 韩建辉, 等. 基于序列匹配的螺旋桨轴频自动提取方法[J]. 振动与冲击, 2018, 37(16): 57–61. doi: 10.13465/j.cnki.jvs.2018.16.009

    YANG Rijie, ZHENG Xiaoqing, HAN Jianhui, et al. An automatic extraction method of propeller shaft frequency based on sequence matching[J]. Journal of Vibration and Shock, 2018, 37(16): 57–61. doi: 10.13465/j.cnki.jvs.2018.16.009
    CHITRE M, POTTER J, and HENG O S. Underwater acoustic channel characterisation for medium-range shallow water communications[C]. OCEANS ‘04 MTS/IEEE TECHNO-OCEAN ’04, Kobe, Japan, 2004: 40–45.
    张安清, 邱天爽, 章新华. α稳定分布的水声信号处理新方法[J]. 电子与信息学报, 2005, 27(8): 1201–1204.

    ZHANG Anqing, QIU Tianshuang, and ZHANG Xinhua. A new underwater acoustic signals processing approach to α-stable distribution[J]. Journal of Electronics &Information Technology, 2005, 27(8): 1201–1204.
    程玉胜, 王易川, 史广智, 等. 基于现代信号处理技术的舰船噪声信号DEMON分析[J]. 声学技术, 2006, 25(1): 71–74. doi: 10.3969/j.issn.1000-3630.2006.01.016

    CHENG Yusheng, WANG Yichuan, SHI Guangzhi, et al. DEMON analysis of underwater target radiation noise based on modern signal processing[J]. Technical Acoustics, 2006, 25(1): 71–74. doi: 10.3969/j.issn.1000-3630.2006.01.016
    张俊林, 刘明骞. 脉冲噪声衰落信道下MPSK信号的符号周期估计[J]. 西安电子科技大学学报: 自然科学版, 2018, 45(5): 38–42, 49. doi: 10.3969/j.issn.1001-2400.2018.05.007

    ZHANG Junlin and LIU Mingqian. Estimation of the symbol period of MPSK signals over fading channels with impulsive noise[J]. Journal of Xidian University:Natural Science, 2018, 45(5): 38–42, 49. doi: 10.3969/j.issn.1001-2400.2018.05.007
    张娟娟. α稳定分布噪声下数字调制信号的分数低阶循环谱分析[D]. [硕士论文], 西安理工大学, 2017.

    ZHANG Juanjuan. Fractional Low-order cyclic spectrum analysis of digital modulation signals under alpha stable distribution noise[D]. [Master dissertation], Xi’an University of Technology, 2017.
  • 期刊类型引用(3)

    1. 江莉,向世召,焦予栋,尚文擎,王燕妮. 一种NLF-对数小波阈值的脉冲干扰降噪方法. 小型微型计算机系统. 2023(12): 2844-2850 . 百度学术
    2. 李新. 功率谱估计在舰船噪声特征提取中的应用仿真. 舰船科学技术. 2022(04): 43-46 . 百度学术
    3. 孟庆松,段建红,黄子豪. 基于LFM信号的MPSK水声信号盲解调. 舰船电子工程. 2022(04): 74-78 . 百度学术

    其他类型引用(4)

  • 加载中
图(10) / 表(1)
计量
  • 文章访问数:  1634
  • HTML全文浏览量:  623
  • PDF下载量:  59
  • 被引次数: 7
出版历程
  • 收稿日期:  2019-11-14
  • 修回日期:  2020-05-04
  • 网络出版日期:  2020-05-16
  • 刊出日期:  2020-10-13

目录

/

返回文章
返回