基于正态云模型区间证据融合的电机轴系故障风险模式辨识方法

    专利查询2022-07-08  121



    1.本发明涉及一种基于正态云模型区间证据融合的电机轴系故障风险模式辨识方法,属于机械设备故障诊断和维护技术领域。


    背景技术:

    2.随着科学技术的发展,电机轴系设备已逐步应用到一些复杂、高速以及精密的机械运动场所当中。由于机械部件之间的联系越来越密切,随着使用时长的增加,电机轴系设备不可避免会发生故障,一旦出现故障,有可能直接导致整个机械系统无法正常运行,从而影响工作效率,造成经济损失,甚至引发相关安全问题。为此利用故障风险模式辨识技术及时准确判断故障风险的发生,对于提升设备的安全性具有重要的实际意义。
    3.在对电机轴系系统进行故障风险模式辨识时,常常因为测量环境的影响以及设备本身运行状态的实时变化,从而导致电机轴系系统所产生的故障风险模式与单一故障特征之间鲜有较为明确的一一对应关系,通过传统的单传感器,单因素监测和单模型的处理与辨识方法已不能对电机轴系系统相对应的故障风险模式做出准确且高效的判断,为了增加故障风险辨识的准确性,通过多源信息融合技术,利用多个传感器采集不同故障特征信息,根据这些信息采用某种规则进行融合,从而获得对故障风险模式更为准确的描述,进而提升辨识的可靠性。


    技术实现要素:

    4.本发明的目的是提出一种基于正态云模型区间证据融合的电机轴系故障风险模式辨识方法,通过多种故障特征信息对不同故障风险模式的描述,根据正态云模型获取区间证据,利用证据推理规则将获得的不同故障特征下的区间证据进行融合,基于融合后的区间型信度,依据辨识准则得到对电机轴系故障风险模式的判定。
    5.本发明提出的基于正态云模型区间证据融合的电机轴系故障风险模式辨识方法,包括以下各步骤:
    6.(1)设定轴系的转速为ρ转/分钟,1000≤ρ≤2500,令轴系系统的故障风险模式组成的故障集合为θ={x,y,z},其中x表示电机轴系中的不平衡故障,y表示电机轴系中的不对中故障,z表示电机轴系中的基座松动故障;
    7.(2)通过分别安装在轴系支撑座水平和垂直方向上的振动加速度传感器和振动位移传感器,按设定时间间隔,连续采集时域振动序列,通过傅里叶变换获取振动加速度1~3倍频下的幅值,以及时域振动位移平均幅值,将其作为故障特征;
    8.设定轴系系统处于故障风险模式u,u∈{x,y,z}下,依次采集振动加速度1倍频、2倍频、3倍频的幅值以及时域振动位移平均幅值这4种故障特征的数据,记为其中i=1,2,3,4依次表示以上4种故障特征,n表示采集到的数据样本个数,50≤n≤500;
    9.(3)设定故障风险模式u在第i种故障特征下的样板模式云模型为构成样板模式云模型的云滴记为其中g=1,2,

    g,g代表生成云滴的总数,1000≤g≤5000,构建4种故障特征下3种故障风险模式的样板模式云模型;
    10.(4)在线获取一组待检样本数据,记为fi,其中n1≤n,表示在线数据样本总数;
    11.(5)用fi代替步骤(3)中的故障风险模式u在故障特征i下的数据,然后重复步骤(3)的整个计算过程,得到在4种故障特征下的待检模式云模型:{c1,c2,c3,c4},ci={d
    i,1
    ,...,d
    i,g
    ,...,d
    i,g
    |d
    i,g
    ={c
    i,g
    ,μ(c
    i,g
    )}};
    12.(6)在故障特征i下,计算样板模式云模型与待检模式云模型ci的匹配度,进而获得待检模式云模型ci支持3种故障风险模式发生的区间证据;
    13.(7)将步骤(6)中得到的区间证据利用证据推理规则进行融合;
    14.(8)根据步骤(7)得到h4次融合结果,并计算融合后所得q的区间证据为
    15.m(q)=[a
    θ
    ,b
    θ
    ]
    ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
    (1)
    [0016][0017][0018]
    根据式(1)-(3)则可以得到融合后的正态云模型区间证据,如表1所示:
    [0019]
    表1融合结果表
    [0020][0021][0022]
    (9)根据表1中所示区间证据给出故障风险模式辨识准则,亦即满足以下2个条件即可判定待检样本fi指向的故障风险模式为q,θ∈{x,y,z,θ}:
    [0023]

    m(θ),(θ≠θ)区间证据的左右端点分别大于其它故障风险模式对应区间证据的左右端点;
    [0024]

    m(θ)的右端点小于0.3。
    [0025]
    本发明的有益效果:
    [0026]
    一、本发明利用正态云模型描述客观对象模糊性和随机性的优势,通过构建待检模式以及故障样板模式的云模型以获得更加准确的区间证据,使得基于区间证据融合所获得的故障风险辨识结果更加符合实际情况。
    [0027]
    二、由于单个的故障特征与电机轴系故障风险模式之间没有较为清晰的映射关系,为此,本发明综合不同故障特征下的信息,采用多种信息融合的方式,依据证据推理规则进行融合,从而增加对故障风险模式辨识的准确性。
    [0028]
    三、由于单值证据对不确定性信息或模糊信息度量的不全面性、粗糙性以及度量过程可能会丢失很多有用的信息等缺点,本发明采用区间型证据,从而保证了不确定信息
    度量的完整性,使得融合辨识结果更具有全面性和可靠性。
    附图说明
    [0029]
    图1是本发明方法的流程框图。
    [0030]
    图2是本发明方法实施例中的待检模式与3种故障样板模式云模型的匹配说明图。
    [0031]
    图3是使用本发明方法的实施例中电机轴系故障风险模式辨识系统结构图。
    [0032]
    图4是本发明实施例中在“振动加速度1倍频幅值”故障特征下的待检模式与3种故障样板模式的云模型匹配图。
    [0033]
    图5是本发明实施例中在“振动加速度2倍频幅值”故障特征下的待检模式与3种故障样板模式的云模型匹配图。
    [0034]
    图6是本发明实施例中在“振动加速度3倍频幅值”故障特征下的待检模式与3种故障样板模式的云模型匹配图。
    [0035]
    图7是本发明实施例中在“时域振动位移平均幅值”故障特征下的待检模式与3种故障样板模式的云模型匹配图。
    具体实施方式
    [0036]
    本发明涉及一种基于正态云模型区间证据融合的电机轴系故障风险模式辨识方法,基于对故障典型数据的分析,在不同的故障特征下构建各种故障样板模式的云模型;根据在线监测电机轴系不同故障特征下的数据,构建故障待检模式云模型;将待检模式与各种故障样板模式的云模型进行匹配,获得该待检模式对不同故障样板模式的匹配度区间;对获取的匹配度区间进行归一化从而得到用于融合的区间证据;利用证据推理规则对区间证据进行融合,并依据一定的辨识准则,进行故障风险模式判定。本发明具体包括以下各步骤:
    [0037]
    (1)设定轴系的转速为ρ转/分钟,1000≤ρ≤2500,令轴系系统的故障风险模式组成的故障集合为θ={x,y,z},其中x表示电机轴系中的不平衡故障,y表示电机轴系中的不对中故障,z表示电机轴系中的基座松动故障。
    [0038]
    (2)通过分别安装在轴系支撑座水平和垂直方向上的振动加速度传感器和振动位移传感器,按设定时间间隔,连续采集时域振动序列,通过傅里叶变换获取振动加速度1~3倍频下的幅值,以及时域振动位移平均幅值,将其作为故障特征。
    [0039]
    设定轴系系统处于故障风险模式u,u∈{x,y,z}下,依次采集振动加速度1倍频、2倍频、3倍频的幅值以及时域振动位移平均幅值这4种故障特征的数据,记为其中i=1,2,3,4依次表示以上4种故障特征,n表示采集到的数据样本个数,50≤n≤500。
    [0040]
    (3)设定故障风险模式u在第i种故障特征下的样板模式云模型为构成样板模式云模型的云滴记为其中g=1,2,

    g,g代表生成云滴的总数,1000≤g≤5000,构建4种故障特征下3种故障风险模式的样板模式云模型。
    [0041]
    (4)在线获取一组待检样本数据,记为fi,其中n1≤n,表
    示在线数据样本总数。
    [0042]
    (5)用fi代替步骤(3)中的故障风险模式u在故障特征i下的数据,然后重复步骤(3)的整个计算过程,得到在4种故障特征下的待检模式云模型:{c1,c2,c3,c4},ci={d
    i,1
    ,...,d
    i,g
    ,...,d
    i,g
    |d
    i,g
    ={c
    i,g
    ,μ(c
    i,g
    )}}。
    [0043]
    (6)在故障特征i下,计算样板模式云模型与待检模式云模型ci的匹配度,进而获得待检模式云模型ci支持3种故障风险模式发生的区间证据。
    [0044]
    (7)将步骤(6)中得到的区间证据利用证据推理规则进行融合。
    [0045]
    (8)根据步骤(7)得到h4次融合结果,并计算融合后所得θ的区间证据为
    [0046]
    m(θ)=[a
    θ
    ,b
    θ
    ]
    ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
    (1)
    [0047][0048][0049]
    根据式(1)-(3)则可以得到融合后的正态云模型区间证据,如表2所示:
    [0050]
    表2融合结果表
    [0051][0052]
    (9)根据表2中所示区间证据给出故障风险模式辨识准则,亦即满足以下2个条件即可判定待检样本fi指向的故障风险模式为θ,θ∈{x,y,z,θ}:
    [0053]

    m(θ),(θ≠θ)区间证据的左右端点分别大于其它故障风险模式对应区间证据的左右端点;
    [0054]

    m(θ)的右端点小于0.3。
    [0055]
    进一步说,所述的1倍频为(ρ/60)hz,2倍频为(ρ/30)hz,3倍频为(ρ/15)hz。
    [0056]
    进一步说,振动加速度传感器和振动位移传感器以δt=16s为时间间隔,连续采集时域振动序列。
    [0057]
    进一步说,步骤(3)具体是:
    [0058]
    (3.1)由步骤(2)可得到故障风险模式u在故障特征i下的数据将其算数平均值作为云模型的期望值ex
    i,u
    [0059][0060]
    (3.2)设定的熵为en,0.05≤en≤0.1,超熵值为he,0.005≤he≤0.01,且
    [0061]
    (3.3)以熵en为期望,超熵he为标准差生成正态随机数然后再以ex
    i,u
    为期望,为标准差,生成正态随机数作为云滴
    [0062]
    (3.4)计算云滴对应的确定度且确定度具有概率似然性质,具体计算方式如下:
    [0063][0064]
    (3.5)令来反映云滴在二维空间中的位置信息,重复步骤(3.3)—(3.4)的计算过程g次,直至产生g个云滴以及它们各自的确定度,则可获得故障风险模式u在故障特征i下的样板模式云模型那么对于4种故障特征下的3种故障风险模式,共计可以生成4
    ×
    3=12种样板模式云模型:
    [0065]
    在1倍频幅值特征下3种故障风险模式的样板模式云模型为
    [0066]
    在2倍频幅值特征下3种故障风险模式的样板模式云模型为
    [0067]
    在3倍频幅值特征下3种故障风险模式的样板模式云模型为
    [0068]
    以及在时域振动位移平均幅值特征下3种故障风险模式的样板模式云模型为
    [0069]
    进一步说,步骤(6)具体是:
    [0070]
    (6.1)构建样板模式云模型的含熵期望曲线为待检模式云模型ci的含熵期望曲线为其中和xi为随机变量且服从高斯分布
    [0071]
    (6.2)在故障特征i下,求取3种样板模式云模型的与待检模式云模型的vi(xi)的交点其在二维空间下的坐标为具体计算方式如下:
    [0072][0073][0074]
    (6.3)由步骤(3.5)已得到在故障特征i下,3种样板模式云模型的云滴在二维空间中的坐标为设定为距离交点最近的k个云滴,且100≤k≤250,这些云滴的确定度μ(c
    i,k
    ),μ(c
    i,k
    )∈[0,1]即为待检模式云模型ci与样板模
    式云模型的匹配度,并构成匹配度集合取最小值为a
    i,u
    ,最大值为b
    i,u
    ,构成匹配度区间
    [0075]ii
    (u)=[a
    i,u
    ,b
    i,u
    ]
    ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
    (8)
    [0076]
    (6.4)将ii(u)进行归一化,得到归一化后的匹配度区间:
    [0077][0078][0079][0080][0081]
    (6.5)计算故障全集θ={x,y,z}的信度区间为
    [0082][0083][0084][0085]
    mi(θ)度量了步骤(6.4)中匹配的总体不确定度;
    [0086]
    (6.6)由式(9)和式(13)构成正态云模型区间证据,如表3所示
    [0087]
    表3通过云模型获取的区间证据表
    [0088][0089][0090]
    为了加深对待检模式云模型与3种故障风险模式的样板模式云模型进行匹配,构成区间证据的理解,这里举例说明。由于含熵期望曲线反映云模型的轮廓特征,是云滴集合的骨架,所有的云滴都在含熵期望曲线附近随机波动,则可根据含熵期望曲线相交情况来确定故障样板模式云模型与待检模式云模型的匹配程度,设在1倍频特征下电机轴系设备的3种故障风险模式的样板模式云模型及待检模式云模型c1的参数(期望、熵、超熵)如表4所示:
    [0091]
    表4构建1倍频特征下的云模型参数表
    [0092][0093]
    通过步骤(3)所述方法可得到图2所示的待检模式云模型和3种故障风险模式的样板模式云模型,通过公式(6)-(7)可得到待检模式云模型与样板模式云模型的含熵期望曲线的交点坐标即图2中“ ”所示位置,然后选取每种故障风险模式的样板模式云模型中距离对应交点最近的200个云滴,则这200个云滴所对应的确定度即为待检模式对故障样板模式的匹配度,选取其中的最大值与最小值构成匹配度区间如下表所示:
    [0094]
    表5 1倍频特征下的匹配区间表
    [0095][0096]
    根据公式(10)-(12)进行归一化,由公式(14)-(15)可得到故障全集θ的信度区间,则在1倍频特征下的区间证据如下表所示:
    [0097]
    表6 1倍频幅值特征下的证据区间表
    [0098][0099]
    根据上述相同的步骤,则可以得到2倍频幅值、3倍频幅值以及振动位移平均幅值特征下的区间证据,则在4种故障特征下的区间证据如表7所示:
    [0100]
    表7区间证据表
    [0101][0102]
    进一步说,步骤(7)具体是:
    [0103]
    (7.1)在区间θ∈{x,y,z,θ}内分别进行h(h≥2)次取点,将所取的点值记作且满足如下条件:
    [0104]

    [0105]

    [0106]
    (7.2)根据步骤(7.1)所述方法,则可得到在故障特征i下的h条点值型证据,记作ei={e
    i,1
    ,...,e
    i,h
    ,...,e
    i,h
    },其中
    [0107]
    (7.3)设定证据e
    i,h
    的可靠性因子为r
    λ
    ,满足0≤r
    λ
    ≤1,λ=1,2,3,4,重要性权重w
    λ
    =r
    λ

    [0108]
    (7.4)通过步骤(7.2)可获得4种故障特征下的4h条点值型证据,选取每种故障特
    征下的任意一条点值型证据进行组合,且每次组合中含有4个故障特征的点值型证据,则共有h4次组合,将这每次组合中的4条点值型证据利用证据推理规则进行融合,得到l个融合结果(l=1,2,

    ,h4),具体如下:
    [0109][0110][0111][0112][0113]
    其中p(θ)为θ的幂集。
    [0114]
    为了加深对区间证据通过证据推理规则融合的理解,这里举例说明。由步骤(1)-(6)可以得到区间证据,假设区间证据如表7所示,在生成的区间证据内任取2点,作为进行证据推理规则融合的点值型证据,则在4种故障特征下共有8条点值型证据,具体如下:
    [0115]
    1倍频幅值特征下的证据:
    [0116]e1,1
    ={(x,0.3953);(y,0.3864);(z,0.1972);(θ,0.0201)}
    [0117]e1,2
    ={(x,0.3959);(y,0.3872);(z,0.1986);(θ,0.0183)}
    [0118]
    2倍频幅值特征下的证据:
    [0119]e2,1
    ={(x,0.5446);(y,0.2147);(z,0.1859);(θ,0.0549)}
    [0120]e2,2
    ={(x,0.5445);(y,0.2127);(z,0.1828);(θ,0.0601)}
    [0121]
    3倍频幅值特征下的证据:
    [0122]e3,1
    ={(x,0.3526);(y,0.2825);(z,0.3560);(θ,0.0089)}
    [0123]e3,2
    ={(x,0.3529);(y,0.2821);(z,0.3559);(θ,0.0092)}
    [0124]
    时域振动位移平均幅值特征下的证据:
    [0125]e4,1
    ={(x,1.0000);(y,0.0000);(z,0.0000);(θ,0.0000)}
    [0126]e4,2
    ={(x,1.0000);(y,0.0000);(z,0.0000);(θ,0.0000)}
    [0127]
    设定可靠性因子r
    λ
    和重要性权重w
    λ
    都为0.9,将不同特征下的证据根据公式(16)-(17)进行证据推理融合,即每次都有4条证据参与融合,则共有24=16种融合结果。具体结果如表8所示:
    [0128]
    表8证据推理规则融合结果表
    [0129][0130][0131]
    根据证据推理规则融合出来的结果,选取每种故障风险模式融合后信度的最大值与最小值构成区间证据,可得融合结果如表9所示:
    [0132]
    表9区间证据结果表
    [0133][0134]
    根据步骤(9)所述辨识方法,可判定故障风险模式为x(电机轴系不平衡故障)。
    [0135]
    以下结合附图,详细介绍本发明方法的实施例:
    [0136]
    本发明方法的流程框图如图1所示,核心部分是:构建故障风险模式的样板模式云模型和待检模式云模型;将样板模式云模型和待检模式云模型进行匹配得到匹配度区间;利用全局归一化方法对匹配度区间进行处理得到区间证据;根据证据推理规则对区间证据进行融合;将融合后的区间证据依据辨识准则对故障风险模式进行辨识。
    [0137]
    以下结合图3中的电机轴系故障风险模式辨识系统的实施例,详细介绍本发明方法的各个步骤。
    [0138]
    1、电机轴系故障风险模式辨识系统的设置实例
    [0139]
    实验设备如图3中的zhs-2型多功能电机柔性轴系系统,通过分别安装在轴系支撑座水平和垂直方向上的振动加速度传感器和振动位移传感器连续采集时域振动序列,经hg-8902采集箱将信号传输至计算机,之后利用labview环境下的hg-8902数据分析软件得到轴系的振动加速度1~3倍频下的幅值,以及时域振动位移平均幅值,将其作为故障特征。
    [0140]
    2、电机轴系故障风险模式及其特征参数的选取
    [0141]
    在试验台上分别设定了故障风险模式“x为电机轴系中的不平衡”,“y为电机轴系
    中的不对中”,“z为电机轴系中的基座松动”则故障集合为θ={x,y,z};设定轴系的转速为:1500转/分钟,1倍频为25hz,n倍频为(n
    ×
    25)hz,轴系运行正常时,其各个振动加速度频率幅值都不超过0.1mm/s2,当故障发生时,不同故障风险模式所对应的频率及其幅值的增加情况也不相同,根据分析3种故障风险模式的振动能量大都集中在1倍频~3倍频上,如果仅对单一某个频率幅值的分析很难判定哪种故障风险的发生,且故障的发生是一个渐变的过程,再加上传感器易受到工作环境的影响,不能依据单独的故障特征就对故障风险结果进行辨识,故选择振动加速度1倍频~3倍频的幅值以及时域振动位移平均幅值(单位mm)作为故障特征参数。
    [0142]
    3、构建在4种故障特征下的3种故障风险模式的样板模式云模型
    [0143]
    利用本发明步骤(2)所述方法,在电机轴系上设置“x”、“y”、“z”这3种故障风险模式,通过分别安装在轴系支撑座水平和垂直方向上的振动加速度传感器和振动位移传感器,以δt=16s为时间间隔,连续采集时域振动序列,通过傅里叶变换获取振动加速度1~3倍频下的幅值,以及时域振动位移平均幅值这4种故障特征下的样本数据,其中每种故障特征下共采集到200个样本数据。依据步骤(3)所述方法,构建4种故障特征下3种故障风险模式的样板模式云模型,设置每一个云模型的云滴数为5000,共需建立12种云模型,分别如图4-图7所示,其中“x云模型”、“y云模型”、“z云模型”分别表示3种故障风险模式的样板模式云模型,构建云模型所用参数值如表10所示:
    [0144]
    表10构建故障样板模式云模型参数表
    [0145][0146]
    4、构建待检模式云模型
    [0147]
    以故障风险模式x发生为例,根据步骤(2)所述故障特征采集方法,对故障特征振动加速度1~3倍频下的幅值,以及时域振动位移平均幅值分别进行40次观测,在线获取到40个待检样本数据,依据步骤(3)所述方法,设置每一个云模型的云滴数为5000,构建4种故障特征下的4个待检模式云模型如图4-图7所示,其参数如表11所示:
    [0148]
    表11构建待检模式云模型参数表
    [0149][0150]
    5、根据待检模式和故障样板模式云模型求取区间证据
    [0151]
    根据步骤(6.1)所述方法,获取到故障样板模式以及待检模式云模型的含熵期望曲线,利用公式(6)和公式(7)得到在不同特征下的待检模式云模型与故障样板模式云模型的含熵期望曲线交点坐标如表12所示:
    [0152]
    表12待检模式与故障样板模式云模型含熵期望曲线交点坐标表
    [0153][0154]
    选取距离交点最近的200个云滴,这些云滴的确定度即为待检模式云模型ci与样板模式云模型的匹配度,取其最小值和最大值构成匹配度区间如表13所示:
    [0155]
    表13匹配度区间表
    [0156][0157]
    将获得的匹配度区间根据公式(10)-(15)进行全局归一化,得到区间证据如表14所示:
    [0158]
    表14区间证据表
    [0159][0160]
    6、利用证据推理规则对获得的区间证据进行融合得到融合后的区间证据
    [0161]
    依照本发明方法步骤(7)~(8),在生成的区间证据内任取2点,作为进行证据推理规则融合的点值型证据,则在4种故障特征下共有8条点值型证据,具体如下:
    [0162]
    1倍频特征下的证据:
    [0163]e1,1
    ={(x,0.3957);(y,0.3867);(z,0.1978);(θ,0.0198)}
    [0164]e1,2
    ={(x,0.4065);(y,0.3872);(z,0.2001);(θ,0.0062)}
    [0165]
    2倍频特征下的证据:
    [0166]e2,1
    ={(x,0.5446);(y,0.2147);(z,0.1859);(θ,0.0549)}
    [0167]e2,2
    ={(x,0.5656);(y,0.2421);(z,0.1903);(θ,0.0017)}
    [0168]
    3倍频特征下的证据:
    [0169]e3,1
    ={(x,0.3526);(y,0.2825);(z,0.3560);(θ,0.0089)}
    [0170]e3,2
    ={(x,0.3585);(y,0.2889);(z,0.3520);(θ,0.0006)}
    [0171]
    时域振动位移平均幅值特征下的证据:
    [0172]e4,1
    ={(x,1.0000);(y,0.0000);(z,0.0000);(θ,0.0000)}
    [0173]e4,2
    ={(x,1.0000);(y,0.0000);(z,0.0000);(θ,0.0000)}
    [0174]
    设定可靠性因子r
    λ
    和重要性权重w
    λ
    都为0.9,将不同特征下的证据根据公式(16)-(17)进行证据推理融合,即每次都有4条证据参与融合,则共有24=16种融合结果。根据证据推理规则融合出来的结果,选取每种故障风险模式融合后信度的最大值与最小值构成区间证据,可得到最终的融合结果如表15所示:
    [0175]
    表15融合后的证据区间表
    [0176][0177]
    7、根据融合结果进行故障风险模式辨识
    [0178]
    按照本发明方法的步骤(9)中的辨识准则,可以根据融合结果判定电机轴系故障风险为“x(电机轴系不平衡)”,与实际情况相符。
    转载请注明原文地址:https://tc.8miu.com/read-2911.html

    最新回复(0)