基于改进EWT和CMPE的配电网电能质量扰动检测方法与流程

文档序号:19126225发布日期:2019-11-13 02:09阅读:278来源:国知局
基于改进EWT和CMPE的配电网电能质量扰动检测方法与流程

本发明一种电能质量扰动信号检测技术领域,具体涉及一种基于改进ewt和cmpe的配电网电能质量扰动检测方法。



背景技术:

可再生能源由于其不确定性对电能质量的影响较大且渗透率逐年增加,在高渗透主动配电网除了出现传统的电能扰动外,由于环境气候等自然条件的约束,dg出力具有随机性,波动性和间歇性可能会导致震荡或是闪变,光伏阵列燃料电池等dg输出直流电,需经电力电子器件作为接口与系统相连,会造成大量的谐波污染,且由于dg低惯性的特征,系统更易受各类扰动影响失稳;因此对高渗透主动配电网中出现分复合扰动检测与分类有着更加重要的意义。

现有含分布式电源配电网电能质量扰动检测技术中还存在若干问题,例如:希尔伯特黄变换是经验模态分解与希尔伯特谱分析的结合,将信号置于时频域分析,但是其除了不具有可靠的理论依据外,还存在计算速度较慢、过包络等问题;局部均值分解适用于电力系统中非平稳功率信号的特征提取,在迭代筛选过程中按频率递减的规律将信号分解成一系列具有物理含义的模态函数,但是其“端点效应”与“模态混叠”的弊端十分影响结果的准确度;变分模态分解虽然将信号分解从递归筛选模式转换成非递归筛选变分模式,解决了上述局部均值分解的问题,但却是以更高的计算复杂度作为代价且依然无法自适应分解。鉴于目前的方法还有缺陷,本发明提出一种改进经验小波与复合多尺度排列熵相结合的特征提取方法,解决了上述不足。



技术实现要素:

本发明提供一种基于改进ewt和cmpe的配电网电能质量扰动检测方法,该方法将经验小波分解ewt与复合多尺度排列熵cmpe相结合,提出了一种基于改进经验小波分解的复合多尺度排列熵的瞬态特征提取方法;并搭建了ieee13节点为基础的主动配电网作为测试系统,对pq扰动信号进行分析。该方法步骤简单,分类准确,可以提高配电网的可靠性。

本发明采取的技术方案为:

基于改进ewt和cmpe的配电网电能质量扰动检测方法,包括以下步骤:

步骤1:采用改进经验小波分解ewt对主动配电网系统的pq扰动信号进行分解,滤除pq扰动信号的噪声,分解得到包含特征信息的ewt分量;

步骤2:将包含特征信息的ewt分量,作为复合多尺度排列熵cmpe算法的输入信号,利用复合多尺度排列熵cmpe算法,对各个包含特征信息的ewt分量进行排列熵计算,计算每种pq扰动信号在各个模态函数下的熵值矩阵;

步骤3:使用pca算法对计算的熵值矩阵进行降维处理,计算主元分量,得出各类pq扰动信号的特征范围;

步骤4:根据步骤3求得的降维处理后的特征值矩阵,作为svm算法的输入量;

步骤5:对含分布式能源的主动配电网系统,进行pq扰动信号识别。

搭建包含光伏、风能、电动汽车的主动配电网系统,作为检测方法有效性的测试系统。搭建的测试系统是一个ieee-13总线配电网,连接到额定功率为5mva,运行电压为4.16kv和0.48kv的电网;采样频率为3.2khz,采样时长为0.2s,通过对输出信号加入25db大小的噪声,以测试算法的鲁棒性。

所述步骤1中,主动配电网系统的pq扰动包含以下类型:

c1风能系统并网、c2风能系统中断、c3风能系统孤岛运行、c4光伏系统并网、c5光伏系统中断、c6光伏系统孤岛运行、c7电动汽车大规模离网、c8电动汽车大规模并网、c9光伏系统和风能系统同时并网、c10光伏系统和风能系统同时中断、c11光伏系统和风能系统同时孤岛运行。

所述步骤1中,采用改进经验小波变换ewt,对扰动信号进行分解,滤除扰动信号内部噪声,分解得到包含特征信息的ewt分量。所述的针对ewt改进在于构造分割准确的fft序列,以减轻间谐波、频谱泄露等因素对分解结果的影响。

所述步骤2中,将原始扰动信号进行ewt分解,其中分解得到的每个模态函数为:包含该种信号内部不同的扰动特征,分解数由信号复杂程度决定。

所述步骤3中,经pca算法降维后,不同信号扰动具有不同的主元分量范围。

所述步骤3中,分类标准,为pca降维后的特征矩阵,输入至svm分类器中即可完成自动分类,其中svm分类方法选用oao方法,核函数选用多项式核函数。

所述步骤4中,pq扰动信号的特征范围,建立在pca第一主元分量p1、第二主元分量p2、第三主元分量p3上。

所述步骤5中,所述含分布式能源的配电网故障状态,包括电压暂降、电压暂升、谐波、闪变、震荡、或凹陷尖刺;pq扰动信号包含暂升、暂降、或谐波的复合信号。

所述pq扰动信号为负序电压信号。

本发明一种基于改进ewt和cmpe的配电网电能质量扰动检测方法,技术效果如下:

1)、搭建含光伏、风能的配电网作为测试电能质量扰动检测分类算法有效性的测试系统,模拟了11类实际配电网运行情况的电能质量扰动信号。

2)、采用ewt对11类扰动信号进行了检测,在模拟扰动的情况下准确检测出扰动幅值、频率及扰动持续时间等有效信息。

3)、将包含特征信息的模态函数,作为cmpe算法的输入信号进行分类,可从ewt分量中提取不同扰动的特征矩阵。

4)、对特征矩阵使用pca降维后,即可从pca的第一主元分量、第二主元分量、第三主元分量的不同范围对不同扰动进行分类。

5)、本发明能够准确的检测分类出每一种扰动,且分类准确,具有鲁棒性,可以提高配电网的可靠性。

附图说明

图1为测试系统结构示意图。

图2为ewt分解图。

图3为ewt能量频谱。

图4为pq扰动cmpe特征图。

图4(a)为c1扰动信号cmpe特征图;

图4(b)为c2扰动信号cmpe特征图;

图4(c)为c3扰动信号cmpe特征图;

图4(d)为c4扰动信号cmpe特征图;

图4(e)为c5扰动信号cmpe特征图;

图4(f)为c6扰动信号cmpe特征图;

图4(g)为c7扰动信号cmpe特征图;

图4(h)为c8扰动信号cmpe特征图;

图4(i)为c9扰动信号cmpe特征图;

图4(j)为c10扰动信号cmpe特征图;

图4(k)为c11扰动信号cmpe特征图。

具体实施方式

基于改进ewt和cmpe的配电网电能质量扰动检测方法,包括:

搭建包含光伏、风能、电动汽车的主动配电网系统,作为检测方法有效性的测试系统。搭建的测试系统是一个ieee-13总线配电网,连接到额定功率为5mva,运行电压为4.16kv和0.48kv的电网;采样频率为3.2khz,采样时长为0.2s,通过对输出信号加入25db大小的噪声,以测试算法的鲁棒性。

该主动配电网系统中,风力机组采用双馈电机并通过变压器tr3直接并网,总容量为1.5mw;光伏发电通过三相电压型pwm变流器并网,总容量为1mw;而电动汽车采用三相桥式整流充电机与电网相连,其总负荷为0.5mw。其中,电动汽车有快速充电、机械充电与常规充电三种,快速充电的充电机功率较大,造成的电能质量问题更加严重,本发明主要考虑这种充电方式。

通过操作断路器cb1、cb2、cb3、cb4进行扰动信号模拟,其中cb1默认闭合,cb2、cb3、cb4默认断开。模拟信号有cb2闭合引起风机组并网事件c1;cb2断开引起风机组中断事件c2;cb1断开cb2闭合引起风机孤岛运行事件c3;cb3闭合引起光伏并网事件c4;cb3断开引起光伏中断事件c5;cb1断开cb3闭合引起光伏孤岛运行事件c6;cb4闭合模拟大规模电动汽车并网c7;cb5断开模拟大规模电动汽车离网c8;cb2、cb3同时闭合模拟风光同时并网事件c9;cb2、cb3同时断开模拟风光同时中断事件c10;cb1断开、cb2和cb3同时闭合模拟风光同时孤岛事件c11。测试系统结构示意图如图1所示。其电能质量扰动包含11种类型,分别是:

c1风能系统并网、c2风能系统中断、c3风能系统孤岛运行、c4光伏系统并网、c5光伏系统中断、c6光伏系统孤岛运行、c7电动汽车大规模离网、c8电动汽车大规模并网、c9光伏系统和风能系统同时并网、c10光伏系统和风能系统同时中断、c11光伏系统和风能系统同时孤岛运行。

基于改进ewt和cmpe的配电网电能质量扰动检测方法,包括以下步骤:

步骤1:采用改进经验小波分解ewt对主动配电网系统的pq扰动信号进行分解,滤除pq扰动信号的噪声,分解得到包含特征信息的ewt分量。ewt分解模态函数为自适应分解,无需手动设置。

步骤2:将包含特征信息的ewt分量,作为复合多尺度排列熵cmpe算法的输入信号,利用复合多尺度排列熵cmpe算法,对各个包含特征信息的ewt分量进行排列熵计算,计算每种pq扰动信号在各个模态函数下的熵值矩阵。

所述包含特征信息的ewt分量,根据扰动类型的不同有所区别。每个ewt分量代表了不同频率的扰动。

排列熵为对于维度为i、长度为n的时间序列x(i)、延迟时间λ、嵌入维度m的时间序列进行相空间重构,得到矩阵:

其中,1≤i≤n-(m-1)λ,x(i+λ)、x(i+(m-1)λ)为矩阵中不同位置的元素,将xi中的m元素按照升序重新排列,可得:

x(i+(j1-1)λ)≤x(i+(j2-1)λ)≤…≤x(i+(jm-1)λ)(10)

其中,j为新序列中m元素的所在列的索引,j1为第一个m元素的所在列索引,j2为第二个m素的所在列索引,jm为第m个也是最后一个m元素的所在列索引。

因此任意矢量xi都可以得出一组符号排列πn:

πn=(j1,j2,…,jm)n=1,2,…,kk≤m!(11)

其中,m为嵌入维度,j为序列中m元素的所在列索引,m!为不同符号的最大数量值,k为n的最大值,其值小于或等于m!,用p(π1),p(π2),…,p(πk),来表示每组符号排列的概率分布,则有每个πn的概率分布可由下式(12)确定:

其中,维度为i、n为的时间序列x(i)长度、λ为延迟时间、m为嵌入维度,num为取数量。

则经归一化的排列熵可由下式(13)表示:

其中,维度为i、n为的时间序列x(i)长度、λ为延迟时间、m为嵌入维度。j为新序列中m元素的所在列索引,m!为m元素最大数量。p(πj)为第j个排列概率,具体计算公式为公式(12)。

首先,在给定的时间序列{x(i),i=1,2,3,…,n}得到多尺度下的粗粒化处理结果,在给定尺度因子τ下的第k维粗粒化时间序列其中可由下式(14)得到:

其中,为在给定尺度因子τ下的第k维粗粒化时间序列,j为序列中m元素的所在列索引。把相应尺度下的cmpe取平均值则可得到在尺度因子τ下的cmpe值,即为:

其中,τ是产生粗粒化时间序列的尺度因子,m为嵌入维数,λ为延迟时间,x为输入量,为在给定尺度因子τ下的第k维粗粒化时间序列,pe的计算公式在上文公式(13)中给出。

步骤3:使用pca算法对计算的熵值矩阵进行降维处理,计算主元分量,得出各类pq扰动信号的特征范围。pq扰动信号的特征范围建立在pca第一主元分量p1、第二主元分量p2、第三主元分量p3上。该范围界定了不同扰动的数值界限,以此可以判断不同的扰动类型。

步骤4:根据步骤3求得的降维处理后的特征值矩阵,其具体计算方式如下:假设需降维的熵值矩阵为n行m列,为一n×m矩阵

则标准矩阵x*

式中,i=1,2…,n;j=1,2,…,m;xj为矩阵的第j行,sj分别为xj的均值与方差。

则相关矩阵r=x*tx*/(n-1),x*t为矩阵x*的转置,求得相关矩阵r的特征值λ1≥λ2≥…≥λm,及相应特征向量u1,u2,…,um。为确定主成分个数,还需计算方差贡献率ηi与累计方差贡献率ηsum

p的赋值取决于ηsum的大小,一般选取的p应保证ηsum大于80%及以上,本发明中p的赋值为3,贡献率ηsum皆大于90%。按照p的数值构建特征向量矩阵um×p=[u1,u2,….up],m为原始矩阵的列数。则主成分分析降维后的矩阵即为

将zn×p作为svm算法的输入量。svm算法,经模拟信号训练后结合pca降维矩阵,可以完成对电能质量扰动的自动分类,且准确度较高。

步骤5:对含分布式能源的主动配电网系统,进行pq扰动信号识别。含分布式能源的配电网故障状态,指电压暂降、电压暂升、谐波、闪变、震荡、凹陷尖刺等,以及上述复合扰动。

所述pq扰动信号为负序电压信号。在总线上捕获的电压信号通过序列分析器分析成序列元件,负序电压信号在识别运行事件的改变上较为有效的。

实施例:

本发明经验小波分解ewt检测扰动信号的实现方式为:

本发明采用改进经验小波变换,对含分布式电源的配电网的11类电能质量扰动信号c1-c11,进行预处理,首先进行降噪,提取包含特征信息的模态函数,作为复合多尺度排列熵cmpe的输入信号,从而实现电能质量分类。

经验小波分解ewt用于将非递归实值信号f(t),分解成k个具有一定稀疏性质的带限内蕴ewt分量,如下所述:

(1)、使用fft变换提取信号的主频f={fi}i=1,2,…n。其中i为fft提取的主频个数。

(2)、确定边界值ω={ωi}i=1,2,…n,i为fft提取的主频个数,n为最大个数。其将连续傅里叶频谱自适应分割为几个部分。初始边界值ω0设为0,而余下的为信号傅里叶频谱中两个相邻极小值fi、fi+1之间的中点。

(3)、确定分割区间后,经验小波定义为每个分割区间上的带通滤波器。经验小波函数φi和经验尺度函数ψi(ω)可定义如下:

式中,β(γ,ω,ωi)=β(1/2γωi(|ω|-(1-γ)ωi)),ω为输入信号频率,为信号傅里叶频谱中两个相邻极小值fi、fi+1之间的中点,其中,γ是确保两个连续状态区间之间的重叠区域最小的重要参数,其数值由计算出的边界值决定。

(4)、经过滤波后的详细因数可由快速傅里叶逆变换得到。

wx(1,n)=ifft(x(ω)φ1(ω))(23)

wx(i,n)=ifft(x(ω)ψi(ω))(24)

式中,ω为输入信号频率,i为主频次数,φi为经验小波函数和ψi(ω)经为验尺度函数,具体公式在上文(21)(22)中给出,

本发明改进经验小波变换具体的改进在于:

由于傅里叶频谱的分割直接影响ewt的结果,不同的频谱部分对应不同特定支撑频率为中心的模态。当出现电压暂降、暂升、中断等扰动信号,可能会在基波附近产生间谐波,从而导致基波分量出现频谱泄露现象,使ewt分量出现偏差。通常来讲,两组连续的谐波或是两组连续的间谐波之间的距离相较于一组谐波与一组间谐波的距离更长一些,为了更准确的分解电压电流信号,以此为基础,本发明提出一种更加准确分割频谱方法。

首先,构建原始扰动信号的基波与谐波的频域序列其中最小分量幅值不得低于基波的2%,最小频率不低于45hz,并以这组频率作为中心频率代入到公式(22)中计算ψi,其中ωi=fi-5hz,ωi+1=fi+5hz,γ=0.01。该带通滤波器的带宽bwi=10+2γfihz。谐波频谱xh(ω)可由下式(25)计算得出:

式中,xh为谐波频谱,ω为输入信号频率,i为主频次数,λ1为频域序列,φi为经验小波函数和ψi(ω)经为验尺度函数,具体公式在上文(21)(22)中给出。

从原始频谱x(ω)中减去谐波频谱xh(ω)得到只包含间谐波残留频谱xr(ω),构造最小分量幅值不得低于基波的2%,最小频率不低于5hz的间谐波序列将λ1、λ2结合并以升序排列即可得到但由于频谱泄露,一些极度不平稳信号的基波分量与整数次谐波分量中依然有分量残留,因此对λ3序列中,相邻不超过±5hz的分量分为一组,以确保实际频率仅被考虑用于滤波器设计,从而提取单频分量。经处理后的频率序列λ={fi}i=1,2,…n,(n≤m1+m2)为信号中存在的最终频域序列,间距超过10hz的频率分量都得以被准确的分割。

以下步骤描述了与ewt相关的因素。

搭建的实验系统中采样频率为3.2khz,经改进ewt分解后的除第一模态以外的每个模态都代表了不同的扰动特征。电动汽车大规模接入系统时光伏系统并网时扰动信号ewt分解如图2所示。其ewt能量频谱如图3所示,以此可以更加直观的观察其幅频特性。

本发明通过cmpe算法实现扰动信号的分类方式为:

1、将ewt分解得到包含主要特征的模态函数,作为cmpe算法的输入信号以提取不同pq扰动的特征矩阵。

排列熵是以比较时间序列相邻数值的复杂度度量算法,pe具有结构简洁、计算快、非线性单调变换不变性良好的优点,因此,pe很适合对ewt模态分解后的信号进行复杂度排列,以凸显信号特征,排列熵为对于维度为i、长度为n的时间序列x(i)、延迟时间λ、嵌入维度m的时间序列进行相空间重构,得到:

其中,1≤i≤n-(m-1)λ,x(i+λ)、x(i+(m-1)λ)为矩阵中不同位置的元素,将xi中的m元素按照升序重新排列,可得:

x(i+(j1-1)λ)≤x(i+(j2-1)λ)≤…≤x(i+(jm-1)λ)(27)

其中,j为新序列中m元素的所在列的索引,j1为第一个m元素的所在列索引,j2为第二个m素的所在列索引,jm为第m个也是最后一个m元素的所在列索引。因此任意矢量xi都可以得出一组符号排列πn:

πn=(j1,j2,…,jm)n=1,2,…,kk≤m!(28)

其中m为嵌入维度,j为序列中m元素的所在列索引,m!为不同符号的最大数量值,k为n的最大值,其值小于等于m!,用p(π1),p(π2),…,p(πk)来表示每组符号排列的概率分布,则有每个πn的概率分布可由下式(29)确定:

其中维度为i、n为的时间序列x(i)长度、λ为延迟时间、m为嵌入维度。

则经归一化的排列熵可由下式(30)表示:

其中,维度为i、n为的时间序列x(i)长度、λ为延迟时间、m为嵌入维度。j为新序列中m元素的所在列索引,p(πj)为第j个排列概率,具体计算公式为公式(29)。

首先在给定的时间序列{x(i),i=1,2,3,…,n}得到多尺度下的粗粒化处理结果,在给定尺度因子τ下的第k维粗粒化时间序列其中可由下式(31)得到:

其中为在给定尺度因子τ下的第k维粗粒化时间序列,j为序列中m元素的所在列索引。

把相应尺度下的cmpe取平均值则可得到在尺度因子τ下的cmpe值,即为:

其中τ是产生粗粒化时间序列的尺度因子,m为嵌入维数,λ为延迟时间,x为输入量,为在给定尺度因子τ下的第k维粗粒化时间序列,pe的计算公式在上文公式(30)中给出。

图4为不同pq扰动cmpe特征图。从图中可以看由不同dg所造成的不同的pq扰动的特征矩阵有所区别。其中,图4(a)的扰动为风机并网引起的电压暂降,图4(b)的扰动为风机中断引起的电压暂升,图4(c)的扰动为风机孤岛运行引起的电压暂升、暂降复合扰动,图4(d)的扰动为光伏并网引起的谐波,图4(e)的扰动为光伏中断引起的闪变,图4(f)的扰动为光伏孤岛运行引起的谐波、闪变复合扰动,图4(g)的扰动为电动汽车大规模并网引起的电压暂降与谐波复合扰动,图4(h)的扰动为电动汽车大规模离网引起的电压暂升与谐波复合扰动,图4(i)的扰动为风光系统同时并网引起的电压震荡,图4(j)的扰动为风光系统同时中断引起的脉冲,图4(k)的扰动为风光系统同时孤岛引起的电压尖刺、凹陷复合扰动;故而,可以将cmpe的熵值为pq扰动的判别因子。

2、为避免维数过多,选用pca算法对其降维。经降维后的主成分分量可以区分不同的信号特征。

降维后主元分量p1、p2、p3具体范围如表3所示。本发明设置采样点数为640个数据点,为了准确将包含不同频率的扰动信号区分开,对每种扰动的主成分分量p1、p2、p3设定不同范围,具体分类标准如表3所示,以此范围可以将复杂信号中的扰动进行逐一识别。

3、将pca降维后的矩阵输入至svm中完成扰动分类。

svm需要大量的数据进行训练,而从仿真系统中想得到如此大量的数据是十分困难的,故本发明选用ieee标准中所给出的模拟数据进行训练,每种类型的扰动生成400个信号,具体模型与信号具体参数范围(扰动时间、频率、幅值、相角),信号基频为50hz,并在±0.5hz的范围内变动。所有信号皆依照iec标准61000-4-7以200ms作为取值时长,为提升对噪声干扰信号的识别准确率,随机加入信噪比为25db~55db的噪声。svm分类方法选用oao方法,核函数选用多项式核函数。可以明显将11类配电网电能质量扰动区分。电能质量扰动信号数学模型如下所述:

①、正常信号:公式x(t)=sin(2πfft-φf),49.5≤ff≤50.5,0≤φf≤180,其中,t为时间,ff为频率,φf为相角。

②、电压中断:x(t)=[1-α(u(t-t1)-u(t-t2))]sin(2πfft-φf),其中,0.9≤α≤1,t≤t2-t1≤9t,其中,t为时间,α为电压中断幅度,u为原始信号幅值,t1为中断开始时间,t2为中断结束时间,ff为原始信号频率,φf为原始信号相角,t为周期。

③、电压暂降:x(t)=[1-α(u(t-t1)-u(t-t2))]sin(2πfft-φf),其中,0.1≤α≤0.9,t≤t2-t1≤9t,其中,t为时间,α为电压暂降幅度,u为原始信号幅值,t1为暂降开始时间,t2为暂降结束时间,ff为原始信号频率,φf为原始信号相角,t为周期。

④、电压暂升:x(t)=[1+α(u(t-t1)-u(t-t2))]sin(2πfft-φf),其中,0.1≤α≤0.8,t≤t2-t1≤9t,其中,t为时间,α为电压暂升幅度,u为原始信号幅值,t1为暂升开始时间,t2为暂升结束时间,ff为原始信号频率,φf为原始信号相角,t为周期。

⑤、谐波:x(t)=sin(2πfft-φf)+∑aisin(2πfit-φi),其中,0.03≤ai≤0.25,其中,t为时间,ff为原始信号频率,φf为原始信号相角,ai为谐波分量幅值,fi为谐波分量频率,φi为谐波分量相角。

⑥、闪变:x(t)=[1+αsin(2πβt)]sin(2πfft-φf),其中,0.1≤α≤0.2,5≤β≤25中,其中,α为闪变幅值,β为闪变频率,t为时间,ff为原始信号频率,φf为原始信号相角。

⑦、震荡:x(t)=sin(2πfft-φf)+αtexp(-(t-t1)/τ)(u(t-t1)-u(t-t2))sin(2πftt),0.1≤αt≤0.8,0.5t≤t2-t1≤3t,300≤ft≤3500,8ms≤τ≤40ms其中,αt为震荡幅值,ft为震荡频率,exp为以e为底的自然对数,τ为震荡区间,u为原始信号幅值,t为时间,t1为震荡开始时间,t2为震荡结束时间,ff为原始信号频率,φf为原始信号相角,t为周期。

⑧、凹陷:

其中,0.1≤α≤0.4,0≤t1,t2≤5t,0.01t≤t2-t1≤0.05t,n为凹陷次数,α为凹陷幅值,sign为符号函数,其功能是取某个数的符号(正或负),u为原始信号幅值,t为时间,t1为震荡开始时间,t2为震荡结束时间,ff为原始信号频率,φf为原始信号相角,t为周期。

⑨、尖刺:

其中,0.1≤α≤0.4,0≤t1,t2≤5t,0.01t≤t2-t1≤0.05t;n为尖刺次数,α为尖刺幅值,sign为符号函数,其功能是取某个数的符号(正或负),u为原始信号幅值,t为时间,t1为震荡开始时间,t2为震荡结束时间,ff为原始信号频率,φf为原始信号相角,t为周期。

表1负载配置表

表2变压器配置表

表3扰动指数分类表

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