時變工作點分析算法在上海光源的應用研究

2020-12-18 08:42賴龍偉陳方舟

核技術 2020年12期

楊星賴龍偉陳方舟

1（中國科學院上海應用物理研究所上海201800）

2（中國科學院大學北京100049）

3（中國科學院上海高等研究院上海201204）

工作點是環形粒子加速器非常重要的測量參數，并可用于β 函數、色品、阻抗等多種關鍵物理參數的測量。因此，準確測量工作點對加速器穩定運行和機器優化有重要意義［1］。工作點測量一般指其小數部分的測量，即是對儲存環逐圈位置信號進行頻域計算，提取歸一化中心頻率。對于固定工作點的測量一般使用快速傅里葉變換（Fast Fourier Transform，FFT）算法。因為FFT算法的頻率分辨率由計算點數N 決定，為1/N，有效計算點數越多精度越高，如為了達到歸一化頻率0.001 精度要求有1 000點計算數據，因此一般通過kicker加橫向激勵的方式增加有效計算點數以提高測量精度，但這會影響束流正常運行，只能在機器研究期間使用［2-3］。

上海光源自2012年開始運行恒流注入（top-up）模式，間隔幾分鐘注入一次，儲存環內的束流在注入過程中會引起比較大的橫向振蕩，從而可利用該振蕩無擾測量儲存環的工作點［4-5］。但注入引起的儲存環橫向振蕩衰減很快，只能采集到10 000 多點的逐圈數據，普通的諧波分析方法可以計算得到工作點的分辨率為0.000 1，但此方法無法追蹤注入期間儲存環工作點隨時間的變化。對于儲存環來說，注入過程是一個非常理想的觀測非線性的瞬態過程，因為儲存環中的磁場存在高階項，導致束流偏心位置不同時看到的磁場強度不同，這被稱為儲存環的Lattice 非線性，該效應可以通過測量不同束流偏心位置相應的工作點來進行觀測分析。如果能精確測定束流橫向位置隨時間的變化，測定瞬時工作點隨時間的變化，通過分析工作點與束流位置之間的依賴關系，就可以獲得儲存環非線性效應的強弱信息。

本文將研究更高分辨的時變頻率分析算法，以便對上海光源儲存環注入過程中工作點的時變性質進行精確測量。

1 頻率分析算法

頻譜分析理論隨著數學理論的進步，有幾個比較明顯的階段。常規傅里葉方法及其適用于計算機分析的FFT 在信號處理中占據了極高的地位，連續傅里葉變換的定義如下：

通過改變頻率得到信號的幅頻曲線，即傅里葉譜。但是實踐表明Fourier 頻譜分析并非對所有類型信號的分析都有效，被分析的系統必須是線性的；信號必須是嚴格周期的或者平穩的，否則譜分析結果將缺乏物理意義。但是非線性系統頻率隨時間的變化是十分自然的一種物理現象。因此隨后誕生了短時傅里葉變換（Short-Time Fourier Transform，STFT）。

STFT 將整個長信號在時間軸上根據需求分成有可能重復，也可以不重復的若干個數據組，對數據組進行FFT 分析，從而得到隨時間變化的頻率分析結果。STFT的缺點是需要選擇合適的窗函數，且窗是固定的，這使得短時傅里葉分析的分辨率也是固定的，不能同時滿足頻率分辨率與時間分辨率的要求。此外，STFT作為傅里葉分析的框架下做出的改良，頻率分析能力的上限仍然受制于采樣點數N，且時間分辨率和頻率分辨率無法同時達到一個較優的水準。

小波變換是常用的時變頻率分析算法［6］，國內已有用小波分析進行束流分析的相關研究［7］，可對信號進行局域化分析，通過選擇在時間軸上具有有限能量的適當的基底信號，并和傅里葉系數一樣增添系數建立函數系，對采集到的數據進行不止一維的細化分析，連續小波變換的定義如下［8］：

式中：ψ 函數即小波母函數；τ 為延時系數；a 為尺度系數，可以轉化為頻率，最終得到是系數(a，τ)的二維數組，這就是與傅里葉分析算法在本質上的不同。小波分析使用能量衰減的、在時域影響有限的“小”波信號作為基底，和傅里葉變換一樣通過改變系數來生成一組正交基，因此小波變換能夠更加精確地在局部描述所分析的信號，進而分離出信號不同時間域上的特征。小波展開的結果即小波譜是二維矩陣，兩個方向分別代表了該分量所在時域和頻域的位置，因此比傅里葉變換增加了一維的信息。它通過伸縮平移運算對信號逐步進行多維細化，最終達到不同頻段處時間細分，能自動適應時頻信號分析的要求，從而可聚焦到信號的任意細節，解決了傅里葉變換的頻率分析能力的上限受于采樣點數N制約的困難，也解決了短時傅里葉變換窗口大小固定無法動態調整時頻分辨率的缺點。

1.1 小波算法

本文使用Morlet小波作為小波母函數對數據進行分析，Morlet小波是以Gauss函數作為包絡的復小波，它的解析式為：

圖1為Morlet小波的實數部分。

圖1 Morlet母小波Fig.1 Morlet mother wavelet

Hilbert變換是將采集到的實數域的信號擴展到復數空間上，使之成為解析信號的一種方法，給定一時間連續信號x(t)，其Hilbert變換x?(t)定義為［9］：

x(t)的解析信號定義為：

x(n)的解析信號定義為：

Hilbert變換具有變換前后信號頻譜幅度不發生變化的性質，實部是原信號，虛部與其正交，模的大小是原信號的包絡。通過該方法，我們就可以得到信號幅度隨時間的變化。

對原信號進行Hilbert 變換后，與以頻率為變量構建的時域上的Morlet 小波進行卷積，即可得到信號在該頻率分量上隨時間的相對幅值變化。對所有頻率分量計算之后，即可得到二維小波譜，小波函數的頻率和瞬時頻率越接近時，復小波譜的模就越大，對每一個時間點上的頻譜尋峰，即可近似得到該時間的瞬時頻率。

1.2 NAFF算法

基本頻率的數值分析算法（Numerical Analysis of Fundamental Frequencies，NAFF）是20 世紀90 年代提出來的算法，可極大提高FFT 計算分辨率。引入加速器領域后，常用于工作點分析［10］。

選取一個定長n，按照定長n 將N 點數據分成N-n+1 個短數據列，對每個短數據列進行快速傅里葉分析，得到峰值的窄帶范圍，之后在這個范圍內，再對短數據列按照傅里葉公式：

進行頻率插值，頻率和原頻率越相近，fn的值越大，以此尋找最接近的頻率值作為該點的頻率。

經驗證明，當只采用32 點逐圈數據進行分析時，工作點計算分辨率就可以達到0.001［11］，因此對于長度有限的時變工作點數據，可以通過分段進行NAFF 計算的方法精確計算信號頻率隨時間的變化。

2 算法仿真評估

為了評估兩種算法的性能，使用蒙特卡羅仿真方法對算法進行評估。

測試信號模擬儲存環注入瞬態過程數據，回旋頻率為694 kHz，初始幅度為1，衰減形式為指數型衰減，歸一化工作點從0.220線性變化到0.221，變化方式為掃頻余弦形式。

2.1 分析長度

對于NAFF算法來說，分析長度即窗口的大小，決定了時頻分析的相對精度，窗口越大，時域分辨率越低，頻域分辨率越高。

對于Morlet 小波算法來說，分析長度即小波分析半長度（Wave），也影響著時頻分析的相對精度，Wave越大，時域分辨率越低，頻域分辨率越高。

圖2 是測量不確定度隨分析長度的變化，縱坐標為對數坐標，可以看到，兩種算法都隨著分析長度的增加而提高了頻率分辨率，當分析長度足夠時，兩種方法的分辨率均可達10-6。

圖2 測量不確定度隨分析長度的變化Fig.2 Variation of measurement error with analysis length

2.2 數據信噪比

束團位置信號在采集過程中總是會包含一定程度的噪聲，因此需要對算法在不同信噪比下的分析能力做測試，對兩類算法的隨機誤差進行評估，得到數據的信噪比和頻率不確定性之間的聯系。

圖3 是不同分析長度下測量誤差隨信噪比的變化。

圖3 測量誤差隨信噪比的變化Fig.3 Variation of measurement error with SNR

可以看到，隨著信噪比的增加，兩種方法的頻率精度都有較大的提升，且不同分析長度表現出了相似的變化規律。當分析長度確定時，Morlet 小波算法隨機測量誤差小于NAFF算法。

2.3 系統誤差

不考慮噪聲的影響，因為ejω0t的正交性，兩種方法在數學上都屬無偏分析，但是實際中對于有限的數據點數來說，依靠兩種方法都會產生無法避免的系統誤差，并且，當頻率接近整數或者整數分之一時，系統誤差會較其他頻率有所提升，所幸在工作點設計時，會刻意避免此類數據，因此系統誤差可以很小。

圖4 對兩種方法的頻率依賴性進行了分析，分析長度固定為200，共100 個樣本，可以看出兩種算法的隨機誤差均沒有明顯的頻率依賴性，可以對一般工作點所及的頻率進行良好的分析。

圖4 測量誤差隨分析頻率的關系Fig.4 Variation of measurement error with analysis frequency

2.4 算法比較

將仿真中使用兩種算法進行處理得到的一般性結果進行比較，結果如圖5所示。

圖5 兩種算法的區別 (a)Wavelet，(b)NAFFFig.5 Difference between two algorithm(a)Wavelet,(b)NAFF

結合之前的分析，可以看出小波方法在頻率分辨率和數值穩定性上有著明顯的優勢，信噪比較低時，小波方法的抖動遠小于NAFF 方法。因此選用小波方法進行儲存環工作點的提取工作。

3 算法應用

3.1 算法需求分析

上海光源儲存環在注入時，會引起比較大的束流橫向振蕩然后逐漸衰減，圖6 為利用儲存環上的束流位置探測器（Beam Position Monitor，BPM）電子學設備（Libera Brilliance）采集到的上海光源儲存環注入期間水平方向的逐圈數據。由于加速器四極磁鐵的誤差，注入期間振蕩頻率即工作點也會發生微小的變化，變化范圍比較小。由于衰減時間較短，采集到的逐圈位置數據點數有限，如圖6 所示只有10 000 點，此時使用FFT 算法只能知道整個過程平均工作點為0.225 3，并無更多信息，無法得到頻率的變化。

圖6 上海光源儲存環注入期間逐圈束流水平位置數據Fig.6 Horizontal turn-by-turn beam position data from SSRF storage ring during injection period

束流信號有明顯的特征，即剛注入結束時，振幅較大（信噪比較高），工作點相對變化較大，此時對于測量的要求是高時間分辨率，低頻率分辨率。當振蕩衰減接近噪底時，振幅較?。ㄐ旁氡龋┹^低，此時對于測量的要求是低時間分辨率，高頻率分辨率。

可見，不同的數據段對算法時間分辨率、頻率分辨率要求不完全相同，如想得到盡可能好的結果，需要針對被測對象特征調整算法參數。原則上需要在工作點分辨率滿足要求的前提下，使分析窗口盡可能?。ㄌ岣邥r間分辨率）。

3.2 獲取橫向位置變化曲線

束團橫向位置的測量存在一定的測量誤差，即BPM 系統所用逐圈（Turn-by-Turn，TBT）的分辨率，對注入前后穩定的位置數據進行處理，得到橫向測量位置的系統誤差，逐圈位置測量分辨率為2.1 μm，對應的隨機噪聲峰峰值（Peak-to-Peak，PP）為7.6 μm，如圖7所示。

圖7 位置測量誤差Fig.7 Measurement uncertainty of position

束團橫向位置最大振幅的獲取通過對采集到的振蕩曲線用Hilbert 變換來獲得。因為工作點算法的局域特性，所以求振動曲線的包絡時，采用多次求包絡取平均的做法，得到和工作點計算結果更為匹配的振幅。

圖8 信號及其包絡擬合Fig.8 Signal and its envelope fitting

從圖8 可以看出，Hilbert 算法很好地擬合了束團位置曲線的絕對值，可以認為該值即束團在當時的振幅。且束團在注入過程中的振幅從300 μm 迅速下降到10 μm左右，和上述測量誤差進行比較，可以認為6 000點之后的采樣信噪比過低，分析結果不可信，因此舍去。

3.3 確定信號分析窗口

不同分析長度下有著不同的時頻分辨率，需要在滿足頻域分辨率的條件下，選擇最小的分析長度。

通過對儲存環注入瞬態過程中束團位置信號的進一步分析，使用初始振幅330 μm，噪聲幅度7.6 μm，采樣點數6 000，歸一化振蕩頻率在0.225 25～0.225 35 的仿真信號進行不同分析長度的模擬比較（圖9）?？梢缘玫皆?0-5分辨率下，分析長度最小為200。

圖9 測量誤差隨分析長度的變化Fig.9 Variation of measurement error with analysis length

3.4 分析結果

使用Morlet小波方法對上海光源儲存環注入過程的動態工作點進行分析。橫軸使用束團振幅，縱軸使用動態工作點。以檢驗該過程中的非線性效應。

選取2013 年、2014 年、2015 年、2018 年、2019年、2020 年（每年一組）數據進行處理并比較，得到圖10。

后期數據點較少的原因是期間上海光源進行了注入系統及橫向反饋系統的優化，阻尼時間大大縮短，采樣達到1 000 點時，束團振幅已經下降到50 μm 以下，因此非線性分析僅可用1 000 圈內數據。

可以看出上海光源儲存環存在明顯的非線性行為，可以通過分析注入瞬態過程中的工作點漂移來評估非線性效應的強弱，通過簡單的線性擬合，得到在2013 年4 月29 日，這一值為-0.000 41，在2020 年7 月8 日，這一值為0.000 7。不同時期有不同的表現，總體上非線性效應有增強的趨勢。

4 結語

通過仿真證明，小波分析和NAFF 算法可以進行時變工作點測量，隨機測量誤差小于10-5，其中小波分析算法性能好于NAFF算法。利用Morlet小波算法對上海光源儲存環注入引起的束流橫向振蕩進行分析，可以觀察到工作點隨時間的變化。因此，可以利用小波分析算法和NAFF 算法在top-up供光運行的注入期間進行無擾的工作點測量，進而分析儲存環非線性行為，為加速器的性能優化新增了一個有效工具。