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

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

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

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

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

1 引 言

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

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

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

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

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

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

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

其均值和方差分別為

標(biāo)準(zhǔn)脈沖輪廓模型,即總流量密度函數(shù)λ(u),由X射線脈沖星的光子密度及來自背景環(huán)境的光子密度構(gòu)成:

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

對于同一顆脈沖星,累積脈沖輪廓與標(biāo)準(zhǔn)脈沖輪廓的關(guān)系可表示如下[2]:

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

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

i為累積脈沖輪廓相對于標(biāo)準(zhǔn)脈沖輪廓的延遲間隔數(shù),脈沖時延估計量為

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

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

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

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

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

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

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

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

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

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

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

3.2 測量矩陣選取

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

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

第二級精估計測量矩陣為

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

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

3.3 脈沖輪廓重構(gòu)與時延估計算法

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

第一級粗估計與重構(gòu)算法如下:

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

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

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

4)輸出稀疏系數(shù)A,利用等式fR=ΨRA進(jìn)行信號重構(gòu).

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

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

第二級精估計算法如下:

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

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

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

3.4 計算復(fù)雜度分析

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

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

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

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

4 仿真實驗與結(jié)果分析

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

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

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

4.1 基于EPN數(shù)據(jù)的仿真實驗

4.1.1 不同累積時間下的重構(gòu)輪廓

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

仿真實驗中,累積脈沖輪廓的子間隔數(shù)為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在不同累積時間下的脈沖輪廓的實驗結(jié)果如圖2所示.

圖2 不同累積時間下的重構(gòu)輪廓 (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可知,觀測時間較短時,壓縮感知恢復(fù)出的累積脈沖輪廓信噪比遠(yuǎn)高于歷元折疊法的結(jié)果.因此,兩級壓縮感知方法能重構(gòu)出高信噪比的累積脈沖輪廓.

4.1.2 采樣率的最優(yōu)選擇

測量矩陣采樣點(diǎn)數(shù)與計算量以及存儲空間密切相關(guān).當(dāng)測量矩陣采樣點(diǎn)數(shù)少時,計算量較小,時延估計精度低;反之,計算量大,時延估計精度高.

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

圖3(a)和圖3(b)分別給出了累積時間為10—100 s,500—2500 s兩個較短與較長時間內(nèi)不同采樣率的估計誤差值,不同采樣率下的運(yùn)行時間如表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 不同采樣率下的運(yùn)行時間Table 1.Running time with different sampling rates.

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

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

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

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

4.2 基于RXTE實測數(shù)據(jù)的仿真實驗

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

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

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

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

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

5 結(jié) 論

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

參考文獻(xiàn)

[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)[張華,許錄平,謝強(qiáng),羅楠2011物理學(xué)報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物理學(xué)報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天文學(xué)報 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中國科學(xué):物理學(xué)力學(xué)天文學(xué)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]