電力網絡的諧振穩定性分析方法研究

徐政，王世佳，邢法財，肖晃慶

(浙江大學電氣工程學院，杭州市 310027)

電力網絡的諧振穩定性分析方法研究

徐政，王世佳，邢法財，肖晃慶

(浙江大學電氣工程學院，杭州市 310027)

隨著電力系統電力電子化程度的不斷加深，近年來出現了多起機理不明的新的振蕩現象。提出了電力網絡諧振穩定性的概念，試圖將上述機理不明的振蕩現象納入到電力網絡的諧振不穩定范疇，從而基于線性網絡理論對眾多復雜振蕩現象進行分析，為借助數學上完全成熟的線性系統理論解決電力系統實際問題提供一條途徑。圍繞判斷諧振穩定性的分析方法展開，通過引入電力網絡的s域節點導納矩陣，將電力網絡的諧振穩定性問題歸結為判斷s域節點導納矩陣行列式的零點在復平面上的分布問題。首先，理論證明了s域節點導納矩陣行列式的零點就是系統的特征值。其次，給出了求解s域節點導納矩陣行列式零點的實部-虛部交叉迭代法。接著，推導了特定諧振模式下的節點電壓振型和參與因子矩陣，這2個指標可用來定位特定諧振模式發生的位置。最后，通過算例展示了所提方法在分析電力網絡諧振穩定性方面的有效性。

電力網絡；狀態空間模型；s域節點導納矩陣；諧振模式；諧振穩定性；節點電壓振型；參與因子矩陣

0 引 言

振蕩是電力系統運行過程中的一種常見現象[1-2]。一般來說，電力系統振蕩可以分為3種類型：(1)發電機軸系的振蕩(扭振)；(2)發電機轉子之間的振蕩；(3)電力網絡內部的固有諧振。隨著大量電力電子裝置接入電力系統，電力網絡的固有諧振問題變得更加嚴重和復雜。由于某些電力電子裝置在一定的頻段內存在負電阻效應，使得原先諧振穩定的電力網絡有可能變得諧振不穩定，從而造成嚴重的后果。例如，2011年至今，河北沽源地區的雙饋風電場與其電容串補裝置相互作用，引發了上百次頻率為 3～10 Hz的次同步振蕩，造成變壓器異常振動和大量風機脫網[3]。福建廈門柔性直流輸電工程在直流側出現過次同步振蕩問題，該振蕩現象為非衰減和非等幅振蕩，振蕩頻率為25 Hz左右[4]。新疆哈密“7·1”直驅風機引起的火電機組次同步振蕩事件[5]，可以認為是由電力網絡諧振與發電機軸系扭振相互作用而導致的機網復合共振問題。

當電力網絡發生諧振時，很容易引起系統的過電壓、過電流，甚至會與發電機軸系扭振相互作用而導致機網復合共振，威脅電力系統的安全穩定運行。因此，亟需對電力網絡的固有諧振結構和諧振穩定性進行透徹的研究。

本文深入研究基于s域節點導納矩陣Y(s)的電力網絡諧振穩定性分析方法。首先，在理論上論證s域節點導納矩陣Y(s)行列式為0的根就是系統的特征值，并用一個簡單算例進行了驗證。其次，提出一種實部-虛部交叉迭代的方法，用以確定s域節點導納矩陣Y(s)行列式為0的根的位置；這種求解方法與采用Newton-Raphson法迭代求解det[Y(s)]零點的方法相比，可以大幅度減少計算量。接著，推導特定諧振模式下的節點電壓振型和參與因子矩陣，這2個指標可用來定位特定諧振模式發生的位置。最后，通過IEEE次同步諧振第一標準測試系統和IEEE 39節點測試系統，展示對實際電網進行諧振穩定性分析的結果。

1 理論基礎

1.1 諧振穩定性的定義

當考慮輸電線路等分布參數元件時，因描述分布參數元件特性的方程是偏微分方程，整個電力網絡已不能用線性定常系統的標準狀態空間模型來描述。另外，若進一步考慮元件參數隨頻率而變化的特性，那么即使對于由集總參數元件構成的電力網絡，也無法用線性定常系統的標準狀態空間模型來描述。因此，當電力網絡包含分布參數元件和頻變參數元件時，所謂的電力網絡諧振穩定性該如何定義就成為一個重要問題。

1999年，文獻[7]首先提出了采用s域節點導納矩陣Y(s)來分析復雜電力網絡小信號穩定性的概念和方法，這里所謂的“復雜電力網絡小信號穩定性”本質上與本文所稱的“電力網絡諧振穩定性”相一致。因此，對于包含分布參數元件和頻變參數元件的復雜電力網絡，諧振穩定性的概念將基于s域節點導納矩陣Y(s)來定義。2001年，文獻[8]對采用s域節點導納矩陣Y(s)進行小信號穩定性分析的方法做了進一步的發展和完善。其要點如下：(1)所謂的“s域節點導納矩陣Y(s)”，國內也稱為“運算導納矩陣”，是早已存在的概念；例如，電容C的運算導納是sC，電感L的運算導納是1/(sL)；簡單地說，將交流穩態分析時元件導納模型中的j用s來替換就構成了對應元件的運算導納，對于分布參數的輸電線路，也有類似的結果[8]；在得到各元件的運算導納模型后，構建運算導納矩陣的步驟與交流穩態分析時構建節點導納矩陣的步驟完全一致。(2)對于包含分布參數元件和頻變參數元件的電力網絡，構建其s域節點導納矩陣Y(s)并不存在特殊困難；因而對于一般性的電力網絡，基于s域節點導納矩陣Y(s)分析其諧振穩定性具有普遍的適用性。(3)令Y(s)的行列式為det[Y(s)]，那么det[Y(s)]=0的根(以下統稱為det[Y(s)]的零點)就是該電力網絡的諧振模式，det[Y(s)]的所有零點就是該電力網絡的所有諧振模式，如果det[Y(s)]的所有零點都位于復平面的左半平面，那么該電力網絡就是諧振穩定的。容易證明，包含分布參數元件的電力網絡，其諧振模式有無限個。

由上面的介紹可知，判斷電力網絡諧振穩定性與判斷det[Y(s)]的所有零點是否都位于復平面的左半平面等價。因此，判斷電力網絡諧振穩定性的最直接方法就是求出det[Y(s)]的所有零點，或者求出det[Y(s)]在指定頻段內的所有零點。更進一步，可以對電力網絡的諧振結構進行分析，電力網絡的諧振結構包含4方面的信息：(1)電力網絡在分析的頻段內存在哪些固有諧振模式；(2)各諧振模式的頻率；(3)各諧振模式的阻尼；(4)各諧振模式的振型。文獻[8]提出了采用Newton-Raphson法迭代求解det[Y(s)]零點的方法。該方法需要同時建立2個s域節點導納矩陣，一個是Y(s)，一個是Y(s)關于s的導數矩陣dY(s)/ds，計算量比較大。

1.2 系統特征值的定義

電力網絡的動態特性可以通過KVL電壓方程、KCL電流方程和儲能元件動態特性方程進行描述；其中前2類方程為代數方程，后1類方程為一階微分方程。當電力網絡由集總參數元件構成，且元件參數不隨頻率變化時，根據上述3類方程可以列寫出單輸入單輸出系統如下的狀態空間模型[9-10]：

(1)

式中：{A,T}共同構成狀態空間的系統矩陣；x為狀態變量向量，包括儲能元件狀態(電容電壓、電感電流等)；輸出y為第k個節點的電壓；輸入u為第j個節點的注入電流；b是一維常數列向量；c是一維常數行向量。

對式(1)進行拉氏變換，有

(2)

因此s域下轉移阻抗zkj(s)的表達式為

(3)

式中：分母det(sT-A)表示矩陣sT-A的行列式，“*”表示相應矩陣的伴隨矩陣。根據線性代數的基本理論可知，det(sT-A)是系統的特征多項式，det(sT-A)=0的根是系統的特征值。系統特征值對應于電力網絡的所有諧振模式；當特征值都位于復平面的左半平面時，該電力網絡就是諧振穩定的。由于狀態空間模型需要在集總參數元件的條件下列寫，此時諧振模式為有限個。

1.3 基于s域節點導納矩陣求取特征值

除狀態空間模型外，可以采用s域節點導納矩陣分析系統特征值[7-8]。s域下系統節點電壓方程為

Y(s)Vnode(s)=Inode(s)

(4)

式中：Y(s)為s域下的節點導納矩陣；Vnode(s)和Inode(s)分別為s域下的節點電壓向量和節點注入電流向量。同樣考慮單輸入單輸出系統，設在節點j上的注入電流為ij(s)，在節點k上的輸出電壓為vk(s)。則方程(4)可以寫成如下形式：

(5)

這里b′是一維列向量，其第j個元素為1，其余元素為0；而c′是一維行向量，其第k個元素為1，其余元素為0。因此s域下轉移阻抗zkj(s)的表達式為

(6)

式中：det[Y(s)]為Y(s)的行列式；Y(s)*為Y(s)的伴隨矩陣。

由于Y(s)的元素不一定都是s的多項式，例如對于電感元件L，其運算導納的表達式是1/(sL)，因此通過對det[Y(s)]外乘s的冪函數sm(m為正整數)，一定能夠將det[Y(s)]化成s的多項式。對比式(3)和式(6)，det(sT-A)與smdet[Y(s)]都是系統的特征多項式。這樣，在s≠0的條件下，det(sT-A)=0與det[Y(s)]=0的根是完全一致的。因此，det[Y(s)]=0的非零根一定是系統的特征值，從而論證了采用狀態空間法和s域節點導納矩陣法在求取系統特征值上具有一致性。

值得指出的是，在構建s域節點導納矩陣Y(s)時，并不要求網絡元件使用集總參數模型。因此對于一般性的電力網絡(包括分布參數元件和頻變參數元件)，相對于狀態空間法，基于s域節點導納矩陣分析諧振穩定性具有更為普遍的適用性。容易證明，包含分布參數元件的電力網絡，其諧振模式有無限個。

1.4 實例驗證

圖1 RLC串聯電路Fig.1 RLC series circuit

圖1所示RLC電路的狀態空間方程為

(7)

進而可以求出系統的特征方程為

(8)

若采用s域節點導納矩陣法，則該電路的s域節點導納矩陣為：

(9)

而

(10)

顯然，在s≠0的條件下，式(8)和式(10)具有相同的零點。這說明采用狀態空間法或s域節點導納矩陣法進行電力網絡的諧振穩定性分析，所得結果是完全一致的。

2 基于Y(s)的特征值求解算法

以上分析可知，電力網絡特征值求取問題可轉化為求解s域節點導納矩陣行列式等于0的根，即求解如下方程的根：

det[Y(s)]=0

(11)

由于s是復平面上的變量，直接求解上述方程是一件非常困難的事情。現有文獻主要采用Newton-Raphson法迭代求解，需要求取Y(s)關于s的導數矩陣[8]。為了能夠快速準確地求出det(Y(s))的零點，本文提出實部-虛部交叉迭代的求解方法。

2.1 算法原理

當Y(s)行列式不等于0時，其行列式的倒數為一個有限值；當Y(s)行列式趨近于0時，其行列式的倒數趨向于無窮大。故定義hm(s)為：

(12)

hm(s)在某個特征值附近的示意圖如圖2所示。需要注意的是，該曲面的頂點在復平面(σ-jω平面)上的投影就是相應的特征值。

當s的實部固定為σi，虛部在一定范圍內變化時，如圖2所示，相當于平面σ=σi與曲面hm(s)相交的曲線；找出該曲線極大值所對應的虛部坐標，記為jωi。接下來，將s的虛部固定為jωi，而實部在一定范圍內變化，如圖3所示，相當于平面jω=jωi與曲面hm(s)的交線；找出這條交線的極大值點對應的實部坐標，記為σi+1。

重復上述過程，直到所得的實部和虛部不再改變，或者達到所設定的精度要求為止。最后求出的實部和虛部的組合就是s域下系統的特征值。在電力網絡諧振穩定性分析時，系統特征值又稱為諧振模式。因此一旦求出系統特征值，則對應此諧振模式的阻尼和頻率就是已知的。

圖2 平面σ=σi與曲面相交情況Fig.2 Intersection of σ=σi and tapered surface

圖3 平面jω=jωi與曲面相交情況Fig.3 Intersection of jω=jωi and tapered surface

2.2 算法流程圖

上述的求解思路只是針對1個特征值，對含有多個特征值的系統，其求解過程類似。電力網絡特征值的詳細求解流程如圖4所示。

3 諧振模式的節點電壓振型和參與因子矩陣

借鑒文獻[11]的做法，可以定義諧振模式sk的節點電壓振型和參與因子矩陣。設電力網絡的第k個諧振模式為sk，則det[Y(sk)]=0。由矩陣理論知，矩陣行列式的值等于矩陣所有特征值的乘積，因此，復常數矩陣Y(sk)必有一個零特征值=0。設Y(sk)為n階

圖4 電力網絡諧振模式的求解流程圖Fig.4 Solving process of electric network resonance modes

Y(sk)=RΛR-1=RΛL

(13)

式中：Λ=diag(λ1,λ2,…,λn)是對角元素為特征值的對角矩陣，不妨設1=0；R=[R1,R2,…,Rn]是Y(sk)的右特征向量矩陣；L=R-1=T是Y(sk)的左特征向量矩陣。且由于Y(sk)為對稱矩陣，因此又有R-1=RT。根據s域節點導納方程式(4)，有

Y(sk)Vnode(sk)=Inode(sk)

(14)

令

Umode=LVnode(sk)

(15)

Jmode=LInode(sk)

(16)

則式(14)可以變換為

(17)

即

(18)

(19)

因此，我們定義R1為對應諧振模式sk的節點電壓振型，表示在諧振模式sk下電網中各節點電壓的相對大小和相位。而根據式(16)，有

J1=L1Inode(sk)

(20)

因此，根據式(18)～(20)有

(21)

(22)

因此定義矩陣

(23)

為參與因子矩陣，P的元素pij表示在諧振模式sk下電網中節點j的注入電流對節點i電壓的相對作用大小。

4 算例及結果分析

4.1 IEEE次同步諧振第一標準測試系統分析

IEEE次同步諧振第一標準測試系統結構如圖5所示，發電機參數見文獻[12]。發電機的阻抗頻率特性采用測試信號法[13]計算得出，如圖6所示。

圖5 IEEE 次同步振蕩第一標準測試系統接線圖Fig.5 Single line diagram of IEEE subsynchronousresonance (SSR) first benchmark model

對該測試系統在5～115 Hz頻率范圍進行固有諧振結構分析，發現存在1個次同步諧振模式，頻率fssr=39.3 Hz，衰減因子σssr=-0.454 51/s。由于衰減因子為負，說明該測試系統是諧振不穩定的。該諧振模式的節點電壓振型如圖7所示，參與因子矩陣如式(24)所示。由節點電壓振型圖和參與因子列表可見，模式SSR的主要參與節點是3號節點，諧振類型是全局性的。

圖6 發電機的阻抗頻率特性Fig.6 Impedance-frequency characteristicsof generator

(24)

圖7 模式SSR的節點電壓振型圖Fig.7 Nodal voltage mode shape of mode SSR

4.2 IEEE 39節點測試系統分析

IEEE 39節點測試系統[14]如圖8所示。對該測試系統在0～1 500 Hz頻率范圍進行固有諧振結構分析，結果如表1所示。表1顯示該系統在0～1 500 Hz頻率范圍存在19個諧振模式，所有諧振模式的衰減因子均大于0，說明系統是諧振穩定的。下面對諧振模式1(131.6 Hz)進行詳細分析。模式1 (131.6 Hz) 的節點電壓振型圖如圖9所示，數值大于0.1的參與因子矩陣元素只有1個，為p28,28=0.104 1。

由節點電壓振型圖和參與因子矩陣可見，模式1的主要參與節點是28號和26號、29號節點，諧振類型為局部諧振，主要是28號節點的無功負荷和所連線路的對地電容之間的諧振。

圖8 IEEE 39節點測試系統Fig.8 IEEE 39-bus test system

4.3 加入風電場的IEEE 39節點測試系統分析

在IEEE 39節點標準測試系統上增加1個風電場節點，考慮風電場通過遠距離輸電線路加裝50%串補接入到IEEE 39節點標準測試系統的25號節點上，如圖10所示。

由雙饋風力發電機聚合而成的風電場額定容量 1 500 MW，風機阻抗模型采用文獻[15]給出的模型，聚合后風電場在0～100 Hz范圍內的阻抗頻率特性如圖11所示，顯示風電場在25～45 Hz頻率范圍內呈現負電阻-電感特性(相頻特性大于90°)。

表1IEEE39節點測試系統的諧振模式
Table1ResonancemodesofIEEE39-bustestsystem

圖9 模式1(131.6 Hz)的節點電壓振型圖Fig.9 Nodal voltage mode shape of mode 1 (131.6 Hz)

圖10 加入到IEEE 39節點標準測試系統上的風電場Fig.10 Wind farm connected to IEEE 39-bus test system

對該系統在0～1 500 Hz頻率范圍進行固有諧振結構分析，結果如表2所示。表2與表1對比顯示，雙饋風機通過串補系統接入原系統25號節點后，部分諧振模式發生了變化，如新增了29.8 Hz的諧振模式，而且29.8 Hz的諧振模式衰減因子為負(發散)，說明雙饋風機通過串補送出可以引起不穩定的次同步諧振。進一步對不穩定的29.8 Hz諧振模式進行詳細分析，得到模式1(29.8 Hz)的節點電壓振

圖11 風機阻抗頻率特性Fig.10 Impedance-frequency characteristic ofwind generator

圖12 模式1(29.8 Hz)的節點電壓振型圖Fig.12 Nodal voltage mode shape of mode 1 (29.8 Hz)

5 結 論

本文通過構建電力網絡的s域節點導納矩陣，提出了電力網絡諧振穩定性的概念，并提出了相應的穩定判據——s域節點導納矩陣行列式的零點就是系統的特征值，若零點均位于左半平面，則電力網絡是穩定的。這一穩定判據對含頻變元件和分布參數元件的電力網絡具有很好的適用性。在求解零點的過程中，本文提出了一種實部-虛部交叉迭代法，這種求解方法與采用Newton-Raphson迭代求解的方法相比，可以大幅度減少計算量。本文還闡述了在固定諧振模式下電力網絡參與因子矩陣和節點電壓振型的意義。通過算例分析可以發現：

(1)元件的負電阻效應是導致電力網絡諧振不穩定的主要原因，同步發電機在次同步低頻段(低于基頻)由于異步發電機效應會表現出負電阻；

(2)雙饋風機在次同步頻段也存在著負電阻效應，它的接入可能會引起電力網絡在次同步頻段的諧振不穩定；

(3)電力網絡的諧振模式通常表現為2種形式：一種是節點對地的諧振模式，主要特征為節點電壓的振型接近于同一方向；另一種是節點對節點的諧振模式，主要特征為節點電壓的振型會存在2個明顯相反的方向。

[1] Kundur Prabha, Balu Neal J, Lauby Mark G. Power system stability and control[M]. McGraw-hill New York, 1994.

[2] 徐政. 交直流電力系統動態行為分析[M]. 北京: 機械工業出版社, 2005.

[3] 呂敬, 董鵬, 施剛，等. 大型雙饋風電場經MMC-HVDC并網的次同步振蕩及其抑制[J]. 中國電機工程學報, 2015, 35(19): 4852-4860.

Lü Jing, DONG Peng, SHI Gang,et al. Subsynchronous Oscillation and Its Mitigation of MMC-Based HVDC With Large Doubly-Fed Induction Generator-Based Wind Farm Integration[J]. Proceedings of the CSEE, 2015, 35(19): 4852-4860.

[4] 李云豐, 湯廣福, 賀之淵，等. MMC型直流輸電系統阻尼控制策略研究[J]. 中國電機工程學報, 2016,36(20): 5492-5503.

LI Yunfeng, TANG Guangfu, HE Zhiyuan,et al. Damping Control Strategy Research for MMC Based HVDC System[J]. Proceedings of the CSEE, 2016, 36(20): 5492-5503.

[5] 李明節, 于釗, 許濤，等. 新能源并網系統引發的復雜振蕩問題及其對策研究[J]. 電網技術, 2017, 41(4): 1035-1042.

LI Mingjie, YU Zhao, XU Tao,et al. Study of Complex Oscillation Caused by Renewable Energy Integration and Its Solution[J]. Power System Technology, 2017, 41(4): 1035-1042.

[6] 王錫凡, 方萬良, 杜正春. 現代電力系統分析[M]. 北京：科學出版社, 2003.

[7] Semlyen A I. s-domain methodology for assessing the small signal stability of complex systems in nonsinusoidal steady state[J]. IEEE Transactions on Power Systems, 1999, 14(1): 132-137.

[8] Gomes S, Martins N, Portela C. Modal analysis applied to s-domain models of AC networks[C]//IEEE Power and Energy Society.Proceedings of 2001 IEEE Power Engineering Society Winter Meeting,2001.1.28—2001.2.1, Columbus, Ohio, USA: 1305-1310.

[9] Lima L T G, Martins N, Carneiro S. Augmented state-space formulation for the study of electric networks including distributed-parameter transmission line models [C]// IPST commission.Proceedings of 1999 International conference on power systems transients,1999.6.20-24, Budapest, Hungary:87-92.

[10] Varricchio S L, Martins N, Lima L T G. A Newton-Raphson method based on eigenvalue sensitivities to improve harmonic voltage performance[J]. IEEE Transactions on Power Delivery, 2003, 18(1): 334-342.

[11] Xu W, Huang Z, Cui Y, et al. Harmonic resonance mode analysis[J]. IEEE Transactions on Power Delivery, 2005, 20(2): 1182-1190.

[12] IEEE subsynchronous resonance task force. First benchmark model for computer simulation of subsynchronous resonance[J]. IEEE Trans on Power Apparatus and Systems, 1977,96(5):1565-1572.

[13] Xu Zheng, Feng Zhouyan. A novel unified approach for analyzing small-signal stability of power systems[C]//IEEE Power and Energy Society.Proceedings of 2000 IEEE Power Engineering Society Winter Meeting,2000.1.23-27, Singapore:963-967.

[14] Ramos R A, Hiskens I, Canizares C, et al. Benchmark systems for small-signal stability analysis and control[R]. IEEE Power and Energy Society Technical Report PES-TR18, 2015.

[15] Miao Z. Impedance-Model-Based SSR Analysis for Type 3 Wind Generator and Series-Compensated Network[J]. IEEE Transactions on Energy Conversion, 2012, 27(4): 984-991.

2017-09-10

徐政(1962)，男，博士，教授，通信作者，主要研究方向：大規模交直流電力系統分析、直流輸電與柔性交流輸電、風力發電技術與風電場并網技術；

王世佳(1991)，男，博士研究生，主要研究方向：交直流電力系統次同步振蕩；

邢法財(1993)，男，博士研究生，主要研究方向：新能源并網的交直流系統穩定性研究；

肖晃慶(1990)，男，博士研究生，主要研究方向：直流輸電與柔性交流輸電。

(編輯 魏希輝)

QualitativeAnalysisMethodofElectricNetworkResonanceStability

XU Zheng, WANG Shijia, XING Facai, XIAO Huangqing

(College of Electrical Engineering, Zhejiang University, Hangzhou 310027, China)

With the increasing utilization of power electronic equipment in power systems, in recent years,a number of new oscillations with unknown mechanisms have emerged. This paper puts forward the concept of the electric network resonance stability, and tries to classify the unclear reason oscillations mentioned above into the electric network resonance instability category. Thus, many complex power system oscillations can be analyzed by the linear network theory, which provides an approach to solve the actual power system problems by the mathematically mature linear system theory. The objective of this paper is to establish a method for analyzing the electric network resonance stability. By introducing the s-domain nodal admittance matrix of the electric network, this paper transforms the discrimination of the electric network resonance stability into the distribution problem in complex plane of zero point of the determinant of the s-domain nodal admittance matrix. Firstly, it is proved that the zero points of the determinant of the s-domain nodal admittance matrix are actually the eigenvalues of the system. Secondly, we use the cross iteration method of the real part and the imaginary part of the zero point to solve the zero points of the determinant of the s-domain nodal admittance matrix. Thirdly, we derive the nodal voltage mode shape and the participation factor matrix corresponding to a particular resonance mode, which can be used to locate the resonant region of this particular resonance mode in the network. Finally, we illustrate the effectiveness of the proposed method for analyzing the resonance stability of electric networks by several studied cases.

electric network; state space model; s-domain nodal admittance matrix; resonance mode; resonance stability;nodal voltage mode shape; participation factor matrix

國家高技術研究發展計劃項目(863計劃)(2011AA05A119)；國家電網公司科技項目(柔性輸電網規劃評估方法及應用關鍵技術研究)

Project supported by The National High Technology Research and Development of China (863 Program)( 2011AA05A119)

TM711

A

1000-7229(2017)11-0001-08

10.3969/j.issn.1000-7229.2017.11.001