基于兩級壓縮感知的脈沖星時延估計方法?

康志偉 吳春艷 劉勁 馬辛 桂明臻

1)(湖南大學信息科學與工程學院,長沙 410082)

2)(武漢科技大學信息科學與工程學院,武漢 430081)

3)(北京航空航天大學儀器科學與光電工程學院,北京 100191)

1 引 言

X射線脈沖星導航作為一種前瞻性的航天器自主天文導航方法,能為近地軌道、深空和星際飛行的航天器提供位置、速度、時間、姿態等豐富的導航信息[1,2].天文臺通過長期的天文觀測,可獲得高信噪比的脈沖星標準脈沖輪廓,器載X射線探測器通過短時間觀測采集X射線脈沖星信號,并將這些信號按照周期累積,可得到脈沖星累積脈沖輪廓,其信噪比較低.X射線脈沖星導航的基本量測量是脈沖時延,它通過比較脈沖星累積脈沖輪廓相對于標準脈沖輪廓的相位得到[3,4].如何快速精確地進行累積脈沖輪廓時延估計是提高脈沖星自主導航性能的關鍵[5].

脈沖星累積脈沖輪廓的時延估計方法主要分為頻域與時域兩類[6?9].頻域方法所需的變換、計算相對復雜,而時域方法較為直觀、計算量小、實時性較好,已得到研究者關注.Emadzadeh和Speyer[10]提出了在時域利用歷元折疊得到累積脈沖輪廓后,再通過設計累積脈沖輪廓與標準脈沖輪廓的最小二乘目標函數,以獲得相位值的時延估計方法.另外,Emadzadeh和Speyer[11]還提出了利用累積脈沖輪廓與標準脈沖輪廓間的互相關函數的最大值來估計時延.李建勛和柯熙政[12]提出了利用到達時間概率函數構建最大似然準則對脈沖星累積脈沖輪廓進行時延估計的方法.

壓縮感知[13]是一種時域的信號處理方法,它基于信號的稀疏性,通過極低的采樣率獲取的信息來重構原始信號,近年來已被應用于累積脈沖輪廓時延估計.蘇哲等[14]提出了一種基于模信轉換器的觀測矩陣,以構建能在短時間內得到高信噪比累積脈沖輪廓的壓縮感知方法.Li等[15]將脈沖探測過程模擬成對累積脈沖輪廓的隨機采樣,以隨機0-1矩陣作為測量矩陣實現了一種能有效提高實時性的壓縮感知方法.Shen等[16]利用Hadamard矩陣正交性好以及其元素為?1與+1時容易硬件實現和易存儲的特點,提出了一種基于Hadamard測量矩陣的累積脈沖輪廓壓縮感知方法.

以上累積脈沖輪廓時延估計主要是基于測量矩陣及其改進的壓縮感知方法,具有較好的時延估計性能.考慮到壓縮感知方法中字典的原子數與字典精度密切相關,原子個數越多則原子間隔越小、估計越精確,但相應的計算量也越大.因此,從波形字典的角度來考慮,實現具有高精度、低計算復雜度的壓縮感知累積脈沖輪廓時延估計方法是一條值得探討的途徑.

針對傳統字典中原子數增加雖可提高估計精度但會增加計算量這一問題,本文基于波形字典的可分解性,構建粗估計與精估計兩級字典,探尋運用兩級字典壓縮感知的脈沖星累積脈沖輪廓時延估計方法,以實現能進行全局估計與局部估計的動態組合,達到具有高精度和低計算復雜度的時延估計性能.

2 X射線脈沖星信號輪廓模型

X射線脈沖星信號輪廓是通過長時間天文觀測獲取脈沖星輻射數據,并對其進行脈沖周期折疊和同步平均而得到,其本質是探測到的光子流量密度與相位或時間的曲線關系[2].脈沖星標準脈沖輪廓是一組具有超高信噪比的數據,累積脈沖輪廓是在某一觀測時間段累積得出的數據.假設觀測時間間隔為(u0,uf),總的觀測時間為Tobs=uf?u0,定義ui為第i個光子的到達時間,遞增的到達時間集合{u1,u2,···,uM}可表示成u06u16u26···6uM6uf,設u0=0,Cu為(0,uf)時間間隔內探測到的光子數量,它是具有時變速率λ(u)>0的非齊次泊松過程.對于固定的時間uf,Cu是泊松過程的隨機變量,在時間間隔(0,u)內探測到m個光子的概率可表示為[11]

其均值和方差分別為

標準脈沖輪廓模型,即總流量密度函數λ(u),由X射線脈沖星的光子密度及來自背景環境的光子密度構成:

式中λb和λs分別代表已知有效的背景流量和脈沖星源的流量;h(φ)為歸一化周期脈沖輪廓,具有以下特性[17]:minφh(φ)=0,φ∈[0,1);φdet(u)為探測器相位.

對于同一顆脈沖星,累積脈沖輪廓與標準脈沖輪廓的關系可表示如下[2]:

式中f(u)為累積脈沖輪廓,λ(u)為標準脈沖輪廓,n(u)為噪聲,c為直流分量,l為放大系數,τ為累積脈沖輪廓相對于標準脈沖輪廓的時間延遲值.

實際提取出的脈沖星標準脈沖輪廓與累積脈沖輪廓是離散序列.標準脈沖輪廓在一個周期P中有N個采樣點,每個采樣點的間隔為?t=P/N,則累積脈沖輪廓與標準脈沖輪廓關系為

i為累積脈沖輪廓相對于標準脈沖輪廓的延遲間隔數,脈沖時延估計量為

3 兩級壓縮感知的時延估計算法

壓縮感知作為一種新采樣理論,利用了信號的稀疏性,在低采樣率的情況下,通過降維來獲取信號的離散樣本以重建信號.壓縮感知主要包括三個部分:信號的稀疏表示、測量矩陣采樣和信號重構.

壓縮感知以其低計算復雜度與累積時間短而優于其他方法,近年來已被用于脈沖星時延估計.字典是一個影響壓縮感知估計精度的重要因素.字典中原子數越多,估計精度越高,但同時也增加了計算復雜度.為此,本文采用兩級字典使其在提高估計精度的同時降低計算復雜度.兩級字典包括粗估計字典與精估計字典,粗估計字典原子間隔大,能提供全相位估計,精估計字典原子間隔小,且個數少,能提供局部估計.利用粗估計字典對累積脈沖輪廓時延值進行預估計,將其相位調整到精估計字典的估計范圍內,再利用精估計字典得到精確時延值.

3.1 兩級字典的脈沖輪廓稀疏表示

字典選擇在壓縮感知稀疏過程中起著至關重要的作用,由(4)式可知,累積脈沖輪廓和標準輪廓經歸一化后,僅相差一個時間延遲量.因此,對X射線脈沖輪廓,可根據標準脈沖輪廓λ(u)設計出粗估計與精估計兩級字典.

第一級由粗估計字典提供全局相位預估計.粗估計字典包含脈沖輪廓的所有相位,采樣點數越多、字典原子間隔越小、時延估計值越精確,但這種全局估計會增加計算量.因此,第一級可先采用子間隔大的粗估計字典對時延值進行預估計,以減少計算量.粗估計字典ΨR[14,18]設計為N1×N1維:

式中ΨR為粗估計字典;ψri(u)為字典中第i個原子;N1為第一級輸入信號長度即字典原子個數;ρ為累積脈沖信號與第一級輸入信號的長度比,ρ=N/N1,其中N為累積脈沖信號長度,N11.

第二級由精估計字典提供對時延值的局部精確估計.精估計字典只包含部分相位的原子,其數量少、子間隔小,從而減小了字典尺寸、降低了計算量.精估計字典矩陣ΨP為N×D維,其數據量僅為傳統字典的D/N,即使原子個數遠小于信號長度,也可對時延值進行局部精確估計.精估計字典表示如下:

式中,ΨP為精估計字典;ψpj(u)為第j個原子;N為信號長度;D為字典原子個數,D

累積脈沖輪廓可由字典原子與稀疏系數表示為

式中,fR(u)是粗估計累積脈沖輪廓,ai是粗估計字典中第i個原子對應的稀疏系數,fP(u)是精估計累積脈沖輪廓,aj是精估計字典中第j個原子的稀疏系數.

3.2 測量矩陣選取

測量矩陣的選取對數據采樣和信號重構十分重要.由于信號在相同稀疏階次情況下,哈達瑪矩陣的恢復效果優于隨機0-1矩陣[16],本文選擇哈達瑪矩陣作為測量矩陣.測量矩陣采樣是一個降維的過程,第一級中可選取哈達瑪矩陣的前M1行,表示如下[15]:

式中ΦR是粗估計測量矩陣;H為哈達瑪矩陣,矩陣,用于選取哈達瑪矩陣的前M1行;yR是測量值;fR為第一級輸入信號;ΘR為感知矩陣;AR是一階稀疏系數矩陣,AR={a1,a2,a3,···,aN1?1}.

第二級精估計測量矩陣為

其中,ΦP是精估計測量矩陣;

為M2×N矩陣,其中M2=Sa2×N,Sa2為第二級信號采樣率;fP第二級輸入信號;ΘP為感知矩陣.

3.3 脈沖輪廓重構與時延估計算法

累積脈沖輪廓的重構是一階稀疏信號的優化問題,可選擇正交匹配追蹤法[19]進行累積脈沖輪廓重構.時延估計由粗估計與精估計兩級構成,圖1為兩級時延估計算法流程圖.

第一級粗估計與重構算法如下:

1)初始化殘差r0=yR,迭代選取的位置Λ1=φ,增量矩陣T1=φ,初始迭代次數i=1;

2)尋找感知矩陣中與殘差投影系數最大,即內積最大的列原子向量,保存最大投影系數對應的位置保存對應的列原子向量Ti=Ti?1∪θt;

3)求yR=Ti·ai的最小二乘解i;

4)輸出稀疏系數A,利用等式fR=ΨRA進行信號重構.

圖1 兩級時延估計算法流程圖Fig.1.Flow chart of two stage delay estimation algorithm.

第一級輸出索引值為t,即累脈沖輪廓相對于標準脈沖輪廓的延遲間隔數,可得預估時延值為τ1=t×ρ(μs),ρ為第二級與第一級輸入信號的比值.精確時延估計值可通過第二級得到.第二級精估計字典為N×D維,只包含標準脈沖輪廓的部分相位,因此,需利用預估時延值τ1的先驗條件,使累積脈沖輪廓f(u)處于精估計字典包含的相位范圍內,即將其移位到精估計字典中間位置D/2,移位值為?=(D/2)?τ1,移位修正后的輸入信號為fP(u)=f(u??),再進行第二級精確時延估計.

第二級精估計算法如下:

1)初始化殘差r0=yP,迭代選取的位置Λ1=φ,增量矩陣T1=φ,初始迭代次數i=1;

2)尋找感知矩陣中與殘差投影系數最大,即內積最大的列原子向量,保留最大投影系數對應的位置保存其對應的列原子向量Ti=Ti?1∪θt1.

上述步驟得出的t1為精估計索引值,最后可得到精確時延估計值τ=(t×ρ+t1?D/2)×1(μs).

3.4 計算復雜度分析

在壓縮感知中,波形字典大小對時延估計精度影響很大.

在信號長度為N,且采樣率為M/N的情況下,如使用一級壓縮感知方法,當輸入信號長度N=32768,字典原子間隔為1μs時,波形字典應為N×N維,即32768×32768.字典需占用較大內存,且運行時間長.

而采用本文方法,一級字典僅為N1×N1維,即1024×1024,二級字典為N×D維,即32768×300,二者所占存儲空間比一級字典小兩個數量級.

因此,本文方法能有效提高測量精度且兼顧合理的運行時間與內存.

4 仿真實驗與結果分析

使用歐洲脈沖星網絡數據庫(European Pulsar Network,EPN)數據[20]與硬X射線探測衛星(Rossi X-ray Timing Explorer,RXTE)實測數據[21]對本文方法進行了仿真實驗.EPN收集了超過1000顆脈沖星的脈沖輪廓;RXTE是美國NASA發射的具有高時間分辨率、高靈敏度、寬能譜范圍的大型硬X射線探測衛星,用于接收脈沖星光子以獲得實測數據.

利用EPN數據獲得標準脈沖輪廓,累積脈沖輪廓由泊松分布產生,以PSR B0531+21脈沖星作為導航星,對本文方法進行了輪廓重構、采樣率選擇、時延估計精度的仿真實驗.另外,通過RXTE獲取的脈沖星B0531+21的光子到達時刻進行不同時段的歷元折疊,以得到標準脈沖輪廓與累積脈沖輪廓,并在不同累積時間下與基于EPN數據的時延估計精度進行對比.

有關參數為:脈沖星B0531+21的周期P=0.0334 s,背景流量為λb=5×10?3ph/(s·cm2),來自脈沖星的源流量為λs=1.54 ph/(s·cm2),X射線探測器有效面積設為800 cm2.為驗證本文方法的實時性,將其在酷睿i5 CPU 2.1G計算機上進行了仿真.

4.1 基于EPN數據的仿真實驗

4.1.1 不同累積時間下的重構輪廓

為驗證本文方法在重構脈沖輪廓方面的有效性,將兩級壓縮感知與歷元折疊方法獲得的脈沖輪廓信噪比進行比較.

仿真實驗中,累積脈沖輪廓的子間隔數為N=32768,第一級粗估計字典為N1×N1維,測量矩陣為M1×N1維,N1=1024,M1=614,N1×N1,采樣率Sa1為M1/N1=0.6;第二級精估計字典為N2×D維,測量矩陣為M2×N2維,N2=N=32768,M2=1024,D=300,采樣率Sa2為M2/N2=0.03125.針對脈沖星B0531+21在不同累積時間下的脈沖輪廓的實驗結果如圖2所示.

圖2 不同累積時間下的重構輪廓 (a)累積時間1 s;(b)累積時間10 s;(c)累積時間100 sFig.2.Reconstructed pro filesat different integrated time:(a)Integrated time of 1 s;(b)integrated time of 10 s;(c)integrated time of 100 s.

圖2(a)—(c)分別給出了累積時間為1,10,100 s的歷元折疊法與本文方法提供的累積脈沖輪廓.由圖2可知,觀測時間較短時,壓縮感知恢復出的累積脈沖輪廓信噪比遠高于歷元折疊法的結果.因此,兩級壓縮感知方法能重構出高信噪比的累積脈沖輪廓.

4.1.2 采樣率的最優選擇

測量矩陣采樣點數與計算量以及存儲空間密切相關.當測量矩陣采樣點數少時,計算量較小,時延估計精度低;反之,計算量大,時延估計精度高.

針對脈沖星B0531+21在不同采樣率下的時延估計精度的仿真實驗中,第一級實驗參數與4.1.1節一致,第二級精估計字典尺寸也與4.1.1節相同,第二級測量矩陣大小為M2×N2維,其中,N2=N=32768,M2為采樣點數,分別有M2=512,1024,2048,4096,對應的采樣率Sa2為M2/N2=0.015625,0.03125,0.0625,0.125.

圖3(a)和圖3(b)分別給出了累積時間為10—100 s,500—2500 s兩個較短與較長時間內不同采樣率的估計誤差值,不同采樣率下的運行時間如表1所列.

圖3 不同時間與采樣率下的時延估計精度 (a)累積時間10—100 s;(b)累積時間500—2500 sFig.3.Time delay estimation accuracy at different time and sampling rates:(a)Integrated time of 10–100 s;(b)integrated time of 500–2500 s.

表1 不同采樣率下的運行時間Table 1.Running time with different sampling rates.

由圖3與表1可知:采樣率為0.03125對應的第二級采樣點數M2=1024,采樣率大于0.03125時,時延估計誤差緩慢減小,運行時間明顯增大,因此當累積時間大于40 s,選擇0.03125為最優的采樣率;但累積時間小于40 s時,不同采樣率下的估計誤差波動較大,因此當累積時間較短時,有必要增大采樣率,選擇采樣率為0.125.

4.1.3 一級壓縮感知與兩級方法的估計精度

圖4給出了脈沖星B0531+21在本文方法與一級壓縮感知方法下的估計精度比較結果.本文方法下字典與測量矩陣參數與實驗4.1.1節一致,一級壓縮感知方法下累積脈沖輪廓的子間隔數為N=1024,第一級字典為N×N維,測量矩陣為M×N維,M=614,采樣率Sa為M/N=0.6.由圖4可知,在時延估計精度上,基于兩級字典的壓縮感知方法優于一級壓縮感知.隨著累積時間的延長,時延估計精度顯著增加.

圖4 一級與兩級壓縮感知估計精度的比較Fig.4.Comparison of time delay estimation between onelevel compressed sensing and two-level compressed sensing.

4.2 基于RXTE實測數據的仿真實驗

為進一步驗證本文時延估計方法的性能,利用RXTE衛星收集到的實測數據進了仿真實驗.采用RXTE中Crab脈沖星觀測數據,數據包為40805-01-05-000[21],計算不同累積時間下的信號時延估計精度.

標準脈沖輪廓由數據包的光子到達時間序列經過歷元折疊獲得,截取數據包的一段光子到達時間序列數據,對其中所有數據加上10μs的時延,再利用新生成的時間序列進行歷元折疊,即可獲得具有確定時延的累積脈沖輪廓.不同累積時間的累積脈沖輪廓可通過截取數據包中不同時段的數據進行處理.

實驗中,第一級粗估計字典為N1×N1維,測量矩陣為M1×N1維,N1=1024,M1=614;第二級精估計字典為N2×D維,測量矩陣為M2×N2維,N2=N=32768,M2=512,其中D=300.

實驗結果列于表2,由表2可知,累積時間分別為10,30,80,100 s時,隨著累積時間的增加,RXTE實測數據得到的時延估計精度與EPN數據得到的時延估計精度都有所提高,二者的誤差值比較接近,這進一步表明了本文方法的有效性.

表2 RXTE實測數據與EPN數據的時延估計對比Table 2.Comparison of time delay estimation between RXTE real data and EPN data.

5 結 論

基于兩級壓縮感知的脈沖星累積輪廓時延估計是一種能保持高估計精度并降低計算量的有效方法.將字典分解為粗估計字典與精估計字典,可避免大數據量的計算,既能利用粗估計字典進行預估得到全局累積脈沖輪廓時延值,又可利用精估計字典對經過預估時延值調整的累積脈沖輪廓進行精確估計,從而實現了全局估計與局部估計的有機結合.實驗結果表明,該方法能保持對累積脈沖輪廓的高估計精度,并大幅降低計算量.

參考文獻

[1]Hanson J E 1996Ph.D.Dissertation(Stanford:Stanford University)

[2]Sheikh S I 2005Ph.D.Dissertation(Maryland:University of Maryland)

[3]Liu J,Wu J,Xiong L,Fang J C,Liu G 2017Chin.J.Electron6 1325

[4]Taylor J H 1992Philos.T.R.Soc.A341 117

[5]Tran N D,Renaux A,Boyer R,Marcos S 2014IEEE Trans.Aeros.Elec.Sys.50 786

[6]Kang Z W,He X,Liu J 2016Optik127 5050

[7]Zhang H,Xu L P,Xie Q,Luo N 2011Acta Phys.Sin.60 049701(in Chinese)[張華,許錄平,謝強,羅楠2011物理學報60 049701]

[8]Fang H Y,Liu B,Li X P,Sun H F,Xue M F,Shen L R,Zhu J P 2016Acta Phys.Sin.65 119701(in Chinese)[方海燕,劉兵,李小平,孫海峰,薛夢凡,沈利榮,朱金鵬2016物理學報65 119701]

[9]Liu J,Fang J C,Wu J,Kang Z W,Ning X L 2014IET Radar Sonar Navig8 1154

[10]Emadzadeh A A,Speyer J L 2011IEEE Trans.Aeros.Elec.Sys.47 2317

[11]Emadzadeh A A,Speyer J L 2010IEEE Trans.Sig.Proc.58 4484

[12]Li J X,Ke X Z 2010Acta Astronom.Sin.51 263(in Chinese)[李建勛,柯熙政 2010天文學報 51 263]

[13]Do T T,Gan L,Nguyen N H,Tran T D 2012IEEE Trans.Sig.Proc.60 139

[14]Su Z,Xu L P,Gan W 2011Sci.Sin.:Phys.Mech.Astron.41 681(in Chinese)[蘇哲,許錄平,甘偉 2011中國科學:物理學力學天文學41 681]

[15]Li S L,Liu K,Xiao L L 2014Optik125 1875

[16]Shen L R,Li X P,Sun H F,Fang H Y,Xue M F 2016Optik127 4379

[17]Golshan A R,Sheikh S I 2007Annual Meeting of Institute of NavigationCambrige,MA,USA,April 23–25,2007 p413

[18]Shi G,Lin J,Chen X Y,Qi F,Liu D H,Zhang L 2008IEEE Trans.Sig.Proc.55 379

[19]Tropp J A,Gilbert A C 2007IEEE Trans.Sig.Proc.53 4655

[20]EPN http://www jb man acuk/～pulsar/Resources/epn/browser html[2017-1-13]

[21]RXTE https://heasarcnasagov/docs/archivehtml[2017-5-24]