一種頻率信號的自適應測量法

劉亞棟 楊洪耕 馬 超 劉 張

（四川大學電氣信息學院 成都 610065）

1 引言

自從20世紀早期交流電廣泛使用開始，電力系統(tǒng)諧波問題就引起了人們的關(guān)注[1]。近年來，一方面，電子開關(guān)型、電弧型等非線性裝置的廣泛使用給電力系統(tǒng)帶來了大量的諧波、間諧波污染；另一方面，大量敏感設備的使用對于諧波、間諧波的測量與補償提出了更高的要求[2]。因此，電力系統(tǒng)諧波、間諧波等頻率信號的治理是當前一個十分重要的研究領(lǐng)域，而精確地測量出諧波、間諧波是對其進行治理的前提[3,4]。

常規(guī)的頻率信號測量采用時域-頻域的方法[3-8]。這種方法的誤差由采樣是否同步，采樣窗函數(shù)性能以及頻率分辨率等因素共同決定[8]。隨著電力系統(tǒng)中非線性、沖擊性負荷的增多，電網(wǎng)頻率不能維持50Hz的額定值，而是以其為中心，在較小范圍內(nèi)波動。對于此環(huán)境中非穩(wěn)定諧波、間諧波的測量而言，一方面，需要盡可能高的頻率分辨率，另一方面，采樣窗口長度要盡量短[8,9]。為此，IEC61000-4-7標準對這種基于頻譜分析的頻率信號測量方法做出了規(guī)范：推薦50Hz系統(tǒng)采樣窗口長度為10個基頻周波，在此頻譜范圍內(nèi)通過諧波、間諧波分組來測量信號各個頻率分量。標準是通過對頻率分辨率與采樣窗口長度兩方面要求綜合考慮之后才確定出該折中方案[3]，基于此頻率分辨率的各種方法在非穩(wěn)定諧波、間諧波測量中很難得到理想精度的結(jié)果[8-11]；標準推薦的利用間諧波組測量間諧波的方法是將兩個相鄰諧波之間的頻譜作為一個間諧波組，從中得到這兩個諧波之間間諧波的總體有效值，而無法確定出每一個間諧波分量的準確信息[12,13]。

為此，本文提出了一種能夠有效測量電力信號中各頻率分量的自適應算法。依據(jù)梯度下降法構(gòu)建出頻率信號參數(shù)的自適應計算模型，采用李亞普諾夫穩(wěn)定性定理分析論證了該模型在待測信號非穩(wěn)定條件下的漸近穩(wěn)定性，對模型中的參數(shù)因子按照動態(tài)調(diào)節(jié)策略進行控制以兼顧自適應測量算法靈敏性和穩(wěn)定性的要求。

2 頻率信號的自適應測量法

2.1 自適應測量算法模型

正弦信號的自適應測量算法就是用預測所得正弦信號去擬合待測信號。在此過程中，根據(jù)前面的趨勢及誤差來預測當前值，再將當前預測值與實測值進行比較以確定新的趨勢與誤差，直到兩者相擬合。這種自適應測量算法能夠?qū)Ψ?、頻率及相位隨時間變化的非穩(wěn)定正弦信號進行跟蹤測量。

設非穩(wěn)定信號 u(t)中任一待測正弦頻率分量y( t) =A( t) sin （∫ω (τ)dτ + δ(t )），其中 A(t)，ω(t)及δ(t)分別為隨時間變化的幅值、角頻率以及初相位等信號參數(shù)，總相位 φ(t )=∫ω (τ)dτ + δ(t )。要從信號u(t)中檢測出該正弦頻率分量，就要準確測量出該分量的信號參數(shù)，以保證測量所得的輸入與輸出的差值信號盡可能的小。為此，定義待測頻率分量的幅值預測值、角頻率預測值ω?(t)以及總相位預測值，則可得待測頻率分量的預測值為

進一步可知測量所得輸入與輸出的差值信號為

定義目標函數(shù)為

式中，mi(i=1,2,3)為正數(shù)。依次定義自適應測量算法的幅值因子Ka、角頻率因子Kf及相位因子Kp為：Ka=2m1，Kf=2m2，Kp=m3/m2，則可得如下3階動態(tài)系統(tǒng)所示的自適應測量算法模型，其狀態(tài)變量分別為幅值預測值、角頻率預測值ω?(t)及總相位預測值。

在具體軟件計算中，對式（2）及式（8）～式（10）進行離散化處理，則有

在式（8）～式（10）所示的自適應測量算法模型中，參數(shù)因子Ka, Kf及Kp用來控制算法的收斂速度與穩(wěn)態(tài)誤差。只要為該算法設置角頻率初始值ω0以及上述 3個參數(shù)因子，就能保證預測信號與輸入信號中角頻率最接近于ω0的頻率分量相擬合。

2.2 自適應測量算法的穩(wěn)定收斂性

在 2.1節(jié)中給出的頻率信號自適應測量算法模型是一個非線性系統(tǒng)，其穩(wěn)定收斂性往往通過該系統(tǒng)的平衡點或周期解的穩(wěn)定性來分析。在此，依據(jù)李亞普諾夫穩(wěn)定性定理對該系統(tǒng)平衡點的穩(wěn)定性進行分析。

假定信號u(t)為有界連續(xù)的周期信號，根據(jù)傅里葉級數(shù)，同時考慮到信號參數(shù)可能隨時間在較小范圍內(nèi)變化，可以將其寫成

式中，N為u(t)中所含頻率分量的個數(shù)；Ai(t)、ωi(t)及 δi(t)分別為隨時間變化的第 i個頻率分量的幅值、角頻率及初相。為便于分析，依據(jù)文獻[14]中的平均理論對u(t)進行線性平滑近似后可得

T為信號周期。則由式（2）及式（17）可得

對式（18）左右兩端關(guān)于時間求微分，并將式（8）～式（10）代入可得

以測量u(t)中第k個頻率分量為例，則其他頻率分量之和即為輸入與輸出的差值信號 e(t)。在平衡點附近近似有：據(jù)此對式（19）進行簡化后根據(jù)文獻[14]中的平均理論可得

2.3 電力頻率信號的自適應測量法

考慮到電力頻率信號中除含有基波，通常還包含少量的諧波（整次或非整次）及噪聲信號。要想分析電力信號中的各個頻率分量，就需要將多個根據(jù) 2.1節(jié)中算法構(gòu)成的自適應測量單元按照一定的結(jié)構(gòu)組合起來，同時根據(jù)對輸入信號進行FFT預處理所得的頻譜范圍針對性地設置每個測量單元的角頻率初始值ω0，以保證對于輸入信號中任一頻率分量都有唯一的自適應測量單位對其進行快速準確的擬合。圖1給出了自適應測量單元的兩種組合結(jié)構(gòu)。

圖1 自適應測量結(jié)構(gòu)Fig.1 Block diagram of adaptive measurement

圖1所示的兩種結(jié)構(gòu)均能實現(xiàn)對輸入信號中各個頻率分量的測量。在串聯(lián)型結(jié)構(gòu)中，每個單元將上一個單元輸入與輸出的差值信號作為輸入，考慮到電力頻率信號中次數(shù)越高的諧波，含量越小，這種結(jié)構(gòu)可以逐次剔除檢測出的低次諧波分量，準確地測量出含量較小的頻率分量。并聯(lián)型結(jié)構(gòu)中的各個單元同時對輸入信號中相應的頻率分量進行自適應測量，具有較高的效率。在實際應用中可根據(jù)需要靈活選用這兩種結(jié)構(gòu)。本文仿真實驗中采用串聯(lián)型結(jié)構(gòu)對算法進行驗證。

3 自適應測量法的動態(tài)參數(shù)因子調(diào)節(jié)策略

3.1 自適應測量算法的參數(shù)因子特性

注意到 2.2節(jié)中關(guān)于自適應測量算法穩(wěn)定收斂特性的證明是一個局部證明，必須合理選取參數(shù)因子 Ka, Kf與 Kp以及角頻率初始值 ω0才能保證算法準確快速收斂于輸入信號中的待測頻率分量。由式（8）～式（10）可知，算法中 Ka主要控制信號幅值的收斂速度及穩(wěn)態(tài)計算誤差，Kf及Kp主要控制頻率和相位的相應特性。進一步由式可得，上述參數(shù)因子取值較大時算法收斂速度較快，但是當信號畸變嚴重時，即使預測信號逼近待測頻率分量，測量所得輸入與輸出的差值信號 e(t)仍較大，此時，較大的參數(shù)因子將進一步增大信號參數(shù)隨時間變化率的計算值，從而在測量結(jié)果中產(chǎn)生較大的穩(wěn)態(tài)誤差；參數(shù)因子取值較小時，各測量單元能夠精確檢測出待測頻率分量，但不能及時跟蹤信號參數(shù)的快速波動，動態(tài)響應速度較慢。

3.2 動態(tài)參數(shù)因子調(diào)節(jié)策略

電力頻率信號的測量算法在理想情況下應該滿足以下要求，當待測信號的幅值、頻率或相位發(fā)生變化時能迅速做出反應，并快速收斂到新的穩(wěn)定狀態(tài)，具有較強的靈敏性；當待測信號比較穩(wěn)定時能有效抑制干擾，準確收斂到待測信號，表現(xiàn)出較好的穩(wěn)定性。

根據(jù) 3.1節(jié)中對自適應測量算法參數(shù)因子特性的分析，要滿足上述測量要求，常規(guī)的思路是綜合算法靈敏性和穩(wěn)定性兩方面要求，折中地設置參數(shù)因子值[15]。目前看來，這種做法并不是理想的解決方案，采用時變參數(shù)因子能夠取得更理想的測量精度和速度[16]。因此，在本算法中，將自適應計算過程中提供的某些回饋量作為衡量標準來調(diào)節(jié)參數(shù)因子，在計算初始階段或者信號參數(shù)跳變時，參數(shù)因子取值相對大一些，以便有較快的收斂速度；計算值逼近待測值時，參數(shù)因子取值相應地減小一些，以保證算法準確收斂于待測頻率信號。

由于輸入信號畸變嚴重時，每個測量單元輸入與輸出的差值信號 e(t)并不能有效反映信號分量預測值與待測值逼近的真實情況，而本算法中以各信號參數(shù)測量誤差的變化趨勢為依據(jù)，動態(tài)調(diào)節(jié)相應的參數(shù)因子。由數(shù)學知識可知，任意函數(shù)的一階導數(shù)能夠反應該函數(shù)的變化趨勢，導數(shù)值越大說明函數(shù)值變化越劇烈。在離散計算情況下，用差分代替微分，分別以幅值、頻率及相位測量誤差的差分值作為回饋值，對相應的參數(shù)因子進行調(diào)節(jié)。由式（10）中頻率與相位的相關(guān)性可知，參數(shù)因子 Kf與 Kp均與相位正相關(guān)，單獨調(diào)節(jié) Kf或者 Kp對相位的計算均有效。因此，為簡化計算過程，本算法中只對幅值因子Ka及頻率因子Kf進行動態(tài)調(diào)節(jié)。

從算法的收斂速度及穩(wěn)態(tài)誤差兩方面考慮，在計算過程中，對參數(shù)因子作如下限制

式中，i=a, f，0＜Kimin＜Kimax。選擇 Kimax時主要考慮信號畸變較嚴重時算法的穩(wěn)態(tài)誤差，Kimin的選擇需要綜合算法的收斂速度及穩(wěn)態(tài)誤差后取一個折中值。具體的參數(shù)因子動態(tài)調(diào)節(jié)策略如下式所示。

式中，A0、ω0分別為幅值與角頻率的待測值；α、β為正常數(shù)。在計算過程中，分別根據(jù)幅值測量誤差與角頻率測量誤差的變化趨勢調(diào)節(jié)相應的參數(shù)因子。以幅值因子為例，在計算初始時刻或者信號幅值跳變時刻，EAk較大，幅值因子取值接近 Kamax，算法具有較快的幅值收斂速度；反之，當幅值測量值逼近待測值時，幅值因子取值接近Kamin，從而有效防止了穩(wěn)態(tài)失調(diào)。

4 仿真分析

在 Matlab7.1仿真環(huán)境下，首先取多個畸變較嚴重的輸入信號，對每一個輸入信號，分別為幅值因子 Ka、頻率因子 Kf以及相位因子 Kp取多組不同的值，在不同的固定參數(shù)因子條件下進行多次實驗。從中選取幾次典型的信號幅值測量結(jié)果，通過對比分析，從這幾次實驗所設置的幅值因子值確定出本文仿真實驗中幅值因子的取值范圍 Kamax及 Kamin。同樣地，根據(jù)多次頻率測量實驗結(jié)果確定出本文仿真實驗中的Kfmax以及Kfmin。確定了Ka與Kf之后繼續(xù)在不同的相位因子Kp條件下進行多次實驗，從中選擇幾次典型測量結(jié)果，進而確定出Kp的一個合理取值范圍，選取該范圍的中間值作為本文仿真實驗中Kp的固定取值。確定了Ka、Kf的范圍以及Kp的值之后，設置不同的正常數(shù)α、β，根據(jù)式（24）和式（25）在信號幅值及頻率跳變情況下對上述幾個畸變較嚴重的輸入信號分別進行多次實驗，通過對比分析，從中選取能夠在盡可能短的動態(tài)過程后準確收斂到信號參數(shù)跳變后新狀態(tài)的正常數(shù)作為本文仿真實驗中α，β 的值。由于本文仿真算例中的參數(shù)因子范圍或取值是預先通過對不同的信號進行大量實驗對比分析確定的，因而這些取值具有一定的普適性。

下面在輸入信號幅值跳變、基波頻率跳變、基頻頻率連續(xù)波動以及其中存在間諧波分量等幾種典型情況下對本文方法進行仿真驗證。為了凸顯本文方法在穩(wěn)定性與靈敏性方面的優(yōu)勢，還將在信號幅值跳變以及基波頻率跳變兩種信號參數(shù)快速波動場合下對采用動態(tài)參數(shù)因子調(diào)節(jié)策略的本文方法及常規(guī)的折中設置固定參數(shù)因子的方法進行實驗對比。在這里，折中設置的Ka與Kf分別取為Kamax、Kamin及 Kfmax、Kfmin的均值，兩種方法中相位因子 Kp的取值相同。

4.1 幅值跳變的諧波信號分析

設仿真所用的負載電流頻率信號i(t)為

式中，基頻額定頻率 ?0=50Hz；δ1、δ2及 δ3為任意常數(shù)；采樣頻率?s=10kHz。

算法中相關(guān)參數(shù)的取值如下：Kamax=200，Kamin=100，Kfmax=9 000，Kfmin= 7 000，Kp=0.15，α=1.5，β=3.5。假設在t=0.2s時，信號i(t)中各個頻率分量的幅值變?yōu)樵档?.8倍，在t=0.4s時又恢復原值。根據(jù)FFT預處理所得的頻譜范圍設置各個自適應測量單元的ω0值依次為 100π、500π、1 400π及 2 200π。限于篇幅，在此只從仿真結(jié)果中選擇具有代表性的波形及幅值/頻率測量曲線來說明本文方法的有效性。圖2所示為信號i(t)及從中檢測出的基頻分量，圖3為測量到的7次諧波幅值及相應的幅值因子曲線。由圖可見，無論在信號幅值穩(wěn)定還是非穩(wěn)定條件下，本文方法均能準確跟蹤信號幅值，快速收斂到待測頻率分量；通過對比分析圖3a中采用動態(tài)參數(shù)因子調(diào)節(jié)策略的本文方法與固定參數(shù)因子方法的測量曲線，可以發(fā)現(xiàn)，本文方法在信號幅值跳變后1個基波周期內(nèi)就能準確收斂到新的穩(wěn)定值，動態(tài)過渡過程更短，具有更高的測量速度，而且不存在明顯的穩(wěn)態(tài)測量誤差，與理論分析的結(jié)果一致。

圖2 幅值跳變情況下自適應測量法的仿真結(jié)果Fig.2 Simulation result of adaptive measuring method when the amplitudes undergo step change

圖3 7次諧波分量的幅值及幅值因子曲線Fig.3 Amplitude and its factor curves of the 7th harmonic component

4.2 基頻頻率跳變的諧波信號分析

負載電流頻率信號與 4.1節(jié)中的 i(t)相同，采樣頻率?s=10kHz。算法中相關(guān)參數(shù)的選取以及各個自適應測量單元的ω0也與上例一致。假設在 t=0.2s時信號i(t)的基頻頻率變?yōu)?0.5Hz，在t=0.4s時頻率又恢復到額定值。圖 4中給出了 i(t)及從中檢測出的基頻分量波形，圖5a中給出了用本文方法與折中設置固定參數(shù)因子的方法檢測到的5次諧波幅值曲線，圖5b為本文方法中相應的幅值因子曲線，圖6a給出了兩種方法測量得到的基頻頻率，圖 6b為本文方法中相應的頻率因子曲線。結(jié)合圖 5及圖6可見：在基頻頻率跳變情形下，本文方法仍能快速準確收斂到待測的各頻率分量，跟蹤測量各分量的幅值及頻率隨時間的變化；相對于折中選取固定參數(shù)因子的方法，動態(tài)參數(shù)因子調(diào)節(jié)策略保證了本文方法在信號頻率跳變后更短過渡過程內(nèi)就能準確收斂到信號新的穩(wěn)定值，而且未出現(xiàn)明顯的穩(wěn)態(tài)誤差，進一步證實了理論分析。

圖4 頻率跳變情況下自適應測量法的仿真結(jié)果Fig.4 Simulation results of adaptive measuring method when the fundamental frequency undergoes a step change

圖5 5次諧波幅值及幅值因子曲線Fig.5 Amplitude and its factor curves of 5th harmonic component

圖6 基頻頻率及頻率因子曲線Fig.6 Frequency and its factor curves of fundamental component

4.3 基頻頻率隨時間波動的跟蹤分析

考慮到實際電力系統(tǒng)中，電網(wǎng)基波頻率通常在較小的范圍內(nèi)連續(xù)波動，為此，進一步對本文方法跟蹤測量頻率信號的能力進行檢驗。假設 4.1節(jié)中負載電流信號i(t)的基頻頻率f0=50+0.5sin(2πt)，即基頻頻率在0.5Hz的范圍內(nèi)正弦波動。則對于本文方法，依次設置各個自適應測量單元的角頻率初始值ω0以及算法中的相關(guān)參數(shù)與4.2節(jié)中相同，對輸入信號中各個頻率分量進行跟蹤測量。圖7中給出了基頻頻率跟蹤曲線。對比分析圖中的兩條頻率曲線，不難發(fā)現(xiàn)，當基波頻率在較小范圍內(nèi)連續(xù)波動時，本文方法能夠在較短的動態(tài)過程之后實現(xiàn)對信號頻率的準確跟蹤測量，具有較理想的動態(tài)測量速度及穩(wěn)態(tài)測量精度。

圖7 隨時間波動的基頻頻率Fig.7 Fluctuations of fundamental frequency with time

4.4 間諧波信號分量分析

假設輸入頻率信號中除了 4.1節(jié)中 i(t)所包含的各次諧波分量，還包含頻率為68Hz及83Hz的間諧波分量，其幅值分別為0.1(pu)及0.08(pu)。此外，為檢驗本文方法在噪聲環(huán)境下的性能，向輸入信號中加入高斯白噪聲以得到信噪比（Signal Noise Ratio，SNR）為20dB的信號。通過對輸入信號取10個周波進行FFT預處理后確定各頻率分量的大致頻譜范圍，據(jù)此設置的對應于兩個間諧波分量的自適應測量單元的角頻率初始值ω0分別為140π及160π。圖7和圖8中分別給出了測量所得的兩個間諧波分量的幅值及頻率曲線。為進一步提高測量精度，在這里，Ka及Kf的初始值較4.1節(jié)中均相對減小，故收斂速度相對較慢。由圖可見，本文方法在信號畸變較嚴重的情況下能夠有效抑制其他頻率分量及噪聲的影響，準確收斂到待測頻率分量，與理論分析結(jié)果一致。

圖8 頻率為68Hz間諧波的幅值與頻率Fig.8 Amplitude and frequency curves of inter-harmonic with 68Hz for frequency

圖9 頻率為83Hz間諧波的幅值與頻率Fig.9 Amplitude and frequency curves of inter-harmonic with 83Hz for frequency

5 結(jié)論

（1）提出的自適應測量法能夠準確測量出輸入頻率信號中各個頻率的諧波、間諧波分量。

（2）在頻率信號穩(wěn)定和非穩(wěn)定條件下，李亞普諾夫穩(wěn)定性定理均能保證該自適應測量法漸近穩(wěn)定地收斂到輸入信號中各個待測的頻率分量。

（3）動態(tài)參數(shù)因子調(diào)節(jié)策略的采用能夠兼顧算法的動態(tài)響應速度和穩(wěn)態(tài)測量誤差兩方面要求，保證算法的靈敏性和穩(wěn)定性。

[1]Lin H C.Inter-harmonic identification using group-harmonic weighting approach based on the FFT[J].IEEE Transaction on Power Electronics, 2008,23(3): 1309-1319.

[2]Clark S L, Famouri P, Cooley W L.Elimination of supply harmonics[J].IEEE Industry Applications Magazine, 1997, 3(2): 62-67.

[3]Apuzzo M D, Arco M D.A time-domain approach for the analysis of nonstationary signals in power systems[J].IEEE Transaction on Power Delivery,2008, 57(9): 1969-1977.

[4]蔡濤, 段善旭, 劉方瑞.基于實值MUSIC算法的電力諧波分析方法[J].電工技術(shù)學報, 2009, 24(12):149-155.Cai Tao, Duan Shanxu, Liu Fangrui.Power system analysis based on real-valued spectral MUSIC algorithm[J].Transactions of China Electrotechnical Society, 2009, 24(12): 149-155.

[5]趙黎麗.基于相關(guān) Hanning窗插值的間諧波分析算法[J].電工技術(shù)學報, 2008, 23(11): 153-158.Zhao Lili.Inter-harmonics analysis based on correlation hanning window and interpolation algorithm[J].Transactions of China Electrotechnical Society, 2008, 23(11): 153-158.

[6]曾博, 滕召勝.納托爾窗改進 FFT動態(tài)諧波參數(shù)估計方法[J].中國電機工程學報, 2010, 30(1): 65-71.Zeng Bo, Teng Zhaosheng.Improved FFT approach for the dynamic state estimation of harmonic parameters based on the Nuttall window[J].Proceeding of the CSEE, 2010, 30(1): 65-71.

[7]呂干云, 方奇品, 蔡秀珊.一種基于粒子群優(yōu)化算法的的間諧波分析方法[J].電工技術(shù)學報, 2009,24(12): 156-161.Lv Ganyun, Fang Qipin, Cai Xiushan.A method for inter-harmonics analysis based on PSO[J].Transactions of China Electrotechnical Society, 2009,24(12): 156-161.

[8]Aiello M, Cataliotti A, Nuccio S.A comparison of spectrum estimation technique for nonstationary signals in induction motor drive measurements[J].IEEE Transaction on Instrumentation and Measurement,2005, 54(6): 2264-2271.

[9]陳涵, 劉會金, 李大路, 等.非均勻采樣和最小二乘法在間諧波檢測中的應用[J].中國電機工程學報,2009, 29(10): 109-114.Chen Han, Liu Huijing, Li Dalu, et al.Application of nonuniformly sampling and least square technique in inter-harmonic measurement[J].Proceeding of the CSEE, 2009, 29(10): 109-114.

[10]Aiello M, Cataliotti A, Nuccio S.A chirp-z transform-based synchronizer for power system measurements[J].IEEE Transaction on Instrumentation and Measurement, 2005, 54(3): 1025-1032.

[11]Kusljevic M D.Simultaneous frequency and harmonic magnitude estimation using decoupled modules and multirate sampling[J].IEEE Transaction on Instrumentation and Measurement, 2010, 59(4): 954-962.

[12]IEC61000-4-7, General guide on harmonics and inter-harmonics measurements and instrumentation,for power supply systems and equipment connected thereto[S].

[13]Bracale A, Carpinelli G, Lauria D, et al.On some spectrum estimation methods for analysis of nonstationary signals in power systems[C].Proceeding of 11th International Conference on Harmonics and Quality of Power, 2004: 260-265.

[14]Guckenheimer J, Holmes P.Nonlinear oscillations,dynamical systems, and bifurcations of vector fields[M].New York: Springer-Verlag, 1991.

[15]Tey L H, So P L, Chu Y C.Improvement of power quality using adaptive shunt active filter[J].IEEE Transaction on Power Delivery, 2005, 20(2):1558-1568.

[16]李輝, 李亦斌, 鄒云屏, 等.一種新的變步長自適應諧波檢測算法[J].電力系統(tǒng)自動化, 2005, 29(2):69-73.Li Hui, Li Yibin, Zou Yunping, et al.A novel adaptive harmonic detecting algorithm based on variable step-size LMS[J].Automation of Electric Power Systems, 2005, 29(2): 69-73.