分享

使用PMU数据检测持续振荡的多延迟自相干方法

 GXF360 2017-12-17

使用PMU数据检测持续振荡的多延迟自相干方法

莫哈马德雷扎·戈尔巴尼帕瓦,周宁

(宾汉顿大学电子与计算机工程系, 纽约州宾汉顿, 美国)

摘 要:电力系统的振荡通常预示着稳定性问题并引起对系统稳定的疑虑。为了能够及时采取有效的补救措施,早期阶段的振荡检测非常重要。以前的文章提出了使用相量测量单元(phasor measurement unit, PMU)数据来检测持续振荡的自相干方法。通过引入多个时延来扩展自相干方法,以提高振荡检测精度。使用由简单模型和标准16 机68 节点模型生成的仿真数据来评估所提出的方法的性能。测试表明所提出的多延迟方法可以有效提高振荡检测的精度。

关键词:自助抽样法;自相干函数;相量测量单元(PMU);振荡

0 Introduction

Oscillations in power systems often herald stability problems and incur stability concerns. Oscillations with growing amplitude can cause system break-ups and even large-scale power outages[1]. Sustained oscillations can suggest equipment mis-operations and decrease the life expectancy of equipment[2]. To address the problems and concerns, it is essential to detect and analyze oscillations at their early stages.

Based on their root causes, oscillations can be divided into two categories: free oscillations and forced oscillations. Free oscillations are caused by natural interactions among different dynamic equipment. In contrast, forced oscillations are caused by external periodic perturbations. Modes and mode shapes are often used to describe free oscillations[3-4]. Many methods have been proposed to study the modes. Valuable overview of those methods can be found in[5-6]. Recently, forced oscillations are gaining more attention among researchers because forced oscillations can be a serious threat to power grid and they can even cause a system break-down[7]. Thus, there is a need to detect the forced oscillations in their early stages.

Forced oscillations can be caused by different devices in a power system. One of the early studies on forced oscillations goes back to 1966 when Ness[8] studied the effects of cyclic load variation in power systems. Sources of forced oscillations were reported as low-speed diesel generators and wind turbines in[7], nuclear accelerator in[9], steel plants, nuclear accelerators, cement mills and aluminum processing plants in[10], poorly designed power system stabilizers (PSS) in[11]. Recent studies associate with forced oscillations of wind plants and co-generating plants[12]. While there are many different sources of forced oscillations, sources are most often found on generator sides.

Many studies were carried out to detect and locate oscillations using phasor measurement unit (PMU) data[12]. To locate oscillation sources, Ma et al.[13] proposed a hybrid dynamic simulation method, which injects PMU data into a power system model to find the location of forced oscillations. The method requires an accurate dynamic model of the power system, which limits its applicability because an accurate dynamic model is often not available. Chen et al.[14] developed an energy based method, which has been successfully applied to locate oscillations with a single source. Yet, its applicability to oscillations with multiple sources is limited because of complexity of the energy flow pattern with multiple sources. Before locating oscillations, one needs to detect the oscillations. Follum et al.[15] modified the Kay’s method[16] to detect the forced oscillations under colored noise. The method focuses on estimating the frequency, amplitude, start time and end time of forced oscillations. Yet, it does not address the problem of detecting oscillations under low signal-to-noise ratio (SNRs) in a timely manner. In[17], a self-coherence method is developed to detect oscillations using one channel of data. The method is shown to work reasonably well under low SNRs while it does not fully explore the possibility of further improving detection accuracy using multiple time-delay. In[18], a bootstrap-based hypothesis test is proposed to draw a threshold on self-coherence spectra to detect an oscillation. In[19], a covariance driven stochastic subspace identification (SSI-COV)[20] method is applied to detect forced oscillations. This method can detect forced oscillations and estimate system modes. Yet, the problem of distinguishing oscillations remains unsolved. The minimum variance distortionless response (MVDR) method with derivative constraints was proposed in[21] to detect oscillations using multiple channels of data. Yet, it does not address the problem of detecting oscillation when only a single channel of data is available. A comprehensive overview on forced oscillations can be found in[22]. The aforementioned methods showed considerable accuracy in detecting and locating the forced oscillations. The studies also showed that in order to detect oscillations at their early stages with sparse PMU data, there is a need to detect the forced oscillations under low SNRs in a timely manner.

To address this need, this paper proposes a multi-delay self-coherence method to detect sustained oscillations under low SNRs. To increase the detection accuracy, this paper extends the self-coherence method[17] by using multiple time-delay signals to average out more noise. The bootstrap method is used to set up a threshold for detecting oscillations. The hypothesis testing is used to evaluate the detection accuracy of the proposed multi-delay self-coherence method.

The rest of the paper is organized as follows. Section 1 reviews the self-coherence method and extends it to a multi-delay self-coherence method.In addition, the hypothesis testing is introduced for detecting oscillations. In Section 2, the proposed method is evaluated by using simulation data generated from a simple model and the 16-machine 68-bus model. Conclusions are drawn in Section 3.

1 Coherence analysis

This section reviews the self-coherence method[17] and proposes the multi-delay self-coherence method to improve the detection accuracy. In addition, hypothesis testing is introduced to evaluate the accuracy of oscillation detection.

1.1 Review on the self-coherence method

To detect sustained oscillations, a self-coherence method was proposed in[17]. The coherence, a.k.a. magnitude squared coherence (MSC), is defined by (1) for two times-series signals (i.e., y(t) and x(t)). In (1), Pxx is the power spectral density (PSD) of x(t) while Pyy is the PSD of y(t). Pxy is the cross-spectral density of x(t) and y(t). The coherence takes a real value between 0 and 1 as in (2). It reflects how well x(t) and y(t) can be correlated at frequency f. A large MSC at frequency f0 indicates a common oscillation component at f0 shared by x(t) and y(t).

(1)

∀f∈R

(2)

To detect oscillations, [17] proposes the self-coherence method to calculate spectrum Cxx as shown in Fig. 1. As illustrated in Fig. 1, to calculate the Cxx, y(t) is set as the time delayed signal of x(t), i.e., y(t)= x(tt). Here, Δt is the time delay in seconds. Note that when there is a sustained oscillation at frequency f0, x(tt) shall be highly correlated to x(t) at f0. Therefore, there shall be a peak (close to 1) in the self-coherence spectra Cxx at f0. When noise dominates x(t), the correlation between x(tt) and x(t) often decreases with the increase of Δt. Therefore, the corresponding Cxx shall take small values (close to 0). Leveraging the above properties, the self-coherence method proposed in[17] identifies oscillations by detecting peaks in Cxx. Under the low SNRs, the self-coherence method detects oscillations by averaging out noise through the correlation between x(t) and its time-delayed signal x(tt).

Fig.1 Self-coherence spectrum Cxx(f) of signal x(t)

1.2 Multiple-delay self-coherence method

To increase detection accuracy by averaging out more noise, this paper proposes a multi-delay self-coherence method (illustrated in Fig. 2), which extends the coherence method by introducing multiple time-delay using the generalized MSC.

Denote the number of time-series signals as W. Time-delayed signals can be presented as x1(t), x2(t),…, xW(t), where


(3)

W=n+1≥2

(4)

The self-coherence spectrum of the multi-delay signals can be calculated by using the generalized MSC defined by (5) and (6) [23]

(5)

(6)

∀f∈R

(7)

In (5),operator λmax(*) indicates the largest eigenvalue. Coherence pairs of the signals with different time-delay are used to construct ∑x(f) in (6). The properties of the generalized MSC suggests that the multi-delay self-coherence spectrum defined by (5) takes a real value between 0 and 1 as in (7)[23]. If the signals are linearly correlated at frequency f0, the multi-delay self-coherence spectrum takes a value close to 1 at f0. If the signals are not correlated, the corresponding the multi-delay self-coherence spectrum takes a value close to 0.

Fig.2 Multi-delay self-coherence spectrumCx1x2…xW(f) of signal x(t)

To estimate the coherence spectrum of each pair of signals in (6), the Welch algorithm is used in this paper. To estimate the coherence of two channels using Welch algorithm, the MATLAB function mscohere is used.

Note that the multi-delay, which is used by the proposed method to increase the detection accuracy, needs more data than the method proposed in[17]. The requirement will increase detection time. Therefore, there is a trade-off between “detection accuracy” and “detection time”. Users need to choose the proper time-delay based on the needs of their applications.

1.3 Hypothesis testing

Because of the randomness of noise in power systems, this paper employs the hypothesis testing method[24] to evaluate the detection accuracy of the proposed method under a statistical framework. In general, hypothesis testing determines the probability that the observed data are obtained from an assumed distribution under a null hypothesis. In hypothesis testing, one tries to determine whether or not to reject a hypothesis.

To apply hypothesis testing in detecting oscillations,one should consider model (8). Here, A is a signal which represents the existing oscillations. Symbol z(t) is the additive noise. Symbol x(t) stands for measurements with N samples.

x(t)= A + z(t) t = 1,2,…, N

(8)

The null hypothesis (H0) is that there is no oscillation. In other words, x(t) = z(t) or equivalently A = 0. The alternative hypothesis (H1) is that random vector x(t) contains signal A. The idea can be expressed in the hypothesis testing by

H0: x(t) = z(t) H1: x(t) = A + z(t)

(9)

Generally, the self-coherence method provides a value between 0 and 1 showing how well the time-series signals are correlated to each other. To issue an oscillation alarm, a threshold needs to be set up on the coherence spectra. Because ambient noise is colored in a power system, the bootstrap method[25], which is a numerical approach, is selected to set up the threshold. The significance level α of the hypothesis test is pre-selected to set up the risk that users need to accept for issuing a false alarm. Note that there are two types of errors in hypothesis testing[24]. Type I error (i.e., false alarm) is made when H0 is true and the hypothesis testing algorithm fails to accept it. Type Ⅱ error (i.e., missing alarm) is made when H0 is false and the hypothesis testing algorithm fails to reject it. Note that a smaller α indicates a higher threshold, lower risk of false alarm, and higher risk of missing alarm.

2 Case studies

This section evaluates the performance of the proposed multi-delay self-coherence method using simulation data and compares it with the original self-coherence method[17]. Moreover, the receiver operating characteristic (ROC) curve is used to evaluate and compare their detection accuracy[16].

2.1 Simple model case

In this subsection, a simple model shown in (10) is used to generate simulation data to evaluate the proposed multi-delay self-coherence method. In (10), x(t) is the system responses, which contain a sustained oscillation at fe =13.125 Hz to mimic forced oscillation from a wind plant[12]. The Gaussian white noise e(t) is passed through a low-pass filter Gx(s) to mimic the colored ambient noise. Table I summarizes the three modes of Gx(s). The standard deviation (std) of the ambient noise is set 1.00. The amplitude of oscillation (i.e., X) is set 0.3, which makes the SNR equal to -10 dB. Sixty minutes of simulation data are generated with sampling rate of 60 samples/s. A sample time-plot of x(t) and its pure sinusoidal component are shown in Fig. 3. It can be observed that the amplitude of the sinusoidal signal is relatively low compared with x(t).

(10a)

(10b)

Fig.3 x(t) and its sinusoidal oscillation component in time domain

Table 1 Modes of the simple ambient noise model

As a benchmark, the self-coherence method[17] is applied to x(t) with its time-delay Δt1 = 5.0 s using the Welch algorithm. Sixty minutes of PMU data is split into 420 overlapping data segments, each of which has 17.1 s of data (i.e., total number of data N=1 024). The length of data segment is chosen to balance the frequency resolution and time resolution of the coherence spectrogram. Longer data segment will increase the frequency resolution but decrease time resolution. Data segment length between 10 s and 30 s was tried and no significant difference was observed. The 17.1 s was chosen to make frequency resolution equal to 0.06 Hz and the number of data equal to the next power of 2 (i.e., N=210). The resulting Cx1x2 from the original self-coherence method is shown as a heat map in Fig. 4(a). In the heat map, the amplitudes of the spectra are color coded (red: high values; white: low values).

To evaluate the proposed multi-delay self-coherence method, the three additional signals are constructed by setting the time delay as Δt1 = 5.0 s, Δt2 = 10.0 s and Δt3 = 15.0 s. The length of data segment is also set 17.1 s to facilitate comparison. The resulting Cx1x2x3x4 from the proposed multi-delay self-coherence is shown in Fig. 4(b).

It can be observed in both Fig. 4(a) and (b) that there is a solid red horizontal line at f =13.125 Hz, which corresponds to the oscillation frequency. The observation suggests that oscillations can be detected by visually inspecting the coherence on the heat map from both methods.

To compare the detection accuracy of the original self-coherence method and the proposed multi-delay self-coherence method, the 420 coherence spectra from the original self-coherence method (Cx1x2) and those from the multi-delay self-coherence method (Cx1x2x3x4) are shown in Figs. 5 and 6, respectively. In addition, their mean values and mean±3*σ are also calculated and shown in the figures. Note that following the estimation method suggested by[26], the std is estimated using (11), where Mad stands for median absolute deviation. The median values are used to reduce the negative impact of outliers and 1.482 6 is a correction factor[26].

σ=1.482 6·Mad

(11a)

(11b)

Fig.4 (a) Cx1x2 from the self-coherence method; (b) Cx1x2x3x4 from the proposed multi-delay self-coherence method (for 60 minutes of simulation data with oscillation frequency at 13.125 Hz)

Observing that the mean values of Cx1x2x3x4 are similar to those of Cx1x2. Also, observing that the std (σ) of Cx1x2x3x4 in Fig. 6 is much lower than the std (σ) of Cx1x2 in Fig. 5. The observation indicates that the detection accuracy using Cx1x2x3x4 can be higher than that using Cx1x2.

To quantify the detection accuracy, the hypothesis testing described in Section 1.3 is applied. First, the bootstrap method proposed by [18] is used to set up a threshold on the coherence spectra to issue an alarm of oscillation. When the calculated coherence spectra exceed the threshold, an alarm is issued. Probability of false alarm and probability of missing alarms are used to evaluate the accuracy of a detection algorithm. To more clearly show the detection accuracy of the proposed method, the SNR of the simple case is reduced to -25 dB. To study the detection accuracy, the significance level of the hypothesis testing is set as α = 3% for both the self-coherence method and the proposed method. The resulting probability of false alarm and missing alarm for the 60 minutes of data are listed in Table 2.

Observe in Table 2 that the probability of false alarm for both methods is very close to 3% (the setup value for threshold α). The observation indicates that α successfully sets up the probability of false alarm for the detection algorithms. Also observe that under the similar probability of false alarm (i.e., around 2.9%), the proposed method has a much lower probability of missing alarm (i.e., 0.47%) than that of the original self-coherence method (i.e., 4.3%). The observation indicates that the proposed method has a higher detection accuracy than the original self-coherence method. Note that α = 3% gives a point on the ROC curve that is picked as an example for comparison. At this point, the proposed method is compared with the original method by using the probability of missing alarm and probability of false alarm in Table 2. In addition, users can choose a value for α that fits best to their applications.

Fig.5 Overlapping Cx1x2 from the self-coherence method

Fig.6 Overlapping Cx1x2x3x4 from the proposed multi-delay self-coherence method

Table 2 Probability of false alarm and missing alarm for the cases with an oscillation at -25 db of 60 minutes of data

To study the impact of α on the detection accuracy of the proposed method, the ROC curves[16] of both the proposed method and self-coherence method are drawn in Fig. 7. The ROC curves show the relationship between the probability of detection (Pd) (i.e., 1-probability of missing alarm) and the probability of false alarm (Pf) as α varies. Fig. 7 shows that by increasing α, Pd increases at the cost of the increase of Pf for both methods. In addition, for the same value of Pf, the Pd of the proposed multi-delay self-coherence method is always higher than that of the original self-coherence method. The observation indicates that the proposed method has a higher detection accuracy than the self-coherence method for different α.

Fig.7 ROC curve for detecting oscillations under SNR = -25 dB. (Cx1x2 is from the original self-coherence method; Cx1x2x3x4 is from the proposed multi-delay self-coherence method)

2.2 16-machine 68-bus model case

To evaluate the applicability of the proposed method to a power system, the 16-machine 68-bus model[27] shown in Fig. 8 is used in this subsection for case study. The Power System Toolbox[28] is used to generate 60 minutes of simulation data. Gaussian white noise is modulated into all the load buses to simulate ambient noise. The exciter voltage reference of generator G14 is modulated by a sinusoidal signal, which has frequency of 12 Hz and magnitude of 0.1. During the modulation, SNR= -8.9 dB. To simulate the PMU measurements, the active power flow from bus 1 to bus 2 is collected at 30 samples/s.

To focus on dynamic responses, a first-order, high-pass Butterworth filter is applied to the measurement to remove its direct current components. The cut-off frequency of the filter is set as 0.01 Hz.Using the same setup as subsection 2.1, Cx1x2 from the self-coherence method as well as Cx1x2x3x4 from the multi-delay self-coherence method are calculated, and are plotted in Fig. 9 and Fig. 10, respectively. Comparing Fig. 9 and Fig. 10, one can observe that the std (σ) of Cx1x2x3x4 is smaller than that of Cx1x2. Similar to subsection 2.1, this observation indicates that for the same probability of false alarm, the proposed method can have a higher detection probability Pd than the original self-coherence method.

Fig.8 Single-line diagram of the 16-machine 68-bus model

Fig.9 Overlapping Cx1x2 from the self-coherence method by applying the Welch algorithm on the data from the 16-machine 68-bus model

Fig.10 Overlapping Cx1x23x4 from the proposed multi-delay self-coherence method by applying the Welch algorithm on the data from the 16-machine 68-bus model

3 Conclusions

This paper proposes a multi-delay self-coherence method to detect sustained oscillations in a power system using PMU data. To increase the detection accuracy, the proposed method extends the original self-coherence method by introducing multiple time-delayed signals to average out more noise. Using hypothesis testing, it is shown through simulation studies that for the same probability of false alarm, the proposed method can achieve higher probability of detection than the original self-coherence method.

4 References

[1]KOSTEREV D N, TAYLOR C W, MITTELSTADT W A. Model validation for the August 10, 1996 WSCC system outage[J]. IEEE Transactions on Power Systems, 1999, 14(3): 967-979.

[2]MAGDY M A, COOWAR F. Frequency domain analysis of power system forced oscillations[J]. IEEE Proceeding on Generation, Transmission and Distribution. 1990, 137(4): 261-268.

[3]KAMWA I, PRADHAN A, JOOS G. Robust detection and analysis of power system oscillations using the Teager-Kaiser energy operator[J]. IEEE Transactions on Power Systems, 2011, 26(1): 323-333, 2011.

[4]KORBA P. Real-time monitoring of electromechanical oscillations in power systems: first findings[J]. IET Generation, Transmission & Distribution, 2007, 1(1): 80-88.

[5]TRUDNOWSKI D J, PIERRE J W. Overview of algorithms for estimating swing modes from measured responses[C]//2009 IEEE Power & Energy Society General Meeting, Calgary: IEEE, 2009.

[6]VANFRETTI L, DOSIEK L, PIERRE J, et al. Application of ambient analysis techniques for the estimation of electromechanical oscillations from measured PMU data in four different power systems[J].European Transactions on Electrical Power, 2011, 21( 4): 1640-1656.

[7]VOURNAS C D, KRASSAS N, PAPADIAS B C. Analysis of forced oscillations in a multi-machine power system[C]//International Conference on Control, Heriot-Watt University, Edinburgh: IET, 1991.

[8]VAN N J. Response of large power systems to cyclic load variations[J]. IEEE Transactions on Power Apparatus and Systems, 1966, PAS-85(7):723-727.

[9]PINNEILO J A, VAN NESS J. Dynamic response of a large power system to a cyclic load produced by a nuclear accelerator[J].IEEE Transactions on Power Apparatus and Systems[J]. 1971, PAS-90(4): 1856-1862.

[10]RAO K, JENKINS L. Studies on power systems that are subjected to cyclic loads[J]. IEEE Transactions on Power Systems, 1988, 3(1): 31-37.

[11]MAGDY M A, COOWAR F. Frequency domain analysis of power system forced oscillations[J]. IEEE Proceedings C-Generation, Transmission and Distribution, 1990, 37(4): 261-268.

[12]SILVERSTEIN A. NASPI technical report: diagnosing equipment health and misoperations with PMU data[R]. 2015.

[13]MA J, ZHANG P, FU H, et al., Application of phasor measurement unit on locating disturbance source for low-frequency oscillation[J]. IEEE Transactions on Smart Grid, 2010, 1(3): 340-346.

[14]CHEN L, MIN Y, HU W. An energy-based method for location of power system oscillation source[J]. IEEE Transactions on Power Systems, 2013, 28(2): 828-836.

[15]FOLLUM J, PIERRE J W. Initial results in the detection and estimation of forced oscillations in power systems[C]//North American Power Symposium (NAPS), Manhatan: IEEE, 2013.

[16]KAY S M. Fundamentals of Statistical Signal Processing: Detection Theory[M]. Upper Saddle River: Prentice Hall PTR, 1998.

[17]ZHOU N, DAGLE J. Initial result in using self-coherence method for detecting sustained oscillation[J]. IEEE Transactions on Power Systems, 2015, 30(1): 522-530.

[18]GHORBANIPARVAR M, ZHOU N. Bootstrap-based hypothesis test for detecting sustained oscillations[C]// IEEE Power Eng. Soc. General Meeting, Denver: IEEE, 2015.

[19]SARMADI S A N, VENKATASUBRAMANIAN V. Inter-Area resonance in power systems from forced oscillations[J]. IEEE Transactions on Power Systems, 2015, 31(1): 378-386.

[20]SARMADI S A N, VENKATASUBRAMANIAN V. Electromechanical mode estimation using recursive adaptive stochastic subspace identification[J]. IEEE Transactions on Power Systems, 2014, 29(1): 349-358.

[21]GHORBANIPARVAR M, ZHOU N, LI X. Coherence function estimation with a derivative constraint for power grid oscillation detection[C]// IEEE GlobalSIP’2016, Washington, D.C., USA: IEEE 2016.

[22]GHORBANIPARVAR M. A survey on forced oscillations in power system[J].Journal of Modern Power System & Clean Energy, 2017(3):1-12.

[23]RAMíREZ D, VíA J, SANTAMARíA I A generalization of the magnitude squared coherence spectrum for more than two signals: definition, properties and estimation[C]//Proc. 2008 IEEE Int. Conf. Acoustics, Speech and Signal Processing (ICASSP 2008), Las Vegas: IEEE, 2008.

[24]MILI L, VAN CUTSEM T. Implementation of the hypothesis testing identification in power system state estimation[J]. IEEE Transactions on Power Systems, 1988, 3(3): 887-893.

[25]DAVISON A C, HINKLEY D V. Bootstrap methods and their application[M].Cambridge.: Cambridge Univ. Press, 1997.

[26]PHAM-GIA T, HUNG T L, The mean and median absolute deviations[J]. Mathematical and Computer Modelling, 2001, 34(7): 921-936.

[27]ROGERS G. Power System Oscillations[M]. Norwell: Kluwer, 2000.

[28]CHOW J H, CHEUNG K W. A toolbox for power system dynamics and control engineering education and research[J]. IEEE Transactions on Power Systems, 1992, 7(4): 1559-1564.

(编辑 刘文莹)

A Multi-Delay Self-Coherence Method for Detecting Sustained Oscillations Using PMU Data

GHORBANIPARVAR Mohammadreza, ZHOU Ning

(Department of Electrical and Computer Engineering, Binghamton University, Binghamton, New York, USA)

ABSTRACT:Oscillations in power systems often herald stability problems and incur stability concerns. Therefore, it is important to detect oscillations at their early stages so that effective remedial actions can be taken in time. Previously, a self-coherence method was proposed to detect sustained oscillations using phasor measurement unit (PMU) data. This paper extends the self-coherence method by introducing multiple time-delay to improve the oscillation detection accuracy. Simulation data from a simple model and the 16-machine 68-bus model are generated and used to evaluate the performance of the proposed method. Using hypothesis testing, it is shown that the proposed multi-delay method can effectively increase the accuracy of oscillation detection.

KEYWORDS:bootstrap; self-coherence; phasor measurement unit (PMU); oscillation

中图分类号:TM 72

文献标志码:A

文章编号:1000-7229(2017)09-0073-08

DOI:10.3969/j.issn.1000-7229.2017.09.011

收稿日期:2016-12-23

作者简介:

莫哈马德雷扎 · 戈尔巴尼帕瓦(1986),男,博士研究生,主要研究方向为数字信号处理、电力系统动态特性分析、振荡检测;

周宁(1970),男,博士,助理教授,博士生导师,主要研究方向为电力系统动态特性分析、小信号稳定性分析、信号处理。

    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多