一种轨道交通牵引供电系统的多模态数字孪生仿真方法和装置与流程

文档序号:25298903发布日期:2021-06-04 12:29阅读:152来源:国知局
一种轨道交通牵引供电系统的多模态数字孪生仿真方法和装置与流程

1.本发明涉及轨道交通技术领域,尤其涉及一种轨道交通牵引供电系统的多模态数字孪生仿真方法和装置。


背景技术:

2.城市轨道交通牵引供电系统是整条线路安全运行的前提和保障,供电系统的安全性和可靠性关系到整个轨道交通系统的安全和正常运营,因此研究供电系统对于整个城市轨道交通的安全运行有着至关重要的意义。供电系统稳态模型可以对电压、电流、容量、能耗、潮流分布等稳态指标进行研究和优化,暂态模型可模拟各稳态之间的暂态过程,以及对谐波、纹波、短路电流等暂态特性进行分析和评估。
3.稳态仿真是指基于供电系统及负载的数学模型,通过迭代计算等手段,得到电路拓扑中各个点稳态情况下的功率、电压、电流值,并进行结果展示的一种仿真方法。稳态仿真通常不依赖于特定的仿真软件,可以采用不同的程序语言进行建模和设计,根据迭代计算的步长不同,可以对一个较长时间段系统的运行状态进行仿真。暂态仿真是建立在对供电系统及负载的完整电路模型和详细控制算法之上,对系统中各点瞬态电压、电流、纹波、谐波等实时波形和数据变化过程进行计算的一种仿真方法。由于稳、暂态仿真的要求不同,对于模型的精确程序也不相同,目前稳态和暂态模型难以在同一模型中实现,这两种仿真通常是在两种不同的仿真平台上进行。


技术实现要素:

4.本发明的实施例提供了一种城轨供电系统的多模态数字孪生建模方法和装置,提高了模型的精度。
5.一种城轨供电系统的多模态数字孪生建模方法,包括:
6.步骤1,获取预先建立的城轨供电系统的稳态模型;
7.步骤2,根据所述稳态模型,形成稳态仿真直流系统拓扑图;根据所述稳态仿真直流系统拓扑图,计算直流侧导纳矩阵;根据所述直流侧导纳矩阵以及牵引计算得到的各列车功率进行直流潮流计算,得到各牵引所的电压和功率,并将所述各牵引所的功率传给所述稳态模型的交流侧;
8.步骤3,采用牛顿

拉夫逊法以及从直流侧传入的所述牵引所的功率,对供电系统交流侧进行交流潮流计算,得到交直流公共节点上的电压和功率,并将所述公共节点上的电压传给所述稳态模型的直流侧;
9.步骤4,在所述交流侧和所述直流侧交互迭代计算,当达到预定仿真时长后,得到稳态计算步长为预定时长下的稳态交流电压电流和稳态直流电压电流
10.步骤5,获取预先建立的所述牵引供电系统的暂态模型;
11.步骤6,将稳态计算得到的所述稳态交流电压、电流和所述稳态直流电压、电流,作为计算初值传入所述暂态模型进行暂态计算,暂态计算总时长为稳态计算步长得到每一暂态步长下的暂态交流电压电流直流电压电流
12.步骤7,比较暂态步长结束时刻时稳态直流电压、电流、稳态直流电压电流与暂态交流电压、电流、暂态直流电压、电流之间是否相等;如果均相等,则输出所述暂态模型;否则,改变所述暂态模型的器件参数或者控制策略,生成更新后的暂态模型,并返回所述步骤6。
13.一种城轨供电系统的多模态数字孪生建模装置,包括:
14.第一获取单元,获取预先建立的城轨供电系统的稳态模型;
15.第一计算单元,根据所述稳态模型,形成稳态仿真直流系统拓扑图;根据所述稳态仿真直流系统拓扑图,计算直流侧导纳矩阵;根据所述直流侧导纳矩阵以及牵引计算得到的各列车功率进行直流潮流计算,得到各牵引所的电压和功率,并将所述各牵引所的功率传给所述稳态模型的交流侧;
16.第二计算单元,采用牛顿

拉夫逊法以及从直流侧传入的所述牵引所的功率,对供电系统交流侧进行交流潮流计算,得到交直流公共节点上的电压和功率,并将所述公共节点上的电压传给所述稳态模型的直流侧;
17.迭代单元,在所述交流侧和所述直流侧交互迭代计算,当达到预定仿真时长后,得到稳态计算步长为预定时长下的稳态交流电压电流和稳态直流电压电流
18.第二获取单元,获取预先建立的所述牵引供电系统的暂态模型;
19.第三计算单元,将稳态计算得到的所述稳态交流电压、电流和所述稳态直流电压、电流,作为计算初值传入所述暂态模型进行暂态计算,暂态计算总时长为稳态计算步长得到每一暂态步长下的暂态交流电压电流直流电压电流
20.比较单元,比较暂态步长结束时刻时稳态直流电压、电流、稳态直流电压电流与暂态交流电压、电流、暂态直流电压、电流之间是否相等;如果均相等,则输出所述暂态模型;否则,改变所述暂态模型的器件参数或者控制策略,生成更新后的暂态模型。
21.由上述本发明的实施例提供的技术方案可以看出,该建模方法对实际城轨供电系统得稳态和暂态建立了全面的数字孪生模型,提高了模型的精度。
22.本发明附加的方面和优点将在下面的描述中部分给出,这些将从下面的描述中变得明显,或通过本发明的实践了解到。
附图说明
23.为了更清楚地说明本发明实施例的技术方案,下面将对实施例描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他
的附图。
24.图1为本发明所述的一种城轨供电系统的多模态数字孪生建模方法的流程示意图;
25.图2为本发明提供的城轨供电系统示意图;
26.图3为本发明提供的数字孪生模型示意图;
27.图4为本发明提供的城轨供电系统稳态仿真直流系统等效模型示意图;
28.图5为本发明提供的城轨供电系统稳态仿真交流系统等效模型示意图;
29.图6为本发明提供的城轨供电系统暂态仿真模型示意图;
30.图7为本发明提供的城轨供电系统稳态、暂态联合仿真框图;
31.图8为本发明提供的稳态仿真和暂态仿真结果示意图。
具体实施方式
32.下面详细描述本发明的实施方式,所述实施方式的示例在附图中示出,其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的,仅用于解释本发明,而不能解释为对本发明的限制。
33.为便于对本发明实施例的理解,下面将结合附图以几个具体实施例为例做进一步的解释说明,且各个实施例并不构成对本发明实施例的限定。
34.如图1所示,为本发明所述的一种城轨供电系统的多模态数字孪生建模方法,包括:
35.步骤1,获取预先建立的城轨供电系统的稳态模型;
36.步骤2,根据所述稳态模型,形成稳态仿真直流系统拓扑图;根据所述稳态仿真直流系统拓扑图,计算直流侧导纳矩阵;根据所述直流侧导纳矩阵以及牵引计算得到的各列车功率进行直流潮流计算,得到各牵引所的电压和功率,并将所述各牵引所的功率传给所述稳态模型的交流侧;
37.步骤3,采用牛顿

拉夫逊法以及从直流侧传入的所述牵引所的功率,对供电系统交流侧进行交流潮流计算,得到交直流公共节点上的电压和功率,并将所述公共节点上的电压传给所述稳态模型的直流侧;
38.步骤4,在所述交流侧和所述直流侧交互迭代计算,当达到预定仿真时长后,得到稳态计算步长为预定时长下的稳态交流电压电流和稳态直流电压电流
39.步骤5,获取预先建立的所述牵引供电系统的暂态模型;
40.步骤6,将稳态计算得到的所述稳态交流电压、电流和所述稳态直流电压、电流,作为计算初值传入所述暂态模型进行暂态计算,暂态计算总时长为稳态计算步长得到每一暂态步长下的暂态交流电压电流直流电压电流
41.步骤7,比较暂态步长结束时刻时稳态直流电压、电流、稳态直流电压电流与暂态交流电压、电流、暂态直流电压、电流之间是否相等;如果均相等,则输出所述暂态模型;否则,改变所述暂态模型的器件参数或者控制策略,生成更新后的暂态模型,并返回所述步骤6。
42.其中,所述步骤2中根据所述稳态仿真直流系统拓扑图,计算直流侧导纳矩阵包括:
43.直流供电系统导纳矩阵:
[0044][0045]
含n个节点的直流潮流修正方程表示为:
[0046][0047]
式中:δp
dc
=[δp1,δp2,

δp
n
]
t
为直流功率偏差向量;u
dc
=[u
dc1
,u
dc2

u
dcn
]
t
为直流节点电压向量;δu
dc
=[δu
dc1
,δu
dc2

δu
dcn
]
t
为直流电压偏差向量;j
dc
为直流雅可比矩阵向量。
[0048]
所述根据所述直流侧导纳矩阵以及牵引计算得到的各列车功率进行直流潮流计算,得到各牵引所的电压和功率的步骤包括:
[0049]
求解直流潮流计算方程为:
[0050][0051]
所述直流潮流计算方法如下:
[0052]
步骤21,根据直流系统拓扑得到直流系统节点导纳矩阵y,设置直流系统节点电压初值列向量u
d(0)

[0053]
步骤22,求出节点功率的不平衡量δp
d(m)
=p
d(m)

p
d(n)
,δp
d(m)
为进行直流潮流计算第m次迭代的注入有功功率的计算值;
[0054]
步骤23,计算节点电压修正量向量δu
d(m)

[0055]
步骤24,对节点电压修正量δu
d(m)
进行判断,若δu
d(m)
大于限值ε1,则根据节点电压修正量对第m+1次节点电压进行重新赋值,u
d(m+1)
=u
d(m)

δu
d(m)
,并回到步骤2;否则,继续下一步骤;
[0056]
步骤25,对节点功率修正量δp
d(m)
进行判断,若δp
d(m)
大于限值ε2,则修改节点电压u
d(m+1)
,并回到步骤2;否则,继续下一步骤;
[0057]
步骤26,计算得到直流系统节点电压u
d
和功率p
d

[0058]
所述步骤3包括:
[0059]
步骤31,输入线路原始参数;
[0060]
步骤32,进行交流系统中各元件的数学建模,并得到交流系统节点导纳矩阵y;
[0061]
步骤33,设置交流侧各节点电压初始值u
i(0)
=1,δ
i(0)
=0;
[0062]
步骤34,求解交流潮流修正方程式中的不平衡量δpi
(m)
、δqi
(m)

[0063]
步骤35,求解修正方程式中雅可比矩阵中的各元素h
ij(m)
、j
ij(m)
、n
ij(m)
、l
ij(m)

[0064]
步骤36,求解修正方程式,求出各节点电压幅值和相角的修正量δu
i(m)
、δδ
i(m)

[0065]
步骤37,对节点电压幅值和相角的修正量δu
i(m)
、δδ
i(m)
进行判断;若电压幅值修正量δu
i(m)
大于ε1,相角修正量δδ
i(m)
大于ε2,则根据节点电压和相角修正量对第m+1次节点电压以及相角进行重新赋值;
[0066]
u
i(m+1)
=u
i(m)
+δu
i(m)
、δ
i(m+1)
=δ
i(m)
+δδ
i(m)
并回到步骤34;否则,继续下一步骤;
[0067]
步骤38,运用求出的各节点电压幅值和相角,继续求解交流系统的功率分布,并得到交流侧节点电压及功率的数值。
[0068]
所述修正方程式为:
[0069][0070]
方程中有功功率不平衡量δp
i
、无功功率不平衡量以及δq
i
以及各元素表达式如下:
[0071][0072]
[0073][0074][0075]
其中:δp
i
—第i个节点有功功率不平衡量;δq
i
—第i个节点无功功率不平衡量;p
i
—第i个节点有功功率;q
i
—第i个节点无功功率;u
i
—第i个节点电压;u
j
—第j个节点电压;g
ij
—ij节点之间的电导;b
ij
—ij节点之间的电纳;δ
ij
—电压u
i
与电压u
j
之间的相角差;
[0076]
本发明还提供一种城轨供电系统的多模态数字孪生建模装置,包括:
[0077]
第一获取单元,获取预先建立的城轨供电系统的稳态模型;
[0078]
第一计算单元,根据所述稳态模型,形成稳态仿真直流系统拓扑图;根据所述稳态仿真直流系统拓扑图,计算直流侧导纳矩阵;根据所述直流侧导纳矩阵以及牵引计算得到的各列车功率进行直流潮流计算,得到各牵引所的电压和功率,并将所述各牵引所的功率传给所述稳态模型的交流侧;
[0079]
第二计算单元,采用牛顿

拉夫逊法以及从直流侧传入的所述牵引所的功率,对供电系统交流侧进行交流潮流计算,得到交直流公共节点上的电压和功率,并将所述公共节点上的电压传给所述稳态模型的直流侧;
[0080]
迭代单元,在所述交流侧和所述直流侧交互迭代计算,当达到预定仿真时长后,得到稳态计算步长为预定时长下的稳态交流电压电流和稳态直流电压电流
[0081]
第二获取单元,获取预先建立的所述牵引供电系统的暂态模型;
[0082]
第三计算单元,将稳态计算得到的所述稳态交流电压、电流和所述稳态直流电压、电流,作为计算初值传入所述暂态模型进行暂态计算,暂态计算总时长为稳态计算步长得到每一暂态步长下的暂态交流电压电流直流电压电流
[0083]
比较单元,比较暂态步长结束时刻时稳态直流电压、电流、稳态直流电压电流与暂态交流电压、电流、暂态直流电压、电流之间是否相等;如果均相等,则输出所述暂态
模型;否则,改变所述暂态模型的器件参数或者控制策略,生成更新后的暂态模型。
[0084]
以下描述本发明的应用场景。
[0085]
针对现有技术中存在的缺陷,本发明提供了一种城轨供电系统多模态数字孪生建模方法,该建模方法对实际城轨供电系统得稳态和暂态建立了全面的数字孪生模型,提高了模型的精度和可靠性。
[0086]
本发明提供一种城轨供电系统多模态数字孪生建模方法,包括:
[0087]
步骤1,建立牵引供电系统稳态模型,包括:主变电所、中压环网、整流机组、能馈装置、牵引网、列车、钢轨等;
[0088]
步骤2,根据已搭建的模型形成稳态仿真直流系统拓扑图,建立各元件数学模型,依据直流系统拓扑图计算直流侧导纳矩阵,根据牵引计算得到的各列车功率进行直流潮流计算得到各牵引所的电压、功率,并将计算得到的各牵引所功率传给交流侧;
[0089]
步骤3,采用牛顿

拉夫逊法以及从直流侧传入的牵引所功率对供电系统交流侧进行交流潮流计算,得到交直流公共节点上的电压、功率,并将公共节点上的电压传给直流侧,交直流侧交互迭代计算,最终得到稳态计算步长下的稳态交流电压电流和稳态直流电压电流并传入暂态系统;
[0090]
步骤4,搭建基于供电系统及负载的完整电路模型和详细控制算法的牵引供电系统暂态模型;
[0091]
步骤5,将稳态计算得到的交流电压电流直流电压电流作为计算初值传入暂态模型进行暂态计算,暂态计算总时长为稳态计算步长得到每一暂态步长下的暂态交流电压电流和暂态直流电压电流
[0092]
步骤6,(假设在t=2

3s之间进行暂态仿真)比较3s时稳态交直流电压、电流和暂态交直流电压、电流是否相等,若相等,则结束,若不等,则改变控制参数或策略,并返回步骤5。
[0093]
所述供电系统稳态模型可以反映系统在各时刻下的电压电流及潮流分布情况。牵引供电系统潮流计算的结果可以用来对仿真对牵引供电系统的不同性能指标、牵引网及中压网络各节点电压和注入电流功率分布的合理性进行判断。
[0094]
所述稳态潮流计算通过迭代计算等手段,得到电路拓扑中各个点稳态情况下的功率、电压、电流值,暂态仿真建立在对供电系统及负载的完整电路模型和详细控制算法之上,对系统中各点瞬态电压、电流、纹波、谐波等实时波形和数据变化过程进行计算。其特征在于,包括交流潮流计算、直流潮流计算、交直流解耦计算,交流潮流计算需要通过直流潮流计算得到牵引所负荷功率,直流潮流计算需要交流潮流计算得到牵引所节点电压,采用基于牛顿拉夫逊的交直流交替迭代的方法,用直流潮流的输出作为交流潮流的输入,实现供电系统的稳态潮流计算。
[0095]
所述供电系统暂态模型可以精确模拟系统中电力电子设备的开断和控制过程及各稳态之间的暂态变化过程,通过暂态计算可以对谐波、纹波、短路电流等暂态特性进行分析和评估;
[0096]
所述暂态仿真,借助matlab的simulink仿真平台搭建供电系统全线暂态模型、设
计控制参数,进而得到一段稳态仿真离散时间点间各节点的电压/电流波形、thd、功率因数,直流纹波等数据。
[0097]
多模态电气模型中供电系统稳态模型可以对一个较长时间段供电系统的运行状态进行仿真,稳态仿真通常以1s为一个步长进行计算,因此无法反映供电系统更短时间内的暂态情况,无法计算出thd、纹波,短路电流等数据。暂态仿真是建立在对供电系统及负载的完整电路模型和详细控制算法之上,对系统中各点瞬态电压、电流、纹波、谐波等实时波形和数据变化过程进行计算。暂态仿真将稳态计算的结果作为初始值,结合暂态模型,对各个稳态之间的暂态过程进行模拟。
[0098]
所述供电系统稳态模型和供电系统暂态模型属于数字孪生多模态模型中物理模型的电气模型。数字孪生多模态电气模型就是建立一套与真实系统在外观和特性上高度对应的仿真模型,通过这套从实际系统映射而来的虚拟系统,可以准确得到当前或其他时刻实际系统的稳态和暂态指标,也可以对系统未来的状态进行预测。
[0099]
图2为本发明提供的城轨供电系统示意图,城市轨道交通牵引供电系统主要由主变电所、牵引变电所、降压所及牵引网组成,供电系统具备交流和直流两种配电形式,其中牵引所将中压环网提供的交流电转换为直流电供给列车使用。城轨供电系统模型的电气主拓扑与图1所示相同。
[0100]
图3为本发明提供的数字孪生模型示意图,所建立的数字孪生模型通过与实际系统进行大量实时数据交互以及状态上传达到同步运行的状态,是实际系统的高度仿真模型。通过大量数据交互及数据处理,数字孪生模型具备不断调整仿真模型的能力,因此提升了模型的建模精度。
[0101]
图4为本发明提供的城轨供电系统稳态仿真直流系统等效模型示意图,直流系统模型中存在牵引变电所节点和列车节点,其中将列车节点作为功率节点,列车功率可由牵引计算获得,将牵引变电所节点作为电压节点,可得直流供电系统导纳矩阵:
[0102][0103]
含n个节点的直流潮流修正方程可以表示为:
[0104][0105]
式中:δp
dc
=[δp1,δp2,

δp
n
]
t
为直流功率偏差向量;u
dc
=[u
dc1
,u
dc2

u
dcn
]
t
为直流节点电压向量;δu
dc
=[δu
dc1
,δu
dc2

δu
dcn
]
t
为直流电压偏差向量;j
dc
为直流雅可比矩阵向量
[0106]
求解直流潮流计算方程:
[0107][0108]
所述直流潮流计算方法如下:
[0109]
步骤1、根据直流系统拓扑得到直流系统节点导纳矩阵y,设置直流系统节点电压初值列向量u
d(0)

[0110]
步骤2、求出节点功率的不平衡量δp
d(m)
=p
d(m)

p
d(n)
,δp
d(m)
为进行直流潮流计算第m次迭代的注入有功功率的计算值;
[0111]
步骤3、计算节点电压修正量向量δu
d(m)

[0112]
步骤4、对节点电压修正量δu
d(m)
进行判断,若δu
d(m)
大于限值ε1,则根据节点电压修正量对第m+1次节点电压进行重新赋值,u
d(m+1)
=u
d(m)

δu
d(m)
,并回到步骤(2);
[0113]
步骤5、对节点功率修正量δp
d(m)
进行判断,若δp
d(m)
大于限值ε2,则修改节点电压u
d(m+1)
,并回到步骤(2);
[0114]
步骤6、计算得到直流系统节点电压u
d
和功率p
d

[0115]
图5为本发明提供的城轨供电系统稳态仿真交流系统等效模型,将外电源作为平衡节点,将牵引降压变电所作为pq节点,采用牛顿

拉夫逊法进行交流潮流计算,根据牛顿

拉夫逊法的数学推导公式,可得修正方程式:
[0116][0117]
所述交流潮流计算如下:
[0118]
步骤1、输入线路原始参数;
[0119]
步骤2、进行交流系统中各元件的数学建模,并得到交流系统节点导纳矩阵y;
[0120]
步骤3、设置交流侧各节点电压初始值u
i(0)
=1,δ
i(0)
=0;
[0121]
步骤4、求解交流潮流修正方程式中的不平衡量δpi
(m)
、δqi
(m)

[0122]
步骤5、求解修正方程式中雅可比矩阵中的各元素h
ij(m)
、j
ij(m)
、n
ij(m)
、l
ij(m)

[0123]
步骤6、求解修正方程式,求出各节点电压幅值和相角的修正量δu
i(m)
、δδ
i(m)

[0124]
步骤7、对节点电压幅值和相角的修正量δu
i(m)
、δδ
i(m)
进行判断,若电压幅值修正量δu
i(m)
大于ε1,相角修正量δδ
i(m)
大于ε2,则根据节点电压和相角修正量对第m+1次节点
电压以及相角进行重新赋值,u
i(m+1)
=u
i(m)
+δu
i(m)
、δ
i(m+1)
=δ
i(m)
+δδ
i(m)
并回到步骤4;
[0125]
步骤8、运用求出的各节点电压幅值和相角,继续求解交流系统的功率分布,并得到交流侧节点电压及功率的数值;
[0126]
将直流潮流计算得到的各牵引所的功率传入交流侧进行交流潮流计算,将交流潮流计算得到的交直流公共节点上的电压传入直流侧进行直流潮流计算,交直流交互迭代计算,最终输出交直流侧电压、电流,传入暂态系统。
[0127]
图6为本发明提供的城轨供电系统暂态模型中四象限变流器示意图,四象限变流器是供电系统中能馈装置的重要组成部分。
[0128]
图7为本发明提供的城轨供电系统稳态、暂态联合仿真框图,将稳态系统仿真求解出的各节点在某固定时间点下的稳态电压和电流值传入暂态系统进行暂态仿真。
[0129]
图8为本发明提供的稳态仿真和暂态仿真结果示意图,稳态仿真步长为且恰好等于暂态仿真总时长,暂态仿真步长取为且恰好等于暂态仿真总时长,暂态仿真步长取为表示t=0,1,2,3

n时刻的稳态电压;u
ts
(2),u
ts
(3)则表示在t=2s和t=3s时刻的暂态电压。假设仅在t=2

3s之间进行暂态仿真,将t=2s时刻将稳态仿真计算值传给暂态仿真,作为暂态仿真的初值,即:u
ts
(2)=u
ss
(2),暂态仿真开始运行,按照暂态仿真步长进行计算,最终并得到t=3s时刻的暂态仿真计算值,并判断u
ts
(3)=u
ss
(3)是否成立,若不成立,则改变控制方法和策略并重新进行暂态计算。
[0130]
综上,本发明涉及一种城轨供电系统多模态数字孪生建模方法,包括:供电系统稳态模型、供电系统暂态模型、稳态潮流计算、暂态仿真;本发明对城轨供电系统的稳态和暂态都建立了数字孪生模型,对系统的形态、特性、动作规律等方面进行了高精度的建模。所建立的数字孪生模型可以根据由实时数据采集和状态条件的上传不断改进模型,准确模拟当前或其他时刻实际系统的状态,也可以通过数据处理与智能算法对系统未来的故障、寿命状态及可能发生的情况进行预测。数字孪生模型突破了传统仿真模型中参数固定的限制,具备根据真实系统的实时数据采集及状态条件不断调整仿真模型的能力,既能模拟系统稳态的运行状况,也能清楚知悉各个稳态之间的暂态变化过程,使得城轨供电系统多模态数字孪生模型大大提升了传统模型的建模精度,对推动牵引供电系统的研究具有重要意义。
[0131]
以上所述,仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应该以权利要求的保护范围为准。
当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1