熊 云 閻映炳 冷用斌 易 星 賴龍偉王寶鵬 張 寧 陳之初

1（中國科學(xué)院上海應(yīng)用物理研究所 上海 201800）

2（中國科學(xué)院研究生院 北京 100049）

上海光源(Shanghai Synchrotron Radiation Facility, SSRF)自2009年5月在Decay模式下投入運(yùn)行以來，有力地推動(dòng)了我國重要前沿領(lǐng)域基礎(chǔ)研究和高新技術(shù)開發(fā)應(yīng)用研究。但是，為更好地滿足光源用戶的需求，SSRF須產(chǎn)生更高亮度、更高穩(wěn)定性的束團(tuán)，這就需要采用恒流注入模式(Top-up Mode)，即宏觀上束流流強(qiáng)不隨時(shí)間衰減，微觀上束流流強(qiáng)好于 1%的恒定水平。同時(shí)，這也需要高精度、高時(shí)效的直流流強(qiáng)檢測(cè)系統(tǒng)對(duì)恒流注入實(shí)施診斷。本文介紹恒流模式下上海光源儲(chǔ)存環(huán)直流流強(qiáng)檢測(cè)系統(tǒng)。

1 原直流流強(qiáng)檢測(cè)系統(tǒng)

原直流流強(qiáng)檢測(cè)系統(tǒng)基于Bergoz NPCT型探頭和 PXI數(shù)據(jù)采集處理平臺(tái)，系統(tǒng)軟件在 Windows平臺(tái)上采用 LabVIEW 圖形化編程語言進(jìn)行開發(fā)，通過Shared Memory IOCcore技術(shù)來實(shí)現(xiàn)EPICS的數(shù)據(jù)接口，完成 LabVIEW 應(yīng)用程序和控制系統(tǒng)之間的數(shù)據(jù)交換[1]。

原系統(tǒng)計(jì)算流強(qiáng)的算法為算術(shù)平均濾波法：以10 kHz數(shù)據(jù)采樣率，累積采樣0.3 s，取其算術(shù)平均為直流流強(qiáng)值。算術(shù)平均濾波算法較為簡(jiǎn)單，適用于對(duì)帶有隨機(jī)干擾的信號(hào)進(jìn)行濾波。

圖1(a)為原系統(tǒng)流強(qiáng)分辨率，束流流強(qiáng)小于300 mA時(shí)，流強(qiáng)分辨率優(yōu)于1.5 μA。圖1(b)為原系統(tǒng)采樣數(shù)據(jù)剔除直流分量后的頻譜圖，束流的直流分量為200 mA，可見束流平均流強(qiáng)為200 mA時(shí)在483.63 Hz處有一個(gè)較大單頻干擾信號(hào)。由于此干擾信號(hào)的存在，用算術(shù)平均濾波算法效果不理想，會(huì)與真值產(chǎn)生偏差，故須優(yōu)化原束流流強(qiáng)算法。

圖1 束流流強(qiáng)分辨率(a)和原始信號(hào)噪聲頻譜(b)Fig.1 The measured resolution of beam current (a) and noise pattern of the original signal (b).

再則，基于 Windows平臺(tái)的直流流強(qiáng)測(cè)量系統(tǒng)，因軟件病毒及網(wǎng)絡(luò)資源泄漏等問題，長(zhǎng)期運(yùn)行穩(wěn)定性欠佳，平均無故障時(shí)間低于1個(gè)月。為提高系統(tǒng)平均無故障時(shí)間，須優(yōu)化原有系統(tǒng)。

該Decay模式下的直流流強(qiáng)測(cè)試系統(tǒng)的硬件分辨率達(dá)μA量級(jí)，完全滿足恒流模式下對(duì)直流流強(qiáng)監(jiān)測(cè)系統(tǒng)的要求，可沿用原的系統(tǒng)硬件配置：(1) 法國Bergoz Instrumentation公司第三代的NPCT175；(2) 直流流強(qiáng)變壓器前端電子學(xué)(DCCT Front-end electronics)；(3) NI PXI-4070數(shù)字萬用表；4) NI PXI-8196 2.0 GHz CPU控制器。

2 軟件優(yōu)化

SSRF用EPICS (Experimental Physics and Industrial Control System)作為的控制系統(tǒng)，其有基于工具集的開發(fā)環(huán)境，我們將其作為軟件開發(fā)平臺(tái)。由于Linux操作系統(tǒng)的高可靠性、高擴(kuò)展性以及優(yōu)異的網(wǎng)絡(luò)性能，考慮系統(tǒng)的穩(wěn)定性以及與EPICS的兼容性需要，我們采用Linux操作系統(tǒng)。

系統(tǒng)的軟件結(jié)構(gòu)如圖2所示。PXI-4070數(shù)據(jù)采集卡處于系統(tǒng)的最底層，也是整個(gè)上層軟件的控制對(duì)象。在安裝了Linux操作系統(tǒng)、采集卡的底層驅(qū)動(dòng)(可通過命令 nilsdev查看是否檢測(cè)到硬件)以及EPICS開發(fā)環(huán)境之后，即可在新建的IOC下根據(jù)用戶手冊(cè)編寫符合實(shí)際需求的功能配置模塊。隨后在IOC中還需對(duì)采集到的信號(hào)進(jìn)行數(shù)據(jù)處理，處理結(jié)果通過Database及Channel Access實(shí)現(xiàn)與其它IOC及中控室OPI的通信。

軟件開發(fā)包括：(1) 束流流強(qiáng)數(shù)據(jù)采集。通過調(diào)用PXI-4070數(shù)字萬用表驅(qū)動(dòng)程序的API函數(shù)，配置PXI-4070數(shù)字萬用表的初始狀態(tài)，采集儲(chǔ)存環(huán)束流流強(qiáng)信號(hào)；(2) 應(yīng)用層軟件。在EPICS框架下編寫相應(yīng)的設(shè)備支持模塊(device support module)，針對(duì)直流流強(qiáng)檢測(cè)系統(tǒng)編制其相應(yīng)的運(yùn)行數(shù)據(jù)庫，并對(duì)采集的原始束流流強(qiáng)數(shù)據(jù)進(jìn)行數(shù)據(jù)處理，得到束流流強(qiáng)、束流壽命以及置信度等信息，具體的數(shù)據(jù)處理方法將在算法部分給出；(3) OPI(操作員界面)的編寫。將運(yùn)行數(shù)據(jù)庫中的束流流強(qiáng)以及束流壽命等需要監(jiān)控的信息顯示在中控室的人機(jī)界面上。軟件配置為：Mandriva Linux 2006；BaseR3.14.10；NI-DMM 2.5。

圖2 系統(tǒng)軟件結(jié)構(gòu)Fig.2 Software structure of DCCT system.

3 算法優(yōu)化

3.1 信號(hào)模型

分析SSRF儲(chǔ)存環(huán)直流流強(qiáng)檢測(cè)系統(tǒng)采集的流強(qiáng)數(shù)據(jù)，流強(qiáng)信號(hào)的組成為：PXI4070與NPCT175的聯(lián)機(jī)噪聲信號(hào)；束流的流強(qiáng)信號(hào)，其隨時(shí)間的變化可描述為：I(t)=I0exp(–t/τ)，其中，τ為束流壽命；I0為初始流強(qiáng)(設(shè)定為200 mA)。于是，束流流強(qiáng)信號(hào)的模型為：Is(i)=I0exp{–i/(fsτ0)}，i=0,1,…,Tfs, 其中τ0為束流真實(shí)壽命，Is為流強(qiáng)采樣值，fs為采樣頻率(設(shè)定為10 kHz)，T為累積采樣時(shí)間(設(shè)定為0.3 s)。由此計(jì)算得到一組無測(cè)量誤差的流強(qiáng)采樣數(shù)據(jù)。

在每個(gè)采樣數(shù)據(jù)上疊加一個(gè)隨機(jī)測(cè)量誤差noise1(i)=A1random(1,TFs)，不同采樣頻率下，噪聲水平不同。因?yàn)楹懔髂Ｊ紻CCT系統(tǒng)設(shè)定的DVM的采樣頻率為fs=10 kHz，則noise1(i)=0.016 random(1,TFs)[1]；在上述數(shù)據(jù)基礎(chǔ)上疊加一個(gè)單頻干擾信號(hào)采樣值 noise2(i)=A2sin(2πif/fs+ψ)，其中A2為干擾信號(hào)振幅，f為單頻干擾信號(hào)頻率，ψ為初始相位。

3.2 束流流強(qiáng)算法

針對(duì)原系統(tǒng)的算術(shù)平均濾波算法，在不同的單頻干擾振幅及單頻干擾頻率條件下，對(duì)計(jì)算值與真實(shí)值的偏差值進(jìn)行了 Matlab仿真，仿真結(jié)果如圖3(a)所示?？紤]原始信號(hào)中的隨機(jī)測(cè)量噪聲以及單頻干擾信號(hào)，當(dāng)累積采樣時(shí)間內(nèi)剛好采集整數(shù)周期的單頻干擾信號(hào)時(shí)，單頻干擾信號(hào)將對(duì)算術(shù)平均值沒有影響。所以對(duì)累積采樣時(shí)間內(nèi)采集的單頻干擾信號(hào)最大整數(shù)周期的原始數(shù)據(jù)進(jìn)行算術(shù)平均，而不是簡(jiǎn)單的對(duì)累積采樣時(shí)間內(nèi)采集的所有數(shù)據(jù)進(jìn)行算術(shù)平均，這樣可以減少與真值偏差。針對(duì)不同頻率以及振幅的單頻干擾信號(hào)，Matlab仿真結(jié)果如圖3(b)所示。對(duì)整數(shù)周期算術(shù)平均算法而言，無論單頻干擾信號(hào)的幅度多大，都可以保證束流流強(qiáng)偏差值小于1 μA。仿真程序中做FFT時(shí)所用點(diǎn)數(shù)為262 144，補(bǔ)0做FFT主要是為了提高計(jì)算分辨率，使求得的單頻干擾頻率更準(zhǔn)確，從而使進(jìn)行算術(shù)平均計(jì)算所需的原始數(shù)據(jù)點(diǎn)數(shù)更準(zhǔn)確。針對(duì)整數(shù)周期算術(shù)平均算法，仿真不同振幅及頻率的單頻干擾信號(hào)對(duì)束流流強(qiáng)分辨率的影響，仿真結(jié)果如圖3(c)所示。不同振幅以及頻率的單頻干擾信號(hào)下，束流流強(qiáng)的分辨率都好于0.38 μA。

圖3 算術(shù)平均算法的偏差值(a)、整數(shù)周期算術(shù)平均算法的偏差值(b)和整數(shù)周期算術(shù)流強(qiáng)分辨率(c)Fig.3 Deviations of arithmetic mean algorithm (a) and the integer cycle arithmetic mean algorithm (b),and resolution of integer cycle arithmetic mean (c).

3.3 束流壽命算法

Top-up模式下直流流強(qiáng)檢測(cè)系統(tǒng)的束流壽命算法為曲線擬合算法，與 Decay模式相同，但由于Top-up模式下束流的特性不同于Decay模式，故用圖4 中的(t1,I1)、(t2,I2)兩點(diǎn)進(jìn)行曲線擬合，于是

針對(duì)兩點(diǎn)曲線擬合求束流壽命的方法，Matlab仿真結(jié)果如圖5(a)所示。當(dāng)單頻干擾信號(hào)的振幅小于1 mA時(shí)，束流壽命的絕對(duì)偏差小于0.04 h。單頻干擾信號(hào)的振幅為1 mA時(shí)，束流壽命的分辨率仿真結(jié)果如圖5(b)所示。不同頻率的單頻干擾信號(hào)下，束流壽命的分辨率好于0.018 h。

圖4 曲線擬合示意圖Fig.4 The curve fitting schemes.

圖5 曲線擬合束流壽命偏差(a)和束流壽命分辨率(b)Fig.5 Beam lifetime error of curve-fitting (a) and resolution of beam lifetime (b).

3.4 置信度算法

引入置信度算法，來確定直流流強(qiáng)檢測(cè)系統(tǒng)每次給出的束流平均流強(qiáng)值以及壽命值的可信度，每次系統(tǒng)根據(jù)當(dāng)前的狀態(tài)來預(yù)測(cè)下一個(gè)數(shù)據(jù)，然后同依據(jù)采樣數(shù)據(jù)計(jì)算得到的束流流強(qiáng)值以及壽命值進(jìn)行比較，得差值，依據(jù)此差值給出實(shí)測(cè)值的可信度。置信度算法中主要是數(shù)據(jù)預(yù)測(cè)算法。我們使用Kalman算法，它是最優(yōu)化自回歸數(shù)據(jù)處理算法，無需保存先前的數(shù)據(jù)，進(jìn)行新的測(cè)量時(shí)，也無需對(duì)原來數(shù)據(jù)進(jìn)行處理。直流流強(qiáng)檢測(cè)系統(tǒng)的 Kalman算法為：

(1)X(k|k–1)=2X(k|k) –X(k–1|k–1)，其中X(k|k–1)是利用上一狀態(tài)預(yù)測(cè)的結(jié)果，X(k–1|k–1)是上一狀態(tài)最優(yōu)結(jié)果；

(2)P(k|k–1)=P(k–1|k–1)+Q，其中P(k|k–1)是X(k|k–1)對(duì)應(yīng)的協(xié)方差，P(k–1|k–1)是X(k–1|k–1)對(duì)應(yīng)的協(xié)方差，Q為過程噪聲的方差；

(3)Kg(k)=P(k|k–1)/[P(k|k–1)+R]，其中Kg為卡爾曼增益(Kalman Gain)，R為量測(cè)噪聲的方差；

(4)X(k|k)=X(k|k–1)+Kg(k)[Z(k)–X(k|k–1)]，其中Z(k)是k時(shí)刻的實(shí)測(cè)值；

(5)P(k|k)=(1–Kg(k))P(k|k–1)，其中P(k|k)是X(k|k)對(duì)應(yīng)的協(xié)方差；

(6)E(k|k)=X(k|k–1)–Z(k)，其中E(k|k)是k時(shí)刻的預(yù)測(cè)值與實(shí)測(cè)值的差值。

給定下列變量的初始值：量測(cè)噪聲的方差R，過程噪聲的方差 Q，第一點(diǎn)的預(yù)測(cè)值X(1|0)，第一點(diǎn)的協(xié)方差P(1|1)，Kalman算法就會(huì)自動(dòng)收斂。利用某次實(shí)測(cè)的數(shù)據(jù)進(jìn)行Matlab仿真，數(shù)據(jù)預(yù)測(cè)算法仿真結(jié)果如圖6(a)所示。置信度的仿真結(jié)果如圖6(b)所示。

圖6 數(shù)據(jù)預(yù)測(cè)算法仿真結(jié)果(a)和置信度仿真結(jié)果(b)Fig.6 Simulated result of data prediction algorithm (a) and simulation result of confidence level (b).

對(duì)比分析圖6(a)與6(b)，可分為以下幾個(gè)階段：

(1) 收斂階段(1–5)：第1個(gè)預(yù)測(cè)值是初值，與實(shí)測(cè)值差值大于160 μA，置信度為0，到第5個(gè)預(yù)測(cè)值時(shí)，與實(shí)測(cè)值差值小于10 μA，預(yù)測(cè)值的置信度為1；

(2) 穩(wěn)定階段(6–454)：束流自然衰減，預(yù)測(cè)值與實(shí)測(cè)值差值小于10 μA，預(yù)測(cè)值的置信度為1；

(3) 突變階段一(455–634)：第 455個(gè)測(cè)量點(diǎn)，儲(chǔ)存環(huán)掉束，預(yù)測(cè)值與實(shí)測(cè)值差值大于160 μA，置信度為0；從456點(diǎn)到634點(diǎn)，儲(chǔ)存環(huán)中并無束流，預(yù)測(cè)值與實(shí)測(cè)值差值小于10 μA，預(yù)測(cè)值的置信度為1；

(4) 突變階段二(635–653)：束流注入過程，束流變換流強(qiáng)變化快，預(yù)測(cè)值跟隨實(shí)測(cè)值變化；

(5) 收斂階段(654–660)：預(yù)測(cè)值與實(shí)測(cè)值差值慢慢減小，直到穩(wěn)定在10 μA之內(nèi)。

重復(fù)上述五個(gè)階段，其中對(duì)于第三個(gè)階段，若無掉束現(xiàn)象，第三個(gè)階段就不存在。

3.5 軟件實(shí)現(xiàn)

中央控制室需顯示的參數(shù)是束流平均流強(qiáng)以及壽命等參數(shù)。Top-up模式下的束流流強(qiáng)測(cè)量以及計(jì)算方法通用于Decay模式，唯束流壽命算法不同。SSRF現(xiàn)尚運(yùn)行于Decay模式，故選擇兩種模式結(jié)合的實(shí)施方案：在現(xiàn)有的軟件模塊上增加新的模塊，實(shí)現(xiàn)Top-up模式下束流壽命的計(jì)算。軟件的流程圖如圖7所示。

圖7 直流流強(qiáng)檢測(cè)系統(tǒng)軟件流程圖Fig.7 Flow chart of DCCT system.

4 結(jié)語

介紹了恒流模式下上海光源直流流強(qiáng)檢測(cè)系統(tǒng)的軟硬件結(jié)構(gòu)。針對(duì)恒流模式下儲(chǔ)存環(huán)束流的特性，對(duì)束流流強(qiáng)以及束流壽命算法進(jìn)行了實(shí)驗(yàn)室仿真。提出了計(jì)算束流流強(qiáng)的整數(shù)周期算術(shù)平均算法以及計(jì)算束流壽命的兩點(diǎn)曲線擬合法。引入了置信度算法對(duì)系統(tǒng)給出的束流平均流強(qiáng)值以及壽命值進(jìn)行評(píng)估，實(shí)驗(yàn)室仿真結(jié)果顯示，該算法準(zhǔn)確、快速、高效。最后，給出了直流流強(qiáng)檢測(cè)系統(tǒng)的軟件流程圖，實(shí)現(xiàn)了上海光源Decay與Top-up運(yùn)行模式的統(tǒng)一。

目前，優(yōu)化后的系統(tǒng)已在上海光源穩(wěn)定地運(yùn)行，下一步的重點(diǎn)工作是對(duì)代碼的優(yōu)化，提高代碼的執(zhí)行效率。對(duì)算法可能存在的缺陷進(jìn)行優(yōu)化，提高算法的效率、健壯性，提升整個(gè)系統(tǒng)的性能。

1 冷用斌, 周偉民, 陳永忠, 等. 上海光源儲(chǔ)存環(huán)直流流強(qiáng)監(jiān)測(cè)系統(tǒng)設(shè)計(jì)[J]. 核技術(shù), 2007, 30(6): 477–480 LENG Yongbin, ZHOU Weimin, CHEN Yongzhong,et al. DCCT system design for SSRF storage ring[J]. Nucl Tech, 2007, 30(6): 477–480

2 Sreenivas T V, Rao P V S. FFT algorithm for both input and out—put pruning[J]. IEEE Transactions on Acoustics Speech and Signal Processing, 1979, ASSP-27(3): 291–292

3 Nagal K. Pruning the decimation in-time FFT algorithm with frequency shift[J]. IEEE Transactions on Acoustics Speech and Signal Processing, 1986, ASSP-34(4): 1008–10101

4 Schütte W. Algorithm for the lifetime/loss measurement with a high precision DC transformer[C], 4th European Particle Accelerator Conference, London, UK, 27 Jun-1 Jul 1994, 1631

5 Sreenivas T V, Rao P V S. High-resolution narrow-band spectra by FFT running[J]. IEEE Transactions on Acoustics Speech and Signal Processing, 1980, ASSP-28(2): 254–257

6 胡廣書. 數(shù)字信號(hào)處理理論、算法與實(shí)現(xiàn). 第2版[M].北京: 清華大學(xué)出版社, 2003: 169–192 HU Guangshu. Digital signal processing theory, algorithm and implementation. 2nd ed[M]. Beijing: Tsinghua University Press, 2003: 169–192

7 Kalman filter http://en.wikipedia.org/wiki/Kalman_Filter [OL]