張 哲,喬登科,汪云濤,薛靜瑋,林 毅,薛安成
(1.華北電力大學(xué) 新能源電力系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 102206;2.國(guó)家電網(wǎng)有限公司華北分部,北京 100053;3.國(guó)網(wǎng)福建省電力有限公司經(jīng)濟(jì)技術(shù)研究院,福建 福州 350012)
相較于傳統(tǒng)同步機(jī),風(fēng)機(jī)慣量小、控制結(jié)構(gòu)復(fù)雜、非線性特性明顯,故高比例風(fēng)電并網(wǎng)會(huì)改變電力系統(tǒng)的動(dòng)態(tài)行為并降低系統(tǒng)穩(wěn)定性[1-2],使系統(tǒng)穩(wěn)定機(jī)理更加復(fù)雜[3]。其中直驅(qū)式永磁同步發(fā)電機(jī)(direct-drive permanent magnet synchronous generator,PMSG)并網(wǎng)所引起的次/超同步振蕩現(xiàn)象受到了學(xué)術(shù)界的廣泛關(guān)注[4]。
含風(fēng)電的電力系統(tǒng)振蕩現(xiàn)象按其數(shù)學(xué)表征可分為負(fù)阻尼振蕩、光滑的強(qiáng)迫振蕩、切換型振蕩以及其他復(fù)雜振蕩[5]。基于負(fù)阻尼振蕩和光滑的強(qiáng)迫振蕩的次/超同步分析法應(yīng)用廣泛[6-10],在分析平衡點(diǎn)附近的小范圍振蕩特性方面具有強(qiáng)大優(yōu)勢(shì)。但該方法忽略了系統(tǒng)中限幅等非線性環(huán)節(jié),無法反映系統(tǒng)大范圍的振蕩特性。然而,大范圍切換型振蕩和復(fù)雜振蕩通常因大擾動(dòng)以及多個(gè)非線性環(huán)節(jié)交互作用而產(chǎn)生[11],如部分負(fù)阻尼發(fā)散后觸及限幅引起的振蕩、正阻尼下非線性環(huán)節(jié)參與的振蕩等,其產(chǎn)生機(jī)理復(fù)雜,研究相對(duì)較少。
目前,對(duì)于非線性環(huán)節(jié)參與引起的大范圍切換型振蕩研究已取得了一定的成果[11-19]。對(duì)于風(fēng)機(jī)非線性環(huán)節(jié)參與的大范圍切換型振蕩研究基本可分為如下2 類:基于時(shí)域模型的特性和機(jī)理分析;基于頻域模型的切換型振蕩近似分析。
在基于時(shí)域模型的特性和機(jī)理分析方面,現(xiàn)有研究主要分析如下2 種類型的振蕩。①考慮非線性環(huán)節(jié)導(dǎo)致系統(tǒng)振蕩,如文獻(xiàn)[12]通過仿真分析表明風(fēng)機(jī)控制策略切換可能會(huì)使系統(tǒng)因失去平衡點(diǎn)而失去穩(wěn)定。②具有穩(wěn)定平衡點(diǎn)的系統(tǒng)因限幅飽和引起的切換型振蕩,如:文獻(xiàn)[13]基于時(shí)域仿真法分析正阻尼PMSG 系統(tǒng)在大擾動(dòng)后由限幅引起的次同步振蕩的影響因素和非光滑分岔特性;文獻(xiàn)[14]基于時(shí)域仿真法分析正阻尼雙饋系統(tǒng)在大擾動(dòng)后由限幅引起的次同步振蕩的影響因素和非光滑分岔特性;文獻(xiàn)[15]基于并網(wǎng)電壓源變流器(voltage source converter,VSC)系統(tǒng)8 階狀態(tài)空間模型以及時(shí)域仿真和降階系統(tǒng),分析了電流限幅持續(xù)飽和引起的切換型振蕩機(jī)理及影響振蕩頻率的因素。
在基于頻域模型的切換型振蕩近似分析方面,現(xiàn)有文獻(xiàn)主要采用描述函數(shù)法分析非線性環(huán)節(jié)對(duì)振蕩的影響。文獻(xiàn)[16]提出了基于大信號(hào)阻抗的并網(wǎng)VSC 諧振與抑制方法,利用描述函數(shù)法推導(dǎo)了包含脈寬調(diào)制(pulse width modulation,PWM)飽和的系統(tǒng)響應(yīng)特性。文獻(xiàn)[17]利用描述函數(shù)法,建立并分析了考慮內(nèi)部限幅環(huán)節(jié)的并網(wǎng)VSC非線性傳遞函數(shù)模型。文獻(xiàn)[18]從閉環(huán)反饋原理的角度探討了限幅引起等幅振蕩的產(chǎn)生機(jī)理,提出了考慮限幅非線性特性和dq軸耦合效應(yīng)的風(fēng)電并網(wǎng)系統(tǒng)頻域傳遞函數(shù)模型。文獻(xiàn)[19]基于單邊限幅的描述函數(shù)和廣義Nyquist 判據(jù),近似分析了不對(duì)稱單邊d軸電流限幅(asymmetric unilaterald-axis current limit,AU-d-CL)參與的并網(wǎng)VSC負(fù)阻尼切換型次同步振蕩現(xiàn)象。
值得注意的是,在現(xiàn)有的PMSG 系統(tǒng)切換型振蕩分析研究中,對(duì)正阻尼的切換型振蕩僅有時(shí)域仿真分析[13],未見頻域近似分析。同時(shí),在頻域分析中,僅對(duì)限幅持續(xù)飽和的系統(tǒng)進(jìn)行了分析,而對(duì)大擾動(dòng)后可能會(huì)出現(xiàn)的其他類型限幅飽和情況未見報(bào)道。另一方面,并網(wǎng)VSC 系統(tǒng)可能出現(xiàn)負(fù)阻尼下AU-d-CL 間歇性飽和參與的次/超同步振蕩[19],也可能出現(xiàn)正阻尼下限幅參與的次同步振蕩[15]。對(duì)于PMSG系統(tǒng),是否也可能出現(xiàn)正阻尼下的AU-d-CL間歇性飽和參與的次/超同步振蕩也未見報(bào)道。
鑒于此,本文先展現(xiàn)大擾動(dòng)后AU-d-CL 間歇飽和所引起的PMSG 并網(wǎng)系統(tǒng)切換型振蕩現(xiàn)象,然后結(jié)合AU-d-CL 的描述函數(shù)和廣義Nyquist 判據(jù)分析該振蕩。具體地,首先給出了PMSG 并網(wǎng)系統(tǒng)13 階簡(jiǎn)化狀態(tài)空間模型,并分析其小干擾穩(wěn)定性以及振蕩發(fā)生的參數(shù)條件;其次發(fā)現(xiàn)了一種新型的切換型振蕩現(xiàn)象,即大擾動(dòng)后PMSG 并網(wǎng)系統(tǒng)因AU-d-CL間歇飽和引起的切換型振蕩;再次結(jié)合PMSG 并網(wǎng)系統(tǒng)網(wǎng)側(cè)變流器的頻域模型以及AU-d-CL的描述函數(shù),給出了含限幅環(huán)節(jié)的PMSG 系統(tǒng)近似分析模型;最后結(jié)合描述函數(shù)和廣義Nyquist 判據(jù),近似分析了不同限幅值和參數(shù)下的新型振蕩頻率。
PMSG 并網(wǎng)系統(tǒng)包含風(fēng)輪機(jī)、永磁發(fā)電機(jī)、機(jī)網(wǎng)側(cè)變流器及其相關(guān)控制系統(tǒng)、濾波電感、鎖相環(huán)等部分,主電路如附錄A圖A1所示。
文獻(xiàn)[20]通過聯(lián)立PMSG 機(jī)械轉(zhuǎn)矩表達(dá)式和軸系運(yùn)動(dòng)方程、dq同步旋轉(zhuǎn)坐標(biāo)系下PMSG 的動(dòng)態(tài)方程、全功率變流器機(jī)網(wǎng)側(cè)控制方程、直流電容方程、鎖相環(huán)動(dòng)態(tài)方程以及主回路電流微分方程,推導(dǎo)出PMSG 并網(wǎng)系統(tǒng)簡(jiǎn)化13 階模型,具體見附錄A 式(A1)所示。13 個(gè)狀態(tài)變量分別為Ω、ids、iqs、uds、uqs、φ、ω、ed、eq、id、iq、x、Udc。其中:Ω為發(fā)電機(jī)機(jī)械角速度;φ為鎖相環(huán)輸出量;ω=ωpll-ω0為鎖相環(huán)輸出角頻率ωpll與同步角頻率ω0之差;x為電壓環(huán)中間變量;uds、uqs和ids、iqs分別為定子電壓和電流的d、q軸分量;ed、eq和id、iq分別為網(wǎng)側(cè)變流器電壓和電流的d、q軸分量;Udc為直流母線電壓。
PMSG 并網(wǎng)系統(tǒng)中平衡點(diǎn)穩(wěn)定性由對(duì)應(yīng)雅可比矩陣J的特征值決定,J由附錄A式(A1)在平衡點(diǎn)處線性化得到。J中第i行第j列元素Ji,j的表達(dá)式為:
式中:fi為狀態(tài)空間方程中第i個(gè)方程;xj為第j個(gè)狀態(tài)變量;xequ為平衡點(diǎn)處狀態(tài)變量的值。在附錄B 表B1 所示典型參數(shù)下,系統(tǒng)的平衡點(diǎn)和特征根如附錄B 表B2 所示。系統(tǒng)含有EQ1—EQ4這4 個(gè)平衡點(diǎn):EQ1為穩(wěn)定運(yùn)行的平衡點(diǎn);EQ2為含有2 個(gè)正實(shí)部特征值的不穩(wěn)定鞍焦點(diǎn);EQ3為含有1對(duì)正實(shí)部共軛特征值的不穩(wěn)定結(jié)焦點(diǎn);EQ4為含有1個(gè)正實(shí)部特征值的不穩(wěn)定鞍焦點(diǎn)。由其對(duì)應(yīng)的狀態(tài)變量(φ、id)的值可知:EQ1為dPout/did>0 時(shí)的穩(wěn)定平衡點(diǎn),Pout為風(fēng)機(jī)輸出功率;EQ2為dPout/did<0 時(shí)的不穩(wěn)定平衡點(diǎn);EQ3和EQ4則是因風(fēng)機(jī)輸出功率過大而導(dǎo)致的不穩(wěn)定靜態(tài)工作點(diǎn)[21]。
對(duì)穩(wěn)定運(yùn)行的平衡點(diǎn)EQ1,含有多種振蕩模式,各振蕩模式的頻率、阻尼比和主要參與變量如附錄B 表B3所示。其中振蕩模式1主要由機(jī)側(cè)變流器相關(guān)控制變量ids、iqs、uds、uqs參與,可稱為機(jī)側(cè)變流器內(nèi)環(huán)振蕩模式;振蕩模式2 由機(jī)側(cè)變流器相關(guān)控制變量Ω、ids、iqs、uds、uqs參與,且主要與Ω相關(guān),故可稱為軸系振蕩模式;振蕩模式3、4 主要由鎖相環(huán)相關(guān)控制變量ω、φ和eq、iq參與,主要與ω、φ相關(guān),可視為鎖相環(huán)振蕩模式;振蕩模式5主要由網(wǎng)側(cè)變流器d軸相關(guān)控制變量ed、id和x參與,主要與ed、id相關(guān),視為d軸振蕩模式。
進(jìn)一步,分析各振蕩模式隨單一參數(shù)變化的變化規(guī)律(改變某一參數(shù)時(shí)維持其他參數(shù)為初始值不變)。參數(shù)對(duì)振蕩模式3的穩(wěn)定性影響結(jié)果如附錄B表B4所示,表中數(shù)值代表振蕩模式變?yōu)椴环€(wěn)定的參數(shù)臨界值。由于振蕩模式3受諸多參數(shù)的影響,將其設(shè)為PMSG 并網(wǎng)系統(tǒng)中參數(shù)改變所引起的主要振蕩模式。而其他振蕩模式受參數(shù)變化的影響較不明顯,其余振蕩模式隨單一參數(shù)變化的變化規(guī)律未列出。
運(yùn)行于穩(wěn)定平衡點(diǎn)的正阻尼系統(tǒng)在受到小擾動(dòng)后系統(tǒng)軌線會(huì)逐漸收斂至平衡點(diǎn);在系統(tǒng)受到大擾動(dòng)后,可通過電流內(nèi)環(huán)設(shè)置的限幅環(huán)節(jié)限制設(shè)備過流過壓。電流達(dá)到限幅值后控制環(huán)節(jié)進(jìn)行動(dòng)態(tài)切換,形成系統(tǒng)的切換型振蕩現(xiàn)象[5]。
設(shè)置PMSG 并網(wǎng)系統(tǒng)d軸電流參考值idref限幅上限值idmax為1 200 A,仿真分析系統(tǒng)在遭受大擾動(dòng)后的動(dòng)態(tài)行為。有/無限幅時(shí)網(wǎng)側(cè)變流器d軸電流id的響應(yīng)如圖1 所示。由圖可知:無限幅環(huán)節(jié)時(shí),系統(tǒng)在受到大擾動(dòng)后恢復(fù)穩(wěn)定(d軸電流穩(wěn)定值id0=1 036.04 A),表明系統(tǒng)為一光滑正阻尼系統(tǒng),與特征值分析結(jié)果相同;加入限幅后,系統(tǒng)在穩(wěn)定平衡點(diǎn)附近為一光滑系統(tǒng),全局上為一非光滑系統(tǒng)。
圖1 有/無限幅環(huán)節(jié)時(shí)id響應(yīng)Fig.1 Response of id with/without limit link
在相同擾動(dòng)下,系統(tǒng)出現(xiàn)振蕩,d軸電流頻譜分析如附錄C 圖C1 所示,振蕩頻率為136.7 Hz。idref和id波形如附錄C 圖C2 所示,大擾動(dòng)后id波形在限幅上限值與振蕩下限值間不斷變化,限幅環(huán)節(jié)間歇飽和。id-idref相圖如圖2 所示,系統(tǒng)軌線形成一不光滑極限環(huán),當(dāng)限幅環(huán)節(jié)飽和時(shí),id-idref軌線被限制在idmax=1 200 A 的邊界面上并在邊界面上滑動(dòng);當(dāng)限幅環(huán)節(jié)未飽和時(shí),id-idref軌線運(yùn)行未受限制,在平面上滑動(dòng),在兩軌線相交處出現(xiàn)不光滑現(xiàn)象(軌線切換),故該振蕩為一切換型振蕩現(xiàn)象。
圖2 id -idref相圖Fig.2 Phase diagram of id -idref
進(jìn)一步,分析改變限幅上限對(duì)系統(tǒng)的振蕩特性的影響。當(dāng)d軸電流參考值限幅上限值idmax分別為1 150、1 175、1 200 A 時(shí),大擾動(dòng)下振蕩頻率分別為139.5、138.9、138.21 Hz。振蕩頻率隨idmax的改變而變動(dòng),與文獻(xiàn)[19]中改變限幅上限值不影響振蕩頻率的結(jié)論有明顯區(qū)別。
本節(jié)分析正阻尼PMSG 并網(wǎng)系統(tǒng)大擾動(dòng)后發(fā)生限幅間歇飽和型切換型振蕩參數(shù)條件。分析過程中,當(dāng)改變單一參數(shù)時(shí),保持其他參數(shù)為附錄B 表B1 所示初始值不變。由仿真分析可得,原PMSG 并網(wǎng)系統(tǒng)因大擾動(dòng)發(fā)生限幅間歇性飽和的切換型振蕩的參數(shù)范圍如表1所示。
表1 振蕩發(fā)生的參數(shù)范圍Table 1 Parameter ranges of oscillation occurrence
無限幅時(shí)系統(tǒng)在大擾動(dòng)后恢復(fù)穩(wěn)定,有限幅時(shí)系統(tǒng)出現(xiàn)了持續(xù)性振蕩現(xiàn)象,可知振蕩由d軸電流參考值限幅間歇飽和引起。振蕩期間因限幅環(huán)節(jié)間歇飽和,直流電壓環(huán)節(jié)不能忽略,不能對(duì)系統(tǒng)進(jìn)行數(shù)學(xué)及物理降階分析[15]。有鑒于此,需探尋限幅間歇飽和的切換型振蕩分析方法。在此引入描述函數(shù)法。
描述函數(shù)法是一種廣泛應(yīng)用于非線性系統(tǒng)的振蕩和穩(wěn)定性分析的頻域近似分析方法[22]。對(duì)于一類含非線性環(huán)節(jié)的系統(tǒng),其在頻域內(nèi)可近似為圖3 所示的單位負(fù)反饋系統(tǒng)。圖中:r(t)、c(t)分別為輸入、輸出信號(hào);N(A)為非線性部分的描述函數(shù),其表達(dá)式見式(2)。
圖3 頻域近似系統(tǒng)Fig.3 Approximation system in frequency domain
式中:A、θ分別為輸出信號(hào)的基波幅值和相角;A0為輸入正弦波幅值。
進(jìn)而可根據(jù)廣義Nyquist 判據(jù)和式(3)分析系統(tǒng)穩(wěn)定性并計(jì)算發(fā)生穩(wěn)定振蕩的振蕩頻率。
基于描述函數(shù)和Nyquist 廣義判據(jù)的振蕩分析方法的前提是得到線性部分和非線性部分的頻域模型。因此,下文推導(dǎo)PMSG 并網(wǎng)系統(tǒng)的單位負(fù)反饋線性系統(tǒng)形式的頻域近似模型和AU-d-CL 的描述函數(shù)。
為簡(jiǎn)化分析過程,系統(tǒng)參數(shù)使用標(biāo)幺值,并假定機(jī)側(cè)變流器輸出功率恒定,可將系統(tǒng)簡(jiǎn)化為8 階VSC 并網(wǎng)系統(tǒng)。經(jīng)推導(dǎo)后可得系統(tǒng)以直流電壓參考值擾動(dòng)ΔUdcref為輸入,直流電壓擾動(dòng)ΔUdc為輸出的閉環(huán)傳遞函數(shù),如式(4)所示。
式中:Hv(s)為電壓環(huán)的比例積分控制器的傳遞函數(shù),Hv(s)=Kp+Ki/s,Kp、Ki分別為比例、積分系數(shù);此處N(A)表示限幅環(huán)節(jié)描述函數(shù);G1(s)的含義與推導(dǎo)過程可參閱文獻(xiàn)[19]。
分離式(4)可得圖3 所示的G(s)部分。G(s)傳遞函數(shù)為:
對(duì)于本文所考慮的AU-d-CL 間歇飽和現(xiàn)象,由頻譜分析可知輸入量idref具有一定的直流分量,輸入幅值增大時(shí)先觸及限幅上限,從而形成不對(duì)稱的單邊限幅環(huán)節(jié)飽和輸出,且振蕩中心(id=1 067 A)未與穩(wěn)態(tài)運(yùn)行值(id0=1 036.04 A)重合??紤]上述因素后可得不對(duì)稱單邊限幅的描述函數(shù)為:
式中:a為限幅上限值;A1為穩(wěn)態(tài)偏差值。
圖4 為G(s)與-1/N(A)的特性曲線。G(s)為角頻率從0到+∞變化的Nyquist曲線。-1/N(A)為負(fù)倒特性曲線/臨界線,表示輸入信號(hào)幅值A(chǔ)從0 到+∞變化的1 條曲線,由于考慮單邊限幅飽和中的輸入側(cè)非周期量的影響,-1/N(A)曲線將隨著A的變化在復(fù)平面上運(yùn)動(dòng),而不僅限于負(fù)半軸上。
圖4 G(s)與-1/N(A)曲線Fig.4 Curves of G(s)and -1/N(A)
在保持系統(tǒng)參數(shù)不變的情況下,可在復(fù)平面上畫出G(s)的Nyquist 曲線及-1/N(A)的負(fù)倒特性曲線。由圖4 可知:隨著-1/N(A)中幅值A(chǔ)的不斷增大,負(fù)倒特性曲線(黑線)與Nyquist 曲線(藍(lán)線)存在多個(gè)交點(diǎn);且負(fù)倒特性曲線穿越至Nyquist 曲線不包圍區(qū)域,這說明系統(tǒng)在該臨界點(diǎn)處存在穩(wěn)定的等幅振蕩。
進(jìn)一步計(jì)算交點(diǎn)處的振蕩頻率fs,即等同于通過Nyquist 曲線得到的交點(diǎn)處的振蕩角頻率ωs。當(dāng)Kiv=50 時(shí),由理論計(jì)算得到的系統(tǒng)振蕩頻率fs為134.85 Hz,與基于仿真頻譜分析獲得的系統(tǒng)振蕩頻率136.7 Hz較為接近。
不同Kiv下的Nyquist 曲線與-1/N(A)曲線的相交示意圖如圖5 所示。理論計(jì)算得到的系統(tǒng)振蕩頻率、基于仿真頻譜分析得到的系統(tǒng)振蕩頻率以及兩者之間的誤差如表2所示。
表2 不同Kiv下的仿真與計(jì)算結(jié)果對(duì)比(idmax=1 200 A)Table 2 Comparison of simulative and calculative results with different Kiv values(idmax=1 200 A)
圖5 不同Kiv下G(s)與-1/N(A)曲線Fig.5 Curves of G(s)and -1/N(A)with different Kiv values
進(jìn)一步,在不改變系統(tǒng)動(dòng)態(tài)行為(振蕩期間限幅間歇性飽和)的前提下,改變限幅上限值,經(jīng)過限幅環(huán)節(jié)前、后idref波形圖如圖6 所示。不同限幅值下振蕩頻率的仿真與計(jì)算結(jié)果對(duì)比變化如表3所示。
表3 不同idmax下的仿真與計(jì)算結(jié)果對(duì)比(Kiv=50)Table 3 Comparison of simulative and calculative results with different idmax(Kiv=50)
圖6 idmax=1 100 A與idmax=1 150 A時(shí)振蕩波形Fig.6 Oscillation waveforms with idmax=1 100 A and idmax=1 150 A
表2、3 表明,在不同控制參數(shù)或限幅上限下,振蕩頻率計(jì)算結(jié)果與仿真結(jié)果之間存在1 Hz 左右的誤差。原因在于分析所用頻域模型為PMSG 網(wǎng)側(cè)變流器頻域模型,忽略了機(jī)側(cè)變流器以及風(fēng)機(jī)、發(fā)電機(jī)等部分;idref振蕩波形并非完美的正弦波形,應(yīng)用描述函數(shù)法時(shí)忽略了存在的微小的諧波分量。以上均會(huì)導(dǎo)致計(jì)算結(jié)果存在誤差。
進(jìn)一步,由表2 可知:不同Kiv下振蕩頻率不同,這是由于改變Kiv使得Nyquist 曲線的位置發(fā)生改變(如圖5 所示),與-1/N(A)曲線的交點(diǎn)在復(fù)平面上不斷變動(dòng)使振蕩頻率發(fā)生變化;不同idmax下振蕩頻率不同,這是因?yàn)閕dmax值的變化改變了不對(duì)稱限幅的描述函數(shù)(限幅上限值a),使-1/N(A)曲線的位置發(fā)生改變,與Nyquist 曲線的交點(diǎn)在復(fù)平面上不斷變動(dòng)使振蕩頻率發(fā)生變化。而在限幅參與的負(fù)阻尼振蕩中,改變限幅上限值將不影響振蕩頻率,這與本文的分析結(jié)果存在本質(zhì)區(qū)別。
本文基于簡(jiǎn)化的PMSG 并網(wǎng)系統(tǒng)狀態(tài)空間模型,發(fā)現(xiàn)了大擾動(dòng)后PMSG 并網(wǎng)系統(tǒng)AU-d-CL 間歇飽和引起的次/超同步振蕩現(xiàn)象,并結(jié)合描述函數(shù)法和廣義Nyquist 判據(jù)對(duì)振蕩頻率進(jìn)行了近似分析。結(jié)論如下:
1)發(fā)現(xiàn)了一種新型的切換型振蕩現(xiàn)象,即大擾動(dòng)后的正阻尼PMSG 并網(wǎng)系統(tǒng)因AU-d-CL 間歇飽和所引起的切換型次/超同步振蕩,并分析了振蕩發(fā)生的參數(shù)條件;
2)結(jié)合AU-d-CL 的描述函數(shù)和廣義Nyquist 判據(jù),近似分析了不同限幅上限值和參數(shù)下的振蕩頻率;
3)解釋了由AU-d-CL 間歇飽和引起的切換型振蕩頻率會(huì)隨限幅上限降低而增加的現(xiàn)象。
特別地,描述函數(shù)法適用前提是系統(tǒng)的線性部分具有低通濾波特性。因此,當(dāng)系統(tǒng)的低通特性改變時(shí),近似準(zhǔn)確度可能發(fā)生變化。
附錄見本刊網(wǎng)絡(luò)版(http://www.epae.cn)。