石立華 司 卿 張 琦 張劉輝 郭東義
(解放軍理工大學 電磁環境效應與電光工程國家級重點實驗室,江蘇 南京210007)
在大多數實際工程條件下,屏蔽室屏蔽效能通常采用連續波測量和頻域屏蔽效能表征.實際測量獲得的實質是一種線性系統的幅頻響應特性A(ω).另外一些應用中,需要直接了解屏蔽室對特定脈沖波的屏蔽能力,而工程現場的脈沖波模擬和測量實施不便.解決這一問題可采用頻響函數估計法,即由幅度響應估計系統傳遞函數,進而估計系統脈沖響應.文獻[1]針對這一問題提出了利用希爾伯特變換的方法從幅頻特性A(ω)恢復出相位信息φ(ω),用于建立系統的頻率響應函數.文獻[2]和文獻[3]在此基礎上構造了沖激響應模型h(n),在時域建立了數字濾波器模型.本文進一步研究利用最佳逼近元法及迭代法由幅頻特性建立濾波器模型的方法.實驗證明,采用此方法對于低階頻響函數和屏蔽室脈沖磁場屏蔽效能可以得到很好的估計效果.
若被測系統為線性系統,可用系統的頻響函數描述系統的特性.因此在討論算法之前應假設屏蔽體本身為線性系統或弱非線性系統.而一般來說,屏蔽效能測量的對象在絕大多數情況下都可以在一定的范圍內視為線性系統,因此這個假設是有效的.
假設測得的屏蔽體幅頻特性為|G(ejω)|,則由于沖激響應是實因果信號,其對數幅度譜與相位之間滿足希爾伯特變換關系[1],兩者之間知道其一,就可以完整恢復系統的頻率響應G(ejω).雖然直接應用該模型可以通過快速傅里葉逆變換(FFT)和快速傅里葉逆變換(IFFT)獲得系統的時域響應,但計算過程比較繁瑣,高頻噪聲影響也較大.因此,這里進一步將G(ejω)進行參數化估計,得到系統的離散傳遞函數模型G(z),利用G(z)與差分方程的直接對應關系,可快速計算系統的時域響應[4].
離散傳遞函數模型具有下述形式

式中正整數n和m為此模型分子、分母的階數.對G(z)在單位圓上取值,即可對應于頻響函數G(ejω)為

式中A(ejω)=|G(ejω)|為系統的幅頻響應.在屏蔽室屏蔽效能的測試中,按照有關測量標準獲得的是k個頻點上的幅頻特性,由希爾伯特變換估計相頻響應φ(ω)后,可按式(1)、(2)得到如下方程

式中,θ為待估計的模型參數,

X為k個頻點上頻率響應組成的矩陣,


式(5)中的列向量可以組成線性空間,張成的空間為

則求最佳模型參數的問題就轉變為下列問題:對于s∈A,找到一個最佳逼近元g使得式(9)成立

容易證明兩組列向量 {x0,x1,…,xm},{y0,y1,…,yn}都是線性無關的.對于屏蔽室模型而言,相位函數φ(ωk)一般為非線性的.在φ(ωk)為非線性的條件下,不妨設{x0,x1,…,xm,y0,y1,…,yn}為A的基,為方便起見寫成{ε0,ε1,…,εm+n+1},則最佳逼近元g可以表示如下

最佳逼近元的系數可以由下式求得[5]

由于X列滿秩,所以XTX為對稱正定矩陣,因此,可通過式(11)求得模型參數.
由于測試頻點較多,因此在運算中難以保持矩陣的稀疏性,所以得到初解后還需要使用迭代法進行修正.構造如下同解方程組

設XTX=(pij)n,對于初解θ(0),利用式(12)可構造如下迭代式

式中:P為迭代矩陣;d為修正系數,利用Gauss-Seidel迭代法[6],P 可由下式求得

各量的意義分別為

d與式(12)右側值有關,在此表達式中其值為0.當k足夠大時或誤差達到容許范圍時可停止迭代并令θ(k)作為最終解.
使用迭代法時要保證迭代式收斂[7],否則會進一步加大誤差.當迭代矩陣滿足式(17)時迭代法是收斂的.如果迭代矩陣不滿足此式則不能使用迭代法.

在迭代法收斂的情況下,如果 {θ(k)}收斂于θ*,則可以通過式(18)估計迭代誤差.

式中e為符合精度要求的誤差系數,它可以根據具體情況人為定義.如果式(18)滿足,則表示估計的系數已滿足要求,否則繼續使用高斯牛頓-賽德爾迭代算法對a、b進行修正,直到滿足此式為止.
一般情況下n和m均未知,這樣就無法直接通過上述方法計算出頻響函數系數.可以通過循環檢索的方法,令n和m從2階開始在這一范圍內窮舉取值,并對每一次取值進行估計,最后利用信號相關性比較的方法計算相關系數最大時的階數,即可得到最佳n和m值.但由于是二階循環,進行大范圍的窮舉取值會大大影響計算速度,而對于屏蔽室模型,經過實驗可以驗證當n和m取值超過10階后系統穩定性急劇下降,此時的計算結果與原系統有很大誤差,因此可以將n和m的取值范圍定為2到10.
為了對算法進行驗證,首先進行如下仿真實驗:
1)且構造一個頻率響應已知的數字濾波器,設定一組時域脈沖波形,令其通過該濾波器,得到一組實際輸出波形;
2)對濾波器按其幅頻響應進行模型估計,求出其頻響函數系數并得到估計模型;
3)令設定的輸入波形通過估計模型,得到估計輸出波形并與實際輸出波形比較結果差異.
在此次仿真實驗中所用的標準波形為雙指數波.其表達式為

其各項參數分別為:峰值場強E0=5×104V/m,修正系數k=1.3,脈沖前、后沿參數α=4×107s-1,β=6×108s-1.波形如圖1所示.
將此脈沖波形作為實際輸入波形并令其通過一個歸一化截止頻率為0.2的巴特沃斯濾波器,得到實際輸出波形,如圖2所示.


將由巴特沃斯濾波器的幅度譜(采樣點為1 024)為基礎進行模型估計.經過對比可以得出當n和m分別為3和4時得到的結果最好,將估計的模型與原模型進行比較,如圖3所示.其相頻特性解卷繞后如圖4所示.再令實際輸入的脈沖波形通過估計模型得到估計的輸出波形,并與實際輸出波形進行比較,其結果如圖5所示.

圖3 估計譜與實際譜


可以看出,通過此算法的相位還原,估計得到的模型與實際模型相比結果是一致的.
將對一屏蔽室進行屏蔽效能的實測,并對實測數據進行模型估計.
為了保證時域測量和頻域測量的等效性,就必須確定屏蔽體是否為線性系統.由于實驗所用的屏蔽材料鋼板為磁性材料,其磁化曲線如圖6所示.隨著磁場強度H的增強,在oa階段磁通密度B的增加幾乎可以看作是線性,而當H達到a點時,其磁通密度的增加就不再與H的增加呈線性關系,而當H 達到b點時處于磁飽和狀態[8].由文獻[9]可知,當磁通密度達到0.6T時鋼板材料才進入飽和狀態,而此次實驗測量系統產生的最大磁通密度在1~2mT之間,處于圖中oa階段,遠低于飽和階段臨界點.因此.可以將屏蔽體看作線性系統.
圖7為實驗測試現場.測試時信號源由脈沖磁場模擬器產生脈沖波,發射天線與接收天線相距60 cm,分別在屏蔽門打開和關閉兩種情況下記錄接收波形,由峰值衰減量計算屏蔽效能.利用時域測量結果換算出頻域屏蔽效能,以此作為建立屏蔽室效能模型的初始已知條件.不同頻點下屏蔽室門中心處屏蔽效能數據經過插值預處理后得到的幅頻曲線如圖8所示.由于測量數據量有限,為了提高精度,對采集數據采用了三次立方插值,使之在不改變原始數據所包含信息的基礎上提高頻率分辨率.
對該數據應用本文算法進行分析建模.結果顯示,當n和m均為8時結果最好,估計模型與實際模型對比如圖9所示.相頻特性如圖10所示.





為了檢驗估計的模型是否能用于對時域脈沖波輸出波形進行估計,在同樣的實驗條件下利用時域脈沖發生器產生一組脈沖波作為激勵,測出通過該屏蔽室后的實際輸出波形,同時將該激勵通過估計的屏蔽室模型,比較估計的輸出波形與實際波形的差異.圖11為實測的脈沖輸入波形與經過屏蔽室后得到的輸出波形.圖12為估計的結果與實測輸出波形比較.從圖中可以看出,兩者比較接近.

圖11 實測輸入/輸出波形

圖12 估計結果與實測結果比較
受測試設備和現場條件的限制,開展屏蔽室對脈沖波屏蔽效能的測試存在一定的困難.本文方法為由連續波測量結果評估脈沖波響應提供了可行方案.對于屏蔽室這樣一個線性模型,其對脈沖磁場的響應近似為一個低通濾波器,可由頻域測量數據建立其離散傳遞函數模型,再由該模型估計其脈沖波響應.該方法特別適用于解決屏蔽室對特定電磁脈沖波屏蔽效果的評估問題.
[1]TESCHE F M.On the use of Hilbert transform for processing measure CW data[J].IEEE Transforms on Electromagnetic Compatibility,1992,34(3):259-266.
[2]石立華,周璧華,陳 彬,等.基于幅-頻曲線的系統時域響應特性評價方法[J].電波科學學報,2000,15(4):465-471.SHI Lihua,ZHOU Bihua,CHEN Bin,et al.Time domain characterization of a system based on the magnitude of its frequency response[J].Chinese Journal of Radio Science,2000,15(4):465-471.(in Chinese)
[3]郭東義,石立華,郭曜華,等.基于頻響函數測量脈沖磁場屏蔽效能的新方法[J].信息與電子工程,2010,2(1):49-52.GUO Dongyi,SHI Lihua,GUO Yaohua,et al.Measurement for shielding effectiveness of pulsed magnetic field based on frequency response function[J].Information and Electronic Engineering,2010,2(1):49-52.(in Chinese)
[4]周壁華,陳 彬,石立華.電磁脈沖及其工程防護[M].北京:國防工業出版社,2003.
[5]TOROKHTI A P,HOWLETT P G.Best approximation of operators in the modeling of nonlinear systems[J].Circuits and Systems I:Fundamental Theory and Applications,2002,49(12):1792-1798.
[6]MOGHNIEH H,LOWTHER D A.The solution of electromagnetic field problems using a sliding window Gauss-Seidel algorithm on a multicore processor[J].IEEE Transactions on Magnetics,2010,46(8):3081-3084.
[7]BIEN Z,HUH K M.Higher-order iterative learning control algorithm[J].Control Theory and Applications,1989,136(3):105-112.
[8]ENDO H,HAYANO S,SAITO Y,et al.Magnetization curve plotting from the magnetic domain images[J].IEEE Transactions on Magnetics,2001,37(4):2727-2730.
[9]KLINKENBUSCH L.On the shielding effectiveness of enclosures[J].IEEE Transforms on Electromagnetic Compatibility,2005,47(3):589-601.