一种评价风速时移性的方法与流程

文档序号:24043717发布日期:2021-02-23 17:36阅读:214来源:国知局
一种评价风速时移性的方法与流程

[0001]
本发明涉及风电场技术领域,尤其涉及一种评价风速时移性的方法。


背景技术:

[0002]
风电并网规模的快速扩大导致风电出力的随机波动性对电网调度、电力系统安全稳定运行及电能质量的影响日益彰显。由于风电机组出力与风速的三次方呈正相关关系,所以对风速的研究至关重要。如何从不同空间位置处随机波动的风速序列中挖掘出风速的时空耦合关系,对风电场(群)爬坡事件、可能出现的极端事件、降低风电并网对电网的冲击等均具有重要意义,也事关风电场效益与电网中风电的并网比例。
[0003]
现有针对不同空间位置处风速序列间时空耦合关系的研究可分为基于模型驱动的物理方法与基于数据驱动的统计方法,存在以下问题及缺陷:
[0004]
(1)物理方法无法依据实际风速的分布情况对模型进行实时优化,且模型精度与模型复杂度呈反比关系,难以适用于地形复杂的风电场与我国普遍存在的成片式开发的风电场群和风电基地。
[0005]
(2)统计方法仅实现了风速序列间相关性的量化,风速间的时移性隐含在空间分布与映射关系内,尚未实现量化。
[0006]
为了解决上述问题,本发明提供了一种基于风过程量化不同时空尺度下风速序列间时移性的方法,并提出了评价风速时移性的指标,即延迟时间dt与速度因子sf,为定量衡量不同空间位置处风速序列间的时间关系与速度变化情况提供了科学依据。所提方法的研究结果可依据实际风况进行不断完善,且适用于任何风电场与任意空间尺度,结果可作为风资源评估的有效衡量指标,为风电场出力特性的研究提供了理论依据。


技术实现要素:

[0007]
本发明的目的是提出一种评价风速时移性的方法,其特征在于,包括以下步骤:
[0008]
步骤1:数据提取及预处理;采集不同空间位置处风况信息,包括风速与风向时间序列数据,并对上述数据进行清洗与预处理;
[0009]
步骤2:时间窗口确定;考虑步骤1所获得数据的时间分辨率,确定风速时移性分析及评价的时间窗口t
w

[0010]
步骤3:起始点位p处风过程集合获取;依据时间窗口同时划分起始点位p处的风速序列与风向序列,获取风过程集合其中分别表示将起始点位p处的风序列按时间窗口t
w
划分后的第1、2、

、n个风过程;
[0011]
步骤4:终止点位q处风过程集合获取;依据点位p与点位q间的相对位置与距离和起始点位p处各风过程内的主风向与风速序列,计算起始点位p处各风过程至终止点位q处的理论延迟时间范围,据此构建终止点位q处相应时变时间窗口的风过程集合;
[0012]
步骤5:时移评价指标计算;计算起始点位p各风过程下的时间延迟dt与速度因子sf,定量评价起始点位p处风速序列至终止点位q处的时移性;
[0013]
步骤6:对所有空间点位重复步骤3至步骤5,得到可定量评价不同空间域度内各风过程时移性的指标:延迟时间dt与速度因子sf;
[0014]
步骤7:风速时移性评价;依据步骤6中计算所得风过程时移性的指标,描绘不同空间点位间风速变化的时移性特征。
[0015]
所述步骤1中的数据提取及预处理包括检验筛选、错误数据剔除及合理的插补与修正。
[0016]
所述步骤3中的风过程包含风速序列与风向序列其中分别为起始点位p处风过程内在t、t+1、

、t+t
w-1时刻的风速,分别为起始点位p处风过程内在t、t+1、

、t+t
w-1时刻的风向。
[0017]
所述步骤4中起始点位p处各风过程主风向的计算方式为:将风向以360
°
/m为间隔分为m个扇区,分别统计起始点位p处各风过程内的风向分布,进而得出各风过程内的主风向;
[0018]
理论延迟时间范围的计算方式为:若起始点位p处风过程内的主风向与点位p至点位q间的相对位置一致,则起始点位p处风过程至终止点位q的理论延迟时间范围定义如下:
[0019][0020]
其中,l
p,q
为两空间点位间的实际物理距离;分别为将起始点位p处风过程内的风速序列进行升序排列后,前a%风速数据的平均值与后b%风速数据的平均值;
[0021]
若起始点位p处风过程内的主风向与点位p至点位q间的相对位置不一致,则认为起始点位p处风过程至终止点位q的理论延迟时间为:
[0022]
所述步骤4中构建终止点位q处相应时变时间窗口的风过程集合包括:
[0023]
步骤41:依据时间窗口t
w
划分终止点位q处的风速序列,得到终止点位q在时间窗口t
w
下的风过程集合;其中终止点位q处的风过程集合可表示为风过程包含的风速序列为的风速序列为分别表示将终止点位q处的风序列按时间窗口t
w
划分后的第1、2、

、n个风过程,分别为终止点位q处风过程内在t、t+1、

、t+t
w-1时刻的风速;
[0024]
步骤42:将终止点位q在时间窗口t
w
下的各风过程分别依据相应理论延迟时间范围进行延伸,构建终止点位q处时变时间窗口的风过程集合;其中终止点位q处时变时间窗
口的风过程集合可表示为:风过程包含的风速序列为:
[0025][0026]
分别为终止点位q处时变时间窗口下的第1、2、

、n个风过程,分别为终止点位q处风过程内在内在时刻的风速。
[0027]
所述步骤5中时移评价指标的求解方法为:
[0028]
步骤51:以起始点位p处风过程内的实际风况及其至终止点位q处的理论延迟时间范围为约束条件,风过程与风过程内在时间窗口t
w
下的风速序列间欧式距离最小为优化目标,采用优化算法对所提出的两个时移性评价指标进行求解,计算起始点位p处风过程至终止点位q处的延迟时间与速度因子
[0029]
步骤52:对起始点位p处的所有风过程重复步骤51,得到起始点位p各风过程下的延迟时间dt与速度因子sf;其中,若点位风过程间的理论延迟时间为零,则认为点位风过程间的速度因子亦为零。
[0030]
所述步骤7中的风速时移性评价包括基于全部风过程的风速时移性评价、基于主风向风过程的风速时移性评价与基于典型风过程的风速时移性评价。
[0031]
基于全部风过程的风速时移性评价的方法为:对所获得的不同点位间各风过程时移性评价指标,延迟时间dt与速度因子sf进行统计分析。
[0032]
所述基于主风向风过程的风速时移性评价包括:
[0033]
步骤71:分别筛选出各起始点位处风过程集合中主风向与点位间相对位置一致的风过程;
[0034]
步骤72:对步骤71中所获得风过程对应的时移性评价指标,延迟时间dt与速度因子sf进行统计分析。
[0035]
所述基于典型风过程的风速时移性评价包括:
[0036]
步骤7a:分别采用聚类方法对各起始点位处风过程集合中主风向与点位间相对位置一致的风过程进行聚类,得到各典型风过程集合;
[0037]
步骤7b:对步骤7a中所获得各典型风过程集合所对应的时移性评价指标,延迟时间dt与速度因子sf分别进行统计分析。
[0038]
本发明的有益效果在于:
[0039]
1、实现不同空间位置处风速序列间时移性的定量分析
[0040]
提出延迟时间dt和速度因子sf两个指标定量评价不同空间位置处各风过程间的时间关系与速度变化情况,定量分析不同空间位置处风速序列间的时移特性,为风资源评
估、风电场设计运行和电力系统调度提供可靠的技术支撑。
[0041]
2、时移性结果可依据实际风况不断完善,结果科学可靠
[0042]
本发明所提供的评价风速时移性的方法属于基于数据驱动的统计方法,随着时间的推移,可用数据不断增加,所得时移性结果将依据实际风况不断完善,逐渐逼近真实情况,结果科学可靠。
[0043]
3、适用于任何风电场与任意空间尺度
[0044]
本发明所提供的风速时移性评价指标通过寻求不同空间位置处风过程间的欧氏距离最小得到,其计算方法的复杂度与地形、风电场(群)和风电基地的空间范围无关,可适用于任何风电场与任意空间尺度的风速时移性评价与分析。
附图说明
[0045]
图1为本发明风速时移性评价方法流程图;
[0046]
图2为实施例1中各空间点位相对位置分布图;
[0047]
图3为实施例1中基于全部风过程的各点位风速序列间延迟时间分布图;
[0048]
图4为实施例1中基于主风向风过程的各点位风速序列间延迟时间分布图;
[0049]
图5为实施例1中基于主风向风过程的各风速段下起始点位1#至其余点位平均延迟时间图;
[0050]
图6为实施例1中起始点位1#处6个典型风过程集合图;
[0051]
图7为实施例1中起始点位1#处各典型风过程至其余点位平均延迟时间图;
[0052]
图8为实施例1中各点位实际风速序列图;
[0053]
图9为实施例1中各点位时移后风速序列图;
[0054]
图10为实施例2中各空间点位相对位置分布图;
[0055]
图11为实施例2中基于全部风过程的各点位风速序列间时移性延迟时间分布图;
[0056]
图12为实施例2中基于全部风过程的各点位风速序列间时移性速度因子分布图;
[0057]
图13为实施例2中基于主风向风过程的各点位风速序列间时移性延迟时间分布图;
[0058]
图14为实施例2中基于主风向风过程的各点位风速序列间时移性速度因子分布图;
[0059]
图15为实施例2中基于主风向风过程的各风速段下起始点位1#至其余点位时移性平均延迟时间图;
[0060]
图16为实施例2中基于主风向风过程的各风速段下起始点位1#至其余点位时移性平均速度因子图;
[0061]
图17为实施例2中起始点位1#处6个典型风过程集合图;
[0062]
图18为实施例2中起始点位1#处各典型风过程至其余点位的时移性平均延迟时间图;
[0063]
图19为实施例2中起始点位1#处各典型风过程至其余点位的时移性平均速度因子图;
[0064]
图20为实施例2中各点位实际风速序列图;
[0065]
图21为实施例2中各点位时移后风速序列图。
具体实施方式
[0066]
本发明提出一种评价风速时移性的方法,下面结合附图和具体实施例对本发明做进一步说明。
[0067]
图1为本发明风速时移性评价的方法流程图;
[0068]
实施例1
[0069]
评价远距离风电基地间风速时移性的方法(无速度因子sf影响时),包括以下步骤:
[0070]
步骤1:数据提取及预处理;分别采集5个空间位置处的风速与风向时间序列数据,数据时间范围为2004年1月1日至2018年10月31日,时间分辨率为1h;对上述数据进行清洗与预处理,包括检验筛选、错误数据剔除及合理的插补与修正。图2为实施例1中各空间点位相对位置分布图;
[0071]
步骤2:时间窗口确定;设定时间窗口为1天,即在1天(24h)的时间窗口下对风速时移性进行分析及评价;
[0072]
步骤3:起始点位1#处风过程集合获取;依据1天的时间窗口同时划分起始点位1#处的风速序列与风向序列,获取风过程集合:其中,风过程包含风速序列:与风向序列:
[0073]
步骤4:终止点位2#处风过程集合获取;根据点位1#与点位2#间的相对位置(274
°
)与距离(103.8km)、起始点位1#处各风过程内的主风向与风速序列,计算起始点位1#处各风过程至终止点位2#处的理论延迟时间范围,据此构建终止点位2#处相应时变时间窗口的风过程集合;
[0074]
所述步骤4中起始点位1#处各风过程主风向计算方式为:将风向以360
°
/36为间隔分为36个扇区,分别统计起始点位1#处各风过程内的风向分布,进而得出各风过程内的主风向;
[0075]
所述步骤4中起始点位1#处各风过程至终止点位2#处的理论延迟时间范围计算方式为:
[0076]
若点位1#至点位2#间的相对位置(274
°
)在起始点位1#处风过程主风向(265
°-
275
°
)的范围内,则起始点位1#处风过程至终止点位2#的理论延迟时间范围为:
[0077][0078]
其中,l
1,2
为起始点位1#与终止点位2#间的实际物理距离,103.8km;分别为将起始点位1#处风过程内的风速序列进行升序排列后,前30%风速数据的平均值与后30%风速数据的平均值。
[0079]
若点位1#至点位2#间的相对位置(274
°
)不在起始点位1#处风过程主风向(45
°-
55
°
)的范围内,则起始点位1#风过程至终止点位2#的理论延迟时间为
[0080]
步骤41:依据1天的时间窗口划分终止点位2#处的风速序列,得到终止点位2#在时间
窗口1天下的风过程集合风过程内包含的风速序列
[0081]
步骤42:将终止点位2#在时间窗口1天下的各风过程分别依据相应理论延迟时间范围进行延伸,构建终止点位2#处时变时间窗口的风过程集合:风过程包含的风速序列为
[0082]
步骤5:时移评价指标计算;由于此实施例中各点位间的距离均超过100km,可认为上游风电基地基本对下游风电基地风速无影响,故仅采用延迟时间dt作为时移性评价指标;计算起始点位1#各风过程下的延迟时间dt,定量评价起始点位1#至终止点位2#处风速序列间的时移性;
[0083]
所述步骤5涉及本发明提出的一个时移性评价指标:
[0084]
时间延迟dt:表示起始点位1#处在时间窗口1天下的风过程与终止点位2#处在相应时变时间窗口下的风过程在时间轴上的关联;
[0085]
步骤51:以起始点位1#处风过程内的实际风况及其至终止点位2#处的理论延迟时间范围为约束条件,风过程与风过程内在1天时间窗口下的风速序列间欧式距离最小为优化目标,采用萤火虫优化算法对所提出的一个时移性评价指标,延迟时间dt进行求解,计算起始点位1#处风过程至终止点位2#处的延迟时间
[0086]
步骤52:对起始点位1#处的所有风过程重复步骤51,得到起始点位1#各风过程下的延迟时间dt。
[0087]
步骤6:对所有空间点位重复步骤3至步骤5(1#-3#,1#-4#,1#-5#,2#-3#,2#-4#,2#-5#,3#-4#,3#-5#,4#-5#),得到不同空间域度各风过程的延迟时间dt。
[0088]
步骤7:风速时移性评价;依据步骤6计算所得时移评价指标,描绘不同空间点位间风速变化的时移性特征。
[0089]
所述步骤7中风速时移性评价包括基于全部风过程的风速时移性评价,基于主风向风过程的风速时移性评价,与基于典型风过程的风速时移性评价。
[0090]
基于全部风过程的风速时移性评价方法为:对所获得的不同点位间各风过程时移性评价指标,时间延迟dt进行统计分析;图3为实施例1中基于全部风过程的各点位风速序列间延迟时间分布图。
[0091]
基于主风向风过程的风速时移性评价包括:
[0092]
步骤71:分别筛选出各起始点位处风过程集合中主风向与点位间相对位置一致的风过程;
[0093]
步骤72:对步骤71中所获得风过程对应的时移性评价指标,时间延迟dt进行统计分析。图4为实施例1中基于主风向风过程的各点位风速序列间延迟时间分布图,图5为实施例1中基于主风向风过程的各风速段下起始点位1#至其余点位平均延迟时间图。
[0094]
从图中可以看出,随着距离的增加,各点位风速序列间的延迟时间逐渐增长,其频率分布也更为分散。从1#点位至其余点位风速序列的延迟时间分别在0h-26h,0h-31h,0h-49h和0h-78h之间变化;各点位风速序列间的平均延迟时间与各点位间的距离和起始点位
处的平均风速具有一定的相关性,但由于风速在空间中从一个点位传播至另一个点位时,其大小会不可避免的受到地形、障碍物等各种因素的影响,其传播时间并不直接等于距离与风速的商。且风过程是由不同大小的风速所组成的时间序列,其传播规律更为复杂。以1#点位处风速序列为例,若直接将1#点位至2#-5#点位间的距离分别与1#点位处的平均风速相除,所得时间分别为3h,7h,10h与13h。而采用本发明中所提供方法进行时移性分析时,1#点位处的风速序列至2#-5#点位处延迟时间的平均值分别为2h,5h,8h与14h。这是由于即使两风过程的平均风速相同,其风况亦会存在较大差别,而不同风过程下其时移特性亦会不同,分别计算各风过程下的延迟时间更加科学合理。
[0095]
基于典型风过程的风速时移性评价包括:
[0096]
步骤7a:分别采用聚类方法对各起始点位处风过程集合中主风向与点位间相对位置一致的风过程进行聚类,得到各典型风过程集合;图6为实施例1中起始点位1#处6个典型风过程集合图。
[0097]
步骤7b:对步骤7a中所获得各典型风过程集合所对应的时移性评价指标,时间延迟dt进行统计分析。图7为实施例1中起始点位1#处各典型风过程至其余点位平均延迟时间图,表1为实施例1中起始点位1#处6个典型风过程集合至其余点位的平均延迟时间。
[0098]
表1实施例1中起始点位1#处6个典型风过程集合至其余点位的平均延迟时间
[0099][0100][0101]
基于上述所得时移性分析结果分别对点位1#处、点位2#处、点位3#处和点位4#处的风速序列进行时间轴上的变换。图8至图9为实施例1中各点位风速时移前后时间序列图。从图中可以看出,不同空间点位处的风速序列间确实存在时移特性,且采用本发明所提供的方法可很好的提取出风速时移性,实现对时移性的定量评价及分析。
[0102]
实施例2
[0103]
评价风电场内风速时移性的方法,包括以下步骤:
[0104]
步骤1:数据提取及预处理;分别采集4个空间位置处的风速与风向时间序列数据,数据时间范围为2014年1月1日至2014年12月31日,时间分辨率为1min;对上述数据进行清洗与预处理,包括检验筛选、错误数据剔除及合理的插补与修正。图10为实施例2中各空间点位相对位置分布图。
[0105]
步骤2:时间窗口确定;设定时间窗口为1h,即在1h(60min)的时间窗口下对风速时移性进行分析及评价;
[0106]
步骤3:起始点位1#处风过程集合获取;依据1h时间窗口同时划分起始点位1#处的
风速序列与风向序列,获取风过程集合:其中,风过程包含风速序列:与风向序列:
[0107]
步骤4:终止点位2#处风过程集合获取;根据点位1#与点位2#间的相对位置(136
°
)与距离(0.87km)、起始点位1#处各风过程内的主风向与风速序列,计算起始点位1#处各风过程至终止点位2#处的理论延迟时间范围,据此构建终止点位2#处相应时变时间窗口的风过程集合;
[0108]
所述步骤4中起始点位1#处各风过程主风向计算方式为:将风向以360
°
/36为间隔分为36个扇区,分别统计起始点位1#处各风过程内的风向分布,进而得出各风过程内的主风向;
[0109]
所述步骤4中起始点位1#处各风过程至终止点位2#处的理论延迟时间范围计算方式为:
[0110]
若点位1#至点位2#间的相对位置(136
°
)在起始点位1#处风过程主风向(130
°-
140
°
)的范围内,则起始点位1#处风过程至终止点位2#的理论延迟时间范围为:
[0111][0112]
其中,l
1,2
为起始点位1#与终止点位2#间的实际物理距离,0.87km;分别为将起始点位1#处风过程内的风速序列进行升序排列后,前30%风速数据的平均值与后30%风速数据的平均值。
[0113]
若点位1#至点位2#间的相对位置(136
°
)不在起始点位1#处风过程主风向(130
°-
140
°
)的范围内,则起始点位1#风过程至终止点位2#的理论延迟时间为
[0114]
步骤41:依据1h时间窗口划分终止点位2#处的风速序列,得到终止点位2#在1h时间窗口下的风过程集合风过程内包含的风速序列为
[0115]
步骤42:将终止点位2#在时间窗口1h下的各风过程分别依据相应理论延迟时间范围进行延伸,构建终止点位2#处时变时间窗口的风过程集合:风过程包含的风速序列为
[0116]
步骤5:时移评价指标计算;计算起始点位1#各风过程下的延迟时间dt和速度因子sf,定量评价起始点位1#至终止点位2#处风速序列间的时移性;
[0117]
所述步骤5涉及本发明提出的两个时移性评价指标:
[0118]
时间延迟dt:表示起始点位1#处在时间窗口1h下的风过程与终止点位2#处在相应时变时间窗口下的风过程在时间轴上的关联;
[0119]
速度因子sf:表示起始点位1#处在时间窗口1h下的风过程与终止点位2#处在相应
时变时间窗口下的风过程在速度轴上的变化情况;
[0120]
步骤51:以起始点位1#处风过程内的实际风况及其至终止点位2#处的理论延迟时间范围为约束条件,风过程与风过程内在时间窗口1h下的风速序列间欧式距离最小为优化目标,采用萤火虫优化算法对所提出的两个时移性评价指标,延迟时间dt和速度因子sf进行求解,计算起始点位1#处风过程至终止点位2#处的延迟时间与速度因子
[0121]
步骤52:对起始点位1#处的所有风过程重复步骤51,得到起始点位1#各风过程下的延迟时间dt和速度因子sf。
[0122]
其中,若点位风过程间的理论延迟时间为零,则认为点位风过程间的速度因子亦为零。
[0123]
步骤6:对所有空间点位重复步骤3至步骤5(1#-3#,1#-4#,2#-3#,2#-4#,3#-4#),得到不同空间域度各风过程的延迟时间dt和速度因子sf。
[0124]
步骤7:风速时移性评价;依据步骤6计算所得时移评价指标,描绘不同空间点位间风速变化的时移性特征。
[0125]
所述步骤7中风速时移性评价包括基于全部风过程的风速时移性评价,基于主风向风过程的风速时移性评价,与基于典型风过程的风速时移性评价。
[0126]
基于全部风过程的风速时移性评价方法为:对所获得的不同点位间各风过程时移性评价指标,时间延迟dt和速度因子sf进行统计分析;图11至图12为实施例2中基于全部风过程的各点位风速序列间时移性分布图。
[0127]
基于主风向风过程的风速时移性评价包括:
[0128]
步骤71:分别筛选出各起始点位处风过程集合中主风向与点位间相对位置一致的风过程;
[0129]
步骤72:对步骤71中所获得风过程对应的时移性评价指标,时间延迟dt和速度因子sf进行统计分析。图13至图14为实施例2中基于主风向风过程的各点位风速序列间时移性分布图。图15至图16为实施例2中基于主风向风过程的各风速段下起始点位1#至其余点位时移性图。
[0130]
从图中可以看出,对于延迟时间而言,随着距离的增加,各点位风速序列间的延迟时间逐渐增长,其频率分布也更为分散。对于速度因子而言,处于下游点位处的风速并不一定小于处于上游点位处的风速。如1#点位处的风速序列至2#-4#点位处速度因子的平均值分别为-0.62m/s,-0.16m/s和0.87m/s。这是因为此实施例中各点位所处地形为丘陵,自然风在不同空间点位处的传播受地形影响较大,且风电机组的尾流效应随着距离的增加、地形的起伏亦会有所减弱。
[0131]
对于延迟时间而言,总体呈现当风速基本相同时,距离越远,延迟时间越长;当距离相同时,风速越大,延迟时间越短。但由于受到各点位间距离较近与数据时间分辨率较低两个因素的影响,出现了当距离相同时,两点位风速序列间的延迟时间并未随着风速的增大而缩短的情况,两点位间距离越近,此情况出现的频率越高。对于速度因子而言,不同空间点位处风速序列间的速度因子随风速的变化情况有所差异,但总体呈现在大风况下自然风在不同空间点位传播时速度变化较大,在小风况下速度变化较小的规律。
[0132]
基于典型风过程的风速时移性评价包括:
[0133]
步骤7a:分别采用聚类方法对各起始点位处风过程集合中主风向与点位间相对位置一致的风过程进行聚类,得到各典型风过程集合;图17为实施例2中起始点位1#处6个典型风过程集合图。
[0134]
步骤7b:对步骤7a中所获得各典型风过程集合所对应的时移性评价指标,时间延迟dt和速度因子sf进行统计分析。图18至图19为实施例2中起始点位1#处各典型风过程至其余点位的时移性图,表2为实施例2中起始点位1#处6个典型风过程集合至其余点位的平均延迟时间与速度因子。
[0135]
表2实施例2中起始点位1#处6个典型风过程集合至其余点位的时移性
[0136][0137]
基于上述所得时移性分析结果分别对点位1#处、点位2#处和点位3#处的风速序列进行时间轴和速度轴上的变换。图20至图21为实施例2中各点位风速时移前后时间序列图。从图中可以看出,不同空间点位处的风速序列间确实存在时移特性,且采用本发明所提供的方法可很好的提取出风速时移性,实现对时移性的定量评价及分析。
[0138]
此实施例仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应该以权利要求的保护范围为准。
当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1