楊 星 賴龍偉 陳方舟

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

2（中國科學(xué)院大學(xué) 北京100049）

3（中國科學(xué)院上海高等研究院 上海201204）

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

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

本文將研究更高分辨的時(shí)變頻率分析算法，以便對上海光源儲存環(huán)注入過程中工作點(diǎn)的時(shí)變性質(zhì)進(jìn)行精確測量。

1 頻率分析算法

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

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

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

小波變換是常用的時(shí)變頻率分析算法［6］，國內(nèi)已有用小波分析進(jìn)行束流分析的相關(guān)研究［7］，可對信號進(jìn)行局域化分析，通過選擇在時(shí)間軸上具有有限能量的適當(dāng)?shù)幕仔盘枺⒑透道锶~系數(shù)一樣增添系數(shù)建立函數(shù)系，對采集到的數(shù)據(jù)進(jìn)行不止一維的細(xì)化分析，連續(xù)小波變換的定義如下［8］：

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

1.1 小波算法

本文使用Morlet小波作為小波母函數(shù)對數(shù)據(jù)進(jìn)行分析，Morlet小波是以Gauss函數(shù)作為包絡(luò)的復(fù)小波，它的解析式為：

圖1為Morlet小波的實(shí)數(shù)部分。

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

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

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

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

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

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

1.2 NAFF算法

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

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

進(jìn)行頻率插值，頻率和原頻率越相近，fn的值越大，以此尋找最接近的頻率值作為該點(diǎn)的頻率。

經(jīng)驗(yàn)證明，當(dāng)只采用32 點(diǎn)逐圈數(shù)據(jù)進(jìn)行分析時(shí)，工作點(diǎn)計(jì)算分辨率就可以達(dá)到0.001［11］，因此對于長度有限的時(shí)變工作點(diǎn)數(shù)據(jù)，可以通過分段進(jìn)行NAFF 計(jì)算的方法精確計(jì)算信號頻率隨時(shí)間的變化。

2 算法仿真評估

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

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

2.1 分析長度

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

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

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

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

2.2 數(shù)據(jù)信噪比

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

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

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

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

2.3 系統(tǒng)誤差

不考慮噪聲的影響，因?yàn)閑jω0t的正交性，兩種方法在數(shù)學(xué)上都屬無偏分析，但是實(shí)際中對于有限的數(shù)據(jù)點(diǎn)數(shù)來說，依靠兩種方法都會產(chǎn)生無法避免的系統(tǒng)誤差，并且，當(dāng)頻率接近整數(shù)或者整數(shù)分之一時(shí)，系統(tǒng)誤差會較其他頻率有所提升，所幸在工作點(diǎn)設(shè)計(jì)時(shí)，會刻意避免此類數(shù)據(jù)，因此系統(tǒng)誤差可以很小。

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

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

2.4 算法比較

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

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

結(jié)合之前的分析，可以看出小波方法在頻率分辨率和數(shù)值穩(wěn)定性上有著明顯的優(yōu)勢，信噪比較低時(shí)，小波方法的抖動遠(yuǎn)小于NAFF 方法。因此選用小波方法進(jìn)行儲存環(huán)工作點(diǎn)的提取工作。

3 算法應(yīng)用

3.1 算法需求分析

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

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

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

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

3.2 獲取橫向位置變化曲線

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

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

束團(tuán)橫向位置最大振幅的獲取通過對采集到的振蕩曲線用Hilbert 變換來獲得。因?yàn)楣ぷ鼽c(diǎn)算法的局域特性，所以求振動曲線的包絡(luò)時(shí)，采用多次求包絡(luò)取平均的做法，得到和工作點(diǎn)計(jì)算結(jié)果更為匹配的振幅。

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

從圖8 可以看出，Hilbert 算法很好地?cái)M合了束團(tuán)位置曲線的絕對值，可以認(rèn)為該值即束團(tuán)在當(dāng)時(shí)的振幅。且束團(tuán)在注入過程中的振幅從300 μm 迅速下降到10 μm左右，和上述測量誤差進(jìn)行比較，可以認(rèn)為6 000點(diǎn)之后的采樣信噪比過低，分析結(jié)果不可信，因此舍去。

3.3 確定信號分析窗口

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

通過對儲存環(huán)注入瞬態(tài)過程中束團(tuán)位置信號的進(jìn)一步分析，使用初始振幅330 μm，噪聲幅度7.6 μm，采樣點(diǎn)數(shù)6 000，歸一化振蕩頻率在0.225 25～0.225 35 的仿真信號進(jìn)行不同分析長度的模擬比較（圖9）。可以得到在10-5分辨率下，分析長度最小為200。

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

3.4 分析結(jié)果

使用Morlet小波方法對上海光源儲存環(huán)注入過程的動態(tài)工作點(diǎn)進(jìn)行分析。橫軸使用束團(tuán)振幅，縱軸使用動態(tài)工作點(diǎn)。以檢驗(yàn)該過程中的非線性效應(yīng)。

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

后期數(shù)據(jù)點(diǎn)較少的原因是期間上海光源進(jìn)行了注入系統(tǒng)及橫向反饋系統(tǒng)的優(yōu)化，阻尼時(shí)間大大縮短，采樣達(dá)到1 000 點(diǎn)時(shí)，束團(tuán)振幅已經(jīng)下降到50 μm 以下，因此非線性分析僅可用1 000 圈內(nèi)數(shù)據(jù)。

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

圖10 不同運(yùn)行時(shí)期數(shù)據(jù)的分析結(jié)果 (a)2013，(b)2014，(c)2015，(d)2018，(e)2019，(f)2020Fig.10 Analysis results of data in different operation periods (a)2013,(b)2014,(c)2015,(d)2018,(e)2019,(f)2020

4 結(jié)語

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