賈運澤,劉 勁,王奕迪,潘 超
(1.武漢科技大學(xué)信息科學(xué)與工程學(xué)院,武漢 430081;2.國防科技大學(xué)空天科學(xué)學(xué)院,長沙 410003;3.湖北經(jīng)濟大學(xué)信息與通信工程學(xué)院,武漢 430205)
脈沖星是一顆高穩(wěn)定度的中子星,而脈沖星導(dǎo)航[1-3]是一項不依賴于地面站的自主導(dǎo)航技術(shù)。同時,脈沖星能夠為近地空間、太陽系乃至星際旅行的航天器提供定軌、守時等功能,滿足航天任務(wù)在不同軌道下的持續(xù)高精度導(dǎo)航需求。脈沖星導(dǎo)航的基本原理為[4]:在脈沖星觀測周期內(nèi),航天器上安裝的X射線探測器會持續(xù)收集到脈沖星發(fā)出的X射線信號,并將這些信號在太陽系質(zhì)心(solar system barycenter,SSB)處折疊成一個固定的脈沖星周期[5-6]。這個過程被稱為歷元折疊(epoch folding,EF),經(jīng)過歷元折疊后會形成一個累積脈沖輪廓,將之與脈沖標(biāo)準(zhǔn)輪廓進行比較可以得到一個相位差,進而可以將相位差轉(zhuǎn)換為脈沖到達(dá)時間(time-of-arrival,TOA),將脈沖TOA 作為測量量即可求得航天器任意時刻的速度與位置信息。由此可以看出,脈沖星TOA 是脈沖星導(dǎo)航中的重要參數(shù),但是因為航天器的運動加上脈沖傳播過程中的噪聲干擾[7-8],會導(dǎo)致接收到的脈沖星信號周期發(fā)生變化,使積聚在固有周期內(nèi)的脈沖星累積輪廓發(fā)生畸變,從而引起脈沖星TOA 漂移[9]。因此,脈沖星周期的觀測對高精度導(dǎo)航有著重要的意義,當(dāng)前如何快速且高精度地獲取脈沖星周期已經(jīng)成為了研究的熱點。
脈沖星周期估計的方法可分為兩類:一類是利用脈沖星TOA 的漂移來估計周期誤差,另一類是利用脈沖星輪廓畸變來反演脈沖星周期,后者為當(dāng)前的主流。其基本原理為:嘗試按不同周期折疊脈沖星信號,得到一系列脈沖畸變輪廓,再找出畸變度最小的脈沖累積輪廓,其對應(yīng)的周期就是固有周期。目前最典型的就是時域χ2統(tǒng)計法[10]和傅立葉頻譜法[11]。除此之外,還有文獻(xiàn)[12]提出的一種基于頻率細(xì)分和連續(xù)倫姆周期圖(CLP)法,文獻(xiàn)[13]提出的基于TOA 信息的脈沖星周期估計法等等。以上方法都有一個共同的問題,那就是多次EF 帶來的高計算量,而對于深空計算設(shè)備來說,高計算量帶來的高負(fù)載是必須要避免的問題之一。
壓縮感知(compressive sensing,CS)[14-15]是一種新興的信號處理方法,它可以通過極少的觀測值重構(gòu)脈沖輪廓,從而解決脈沖星數(shù)據(jù)量過大這一問題。CS對于稀疏信號有較好的恢復(fù)能力,而脈沖信號就是一個稀疏信號,因此,將CS應(yīng)用到脈沖信號處理上已經(jīng)成為一個熱點話題。在CS過程中,測量矩陣的構(gòu)造和字典的設(shè)計是關(guān)鍵,如孫海峰等[16]將哈達(dá)瑪矩陣作為測量矩陣重構(gòu)脈沖星信號,吳春艷等[17]提出了基于粗估計與精估計兩級字典的時延估計法。劉勁等[18]將CS應(yīng)用在脈沖星周期估計中,提出一種基于FFT-CS的脈沖星周期快速估計方法,并將其與蝶形算法結(jié)合,減小了計算量并提高了精度。除此之外,CS 在其他方面也得到了應(yīng)用,如TOA 估計[19]、航天器定位[20]等。
為了進一步減少計算量并提高計算精度,本文提出了一種基于DCT-CS的脈沖星周期快速估計方法。通過構(gòu)造不同程度的畸變度字典來直接估計脈沖星周期。脈沖星信號的能量主要集中在低頻部分,將信號通過DCT變換為測量矩陣,然后經(jīng)過超分辨率稀疏恢復(fù)直接估計出X射線脈沖星周期。
本文提出的基于DCT-CS 的脈沖星周期超快速估計方法需要對脈沖星畸變輪廓進行DCT 計算來得到脈沖星周期。
在脈沖星導(dǎo)航的實際應(yīng)用中,由于航天器運動會造成多普勒效應(yīng),導(dǎo)致觀測信號周期改變。此時如果仍然使用提前預(yù)估的脈沖星周期而不是脈沖星的固有周期進行歷元折疊,就會導(dǎo)致累積脈沖輪廓產(chǎn)生輪廓畸變[21]。
離散余弦變換(discrete cosine transform,DCT)是從離散傅里葉變換(discrete Fourier transform,DFT)推導(dǎo)而來,是輸入函數(shù)為偶函數(shù)情況下的一種特殊的DFT,相當(dāng)于一個長度為其2倍的實偶序列的DFT。相比于DFT,DCT 主要有以下三個優(yōu)點:1)DCT 中不存在虛部,計算主要為實數(shù)變換,具有更好的計算效率。2)在不引入間斷的同時對信號施加了周期性。在DFT 中,當(dāng)時間信號被截斷并假定其周期性時,會在時域中引入不連續(xù)性,而DCT 即使假設(shè)信號存在對稱性,也不會在信號中引入不連續(xù)性和相關(guān)操作。3)DCT 比DFT 具有更好的能量聚集度,能夠?qū)⒛芰烤奂诘皖l部分后進行處理。一維DCT 的表達(dá)式如式(1)所示。
式中,f(x)為輸入信號;F(k)為輸出信號;k為輸出信號點位數(shù);Y為信號長度。
圖1為標(biāo)準(zhǔn)脈沖輪廓和畸變輪廓進行DCT 的結(jié)果對比。圖1(a)為標(biāo)準(zhǔn)輪廓進行DCT 的結(jié)果,圖1(b)為使用最大輪廓畸變度M=21的累積輪廓進行DCT 的結(jié)果,其中輪廓畸變度為畸變寬度與一個脈沖周期內(nèi)的周期間隔數(shù)N之比,表示輪廓畸變的程度。DCT 系數(shù)為脈沖輪廓信號經(jīng)DCT 后得到的數(shù)值,其基本特性為直流和低頻系數(shù)數(shù)值大,高頻系數(shù)數(shù)值小,幅值表示信號強度。由圖1可以看出,隨著輪廓的累積,信號強度會逐漸增大,同時信號能量會向周邊分散。因此,相比于直接對標(biāo)準(zhǔn)輪廓進行DCT 計算,對累積輪廓進行DCT 能更好地利用脈沖信號的低頻部分。
圖1 標(biāo)準(zhǔn)輪廓與畸變輪廓DCT結(jié)果Fig.1 DCT results of standard contour and distorted contour
提出了一種基于DCT-CS 的脈沖星周期快速估計算法。在傳統(tǒng)脈沖星周期估計方法使用的歷元折疊方法中,脈沖信號將會以不同的周期進行折疊并生成一系列的畸變輪廓,通過求解最小失真度來求解脈沖星周期。而本算法不需要進行嘗試性歷元折疊,該算法的核心思想為通過檢測脈沖星輪廓的畸變度來直接估計脈沖星周期誤差。
本算法包括以下四個部分:1)構(gòu)建畸變輪廓字典。2)設(shè)計低頻離散測量矩陣。3)獲取頻率偏移累積輪廓。4)使用最大元素超分辨率稀疏恢復(fù)。算法示意圖如圖2所示。其中,測量矩陣和字典的構(gòu)建是提前在地面完成的;脈沖星輪廓的積累和超分辨率稀疏恢復(fù)在航天器上進行。脈沖星輪廓的積累與脈沖星信號的采集是同步的。在收集完X射線脈沖星光子后,才會開始計算匹配矩陣和使用最大元素超分辨率稀疏恢復(fù)。該算法的目標(biāo)是在脈沖星信號采集后快速獲得脈沖星周期。因此,減少算法第四部分的計算負(fù)荷是主要任務(wù)。
圖2 基于DCT-CS的脈沖星周期超快速估計算法Fig.2 Ultra fast estimation algorithm for pulsar period based on DCT-CS
在本算法中也用到了歷元折疊,為了與傳統(tǒng)脈沖星周期估計方法中的歷元折疊區(qū)分,將本算法中的歷元折疊命名為基本歷元折疊,將傳統(tǒng)算法中的歷元折疊命名為嘗試性歷元折疊?;練v元折疊相比于嘗試性歷元折疊有以下區(qū)別:1)計算量不同,基本歷元折疊在整個脈沖星周期估計算法中只進行了一次,而傳統(tǒng)脈沖星周期估計方法中需要進行多次嘗試性歷元折疊。2)計算時間不同,基本歷元折疊可以在脈沖星信號的采集過程中同步進行,而嘗試性歷元折疊需要在脈沖星信號的采集之后進行,這大大節(jié)省了脈沖星周期估計的計算時間。
下面將對算法流程的四個部分進行逐個介紹。
(1)構(gòu)造畸變輪廓字典
不同頻率偏移引起的累積脈沖星輪廓的畸變程度不同,所以構(gòu)造畸變輪廓字典是必要的。
設(shè)h(φ)為歸一化的標(biāo)準(zhǔn)脈沖輪廓,其滿足以下條件
式中,φ為脈沖相位;h(φ)為一個周期為1的函數(shù),則有以下公式
畸變脈沖輪廓可視為標(biāo)準(zhǔn)脈沖輪廓與門函數(shù)的循環(huán)互相關(guān),因此,構(gòu)造了M個不同畸變度的脈沖輪廓φm,表達(dá)式為
式中,?循環(huán)相關(guān);m是畸變脈沖輪廓序號;^h為反向標(biāo)準(zhǔn)脈沖輪廓,表達(dá)式如式(5)所示
Gm(N×1)為傳播矢量,表達(dá)式如式(6)所示
從式(4)和式(6)可看出,Gm相當(dāng)于M個寬度不同的矩形窗。脈沖星輪廓畸變度定義為畸變脈沖輪廓序號數(shù)m與一個脈沖周期內(nèi)的周期間隔數(shù)N之比,即m/N。隨著m增大,矩形窗變寬,脈沖累積輪廓的畸變度變大。
將M個φm合成畸變輪廓合成一個畸變字典ψ(N×M),其可以表示為
式中,m={1,2,3…M-1};M為脈沖星最大輪廓畸變度,文中為21。
(2)設(shè)計低頻離散測量矩陣
因為脈沖星累積輪廓的相位是未知的,所以脈沖星周期估計方法必須不受相位影響。因為DCT系數(shù)與相位無關(guān),所以本文選用DCT 來求得測量矩陣,由于脈沖星輪廓的能量集中在低頻部分,所以取用低頻部分提取測量矩陣。
文中的低頻測量矩陣由字典中的一個元素φm經(jīng)DCT 得出,則低頻離散測量矩陣Φ的表達(dá)式為
假設(shè)對應(yīng)字典因素φm的頻率偏移為Δfm,可以表示為
式中,N為脈沖輪廓間隔數(shù)。
可以構(gòu)造感知矩陣Θ,即測量矩陣與輪廓字典的乘積,可表示為
(3)獲取頻率偏移累積輪廓
累積脈沖輪廓是根據(jù)脈沖固有周期將脈沖星輻射光子進行折疊而形成的,這一過程稱為歷元折疊。但是,多普勒效應(yīng)以及脈沖星周期躍變會使航天器接收到的脈沖周期存在誤差ΔT,則歷元折疊以T0+ΔT為周期累積脈沖輪廓,T0為預(yù)先設(shè)置的固有周期。脈沖星周期誤差與頻率偏移的關(guān)系為
式中,Δf是頻率偏移;f0是固有頻率。
如果不考慮相位,不能僅通過累積輪廓來區(qū)分Δf和-Δf,為了解決該問題,在此處引入脈沖星周期偏移量Toff?ΔT,-Toff+ΔT和-Toff-ΔT的符號是相同的,但是其振幅不同,由此可以區(qū)分出頻率的正負(fù)符號。在基本歷元折疊中,采用T0-Toff+ΔT為周期進行輪廓累積并進行如下操作:
收集一個觀測周期內(nèi)的所有光子信號,并將其按照[0,T0-Toff+ΔT]的脈沖星周期進行折疊,這就意味著其在一個觀測周期內(nèi)的計算模量為
然后,將周期持續(xù)時間分為等長的幾塊,計算每塊中的光子量為
最后,對得到的光子進行歸一化處理,得到累積脈沖輪廓。
(4)使用最大元素超分辨率稀疏恢復(fù)
設(shè)觀測矢量y的表達(dá)式為
匹配矢量的第m項zm,(m={0,1,2,…,M-1})可表示為觀測矢量和感知矢量的乘積
式中,am為感知矩陣Θ的第m列。
設(shè)匹配矩陣為S,其第m列Sm可表示為
式中,L為感知矩陣測量次數(shù),Θ(Μ,1∶L)代表感知矩陣第M行的1到L列。
然后得到S最大值的行下標(biāo)和列下標(biāo)
最后,脈沖星周期可估計為
本章中將計算整個算法的計算復(fù)雜度,算法流程如圖2所示。其中,脈沖星累積過程可在脈沖星信號的觀測時間內(nèi)進行,因此此處只計算提取低頻部分、匹配字典和超分辨率匹配估計的計算量,即乘-加操作次數(shù)(multiply-accumulate operations,MAC)。
1)提取低頻部分:在此部分對字典中的單個元素φm進行了DCT來提取信號低頻部分,所以計算量集中在DCT 部分,對于一個長度為Y的信號,一維DCT的計算量為Y2。設(shè)脈沖星累積輪廓間隔數(shù)為N,因此在此部分中計算量為N2MAC。
2)匹配字典:在此部分中進行了一次DCT 計算,如公式(10)所示,所以計算量依然為N2MAC。
3)超分辨率匹配估計:該部分的計算量主要由匹配矢量和匹配矩陣決定,其中匹配矢量如式(16)所示,共需要(M×L)MAC。匹配矩陣包含式(17)和式(18)兩方面,設(shè)匹配矩陣大小為L,故其中式(15)包含L次計算,式(18)包含L次計算,要計算出max(Sm)需要計算一次式(17)和M次式(18),共需(2L+M×L)MAC。
則算法整體需2N2+2L+2M×L次計算。在預(yù)先設(shè)定N=33 000,L=5 000,M=21的情況下,共需約2.2×109MAC。
為了對比體現(xiàn)本算法的優(yōu)勢,下面計算同樣情況下使用FFT 算法進行脈沖星周期估計所需的計算量。同樣計算提取信號低頻部分,匹配字典和超分辨率匹配估計三部分的計算量。在提取信號低頻部分中,計算量可近似看作FFT 部分的計算量,其計算量為N×log2NMAC。而在匹配字典部分,需要先進行FFT 再進行IFFT,所以需要進行2N×log2NMAC。最后的超分辨率匹配估計,其中匹配矢量共計算M×L次,由于FFT 需要進行二維計算,所以式(18)的計算量變?yōu)?L次,則計算出max(Sm)需要進行一次式(17),M次式(18)運算,共需要計算(3L+M×L)MAC。則采用FFT 進行脈沖星周期估計共需:N×log2N+2N×log2N+M×L+3L+M×L=3N×log2N+3L+3M×L,約為1.8×106MAC。
從上述計算量分析可知,基于DCT-CS的脈沖星周期估計方法在前期的低頻測量矩陣獲取與歷元折疊過程中的計算量較大,其原因在于DCT 的計算量與脈沖輪廓間隔數(shù)的平方,即N2成正比。當(dāng)已知航天器的速度和位置信息后,可減少搜索點的數(shù)量,屆時可大幅減少計算量。而在其后的脈沖星周期估計算法部分里,基于DCT-CS的脈沖星周期估計方法的計算量相較于FFT-CS較少。其原因在于DCT 沒有虛數(shù)部分且算法流程更為簡單,所以計算量較少。從結(jié)果來看,基于DCT-CS的脈沖星周期估計方法從計算量角度與基于FFT-CS的周期估計方法各有優(yōu)劣,需要通過實驗驗證分析。
本次仿真實驗中使用了歐洲脈沖星網(wǎng)絡(luò)(European pulsar network,EPN)的數(shù)據(jù),此數(shù)據(jù)庫收集了超過一千顆脈沖星的標(biāo)準(zhǔn)脈沖輪廓數(shù)據(jù),而且,數(shù)據(jù)都是采用ASCII字符存儲,使用者可以直接免費下載其純文本文件,非常方便。本實驗選取的脈沖星為蟹狀星云脈沖星PSR B0531+21,其標(biāo)準(zhǔn)脈沖輪廓圖如圖3所示,之所以選這顆脈沖星,是因為其具有較大的光子流量密度,便于觀測和實驗。
圖3 脈沖星PSR B0531+21標(biāo)準(zhǔn)輪廓Fig.3 Standard profile of pulsar PSR B0531+21
設(shè)時間分辨率為1μs。脈沖輪廓間隔數(shù)N等于脈沖星周期除以時間分辨率的商,其值為33 000。脈沖星輻射光子流量密度為1.54 ph/(cm2·s),背景噪聲光子流量密度是0.005 ph/(cm2·s)。X射線探測器的有效探測面積為1 m2,觀測時長為1 000 s,Toff為2.7×10-10s。測量矩陣大小為L×N,其中L為5 000,N為33 000。字典大小為M×L,其中M為21。實驗設(shè)備為一臺處理器為AMD7840@3.80 GHz,RAM 16 G的電腦。
本節(jié)研究觀測時長和探測器面積對脈沖星周期估計誤差的影響,同時通過對比分析不同觀測時長和探測器面積對結(jié)果的影響程度。
圖4 給出了不同觀測時長和探測器面積下的仿真結(jié)果,圖中橫縱坐標(biāo)軸為對數(shù)坐標(biāo)值。由圖4可以看出脈沖星周期估計誤差隨觀測時長和探測器面積增大而減少,且觀測時長對脈沖星周期估計誤差的影響強于探測器面積。由此可以得出,無論是面積的擴大還是觀測時間的延長,都有利于脈沖星周期精度的提高。而且,觀測時間的延長比觀測面積的增大更能提高周期估計的精度。
圖4 時間面積對周期估計誤差的影響Fig.4 The influence of time and area on period estimation error
其原因為根據(jù)脈沖星周期估計誤差的克拉美羅下界(Cramer-Rao lower bound,CRLB)推導(dǎo)得出,脈沖星周期估計誤差與探測器面積和觀測時長的立方成反比[18],其性質(zhì)可以表示為
式中,Te為脈沖星周期估計誤差;Tobs和A分別為觀測時長和探測器面積。由式(23)和式(24)可以看出提高觀測時長和增大探測器面積均可以提高脈沖星周期估計的精度,且提高觀測時長的影響比增大探測面積更大,而實驗結(jié)果與該性質(zhì)相符。
本節(jié)研究了測量次數(shù)對脈沖星周期估計誤差的影響,由于脈沖星輪廓具有較低的信噪比,且脈沖星輪廓的能量集中在低頻部分,因此本算法通過DCT 提取矩陣中選擇低頻部分,而不是傳統(tǒng)的隨機提取。
圖5表明了測量次數(shù)與脈沖星周期估計誤差的關(guān)系,圖中縱坐標(biāo)軸為對數(shù)值。可以看出,隨著測量次數(shù)的增加,脈沖星周期估計誤差逐漸減小,直至某一點到達(dá)最小值,對于不同的探測器面積,達(dá)到最小值的測量次數(shù)均為約5 000,大于5 000時誤差值趨近于定值。其原因為脈沖星輪廓的能量集中在低頻部分,因此小尺寸的測量矩陣就包含了絕大部分能量,后續(xù)再增加測量次數(shù)能得到的能量很小且會引入噪聲。因此本次實驗中采用L=5 000作為測量矩陣的測量次數(shù)。
圖5 測量次數(shù)對周期估計誤差的影響Fig.5 The influence of observation number on period estimation error
光通量和背景噪聲是脈沖星周期估計中的重要參數(shù),本節(jié)研究了光通量和背景噪聲對脈沖星周期估計誤差的影響。
圖6給出了3 000次蒙特卡羅實驗中不同光通量和背景噪聲下脈沖星周期估計的誤差結(jié)果,圖中橫、縱坐標(biāo)軸為對數(shù)坐標(biāo)值。由圖6可以看出,脈沖星周期估計誤差隨著脈沖星光子通量的增加而減小,而隨著背景噪聲的增加而增大,脈沖星周期估計誤差隨著脈沖星光子通量的增加而減小,而隨著噪聲級的增加而增大。因此,無論是增大脈沖星光子通量還是減小背景噪聲,都有利于提高脈沖星周期估計精度。
圖6 光通量和噪聲對周期估計誤差的影響Fig.6 Impacts of flux and noise on period estimation error
本算法中使用DCT 替換了傳統(tǒng)的FFT 來進行脈沖星周期估計,現(xiàn)對比這兩種算法的優(yōu)劣。
表1給出了不同觀測時長和探測器面積的情況下,采用FFT 和DCT 進行脈沖星周期估計的誤差結(jié)果和進行歷元折疊需要的計算時間。由表1可以看出,DCT 的計算時間少于FFT 且與參數(shù)無關(guān),這是因為本算法的計算時間與光子量無關(guān)。而在周期估計誤差方面,在觀測時長較短和探測器面積較小的情況下,采用FFT 的誤差結(jié)果會遠(yuǎn)大于采用DCT 得到的結(jié)果,但隨著觀測時長的增加和探測器面積的增大,兩者的估計誤差結(jié)果會逐漸接近,最終采用FFT 得到的結(jié)果會優(yōu)于采用DCT 方法得到的結(jié)果。其原因在于FFT 較DCT 多出來一個虛數(shù)部分,計算精度損失低于DCT,具有較高的魯棒性,在高分辨率的情況下,達(dá)到的結(jié)果會好于DCT。
表1 兩種周期估計方法的對比Tab.1 Comparison of two cycle estimation methods
從理論角度出發(fā),DCT 對比FFT 有以下幾個優(yōu)勢:1)DCT 相比于FFT 有更好的能量聚集度,能更好地處理脈沖信號。2)DCT 沒有虛部,相比于FFT 擁有更小的計算復(fù)雜度。從實際角度出發(fā),考慮到航天器高速飛行過程中產(chǎn)生的多普勒效應(yīng)的影響,需盡快完成脈沖星周期估計的計算過程,才能保證算法的實時性。根據(jù)表1的實驗結(jié)果表明,基于DCT 的脈沖星周期估計方法需要的計算時間更短。而基于FFT 的周期估計方法只有在高分辨率的情況下的誤差估計結(jié)果才會優(yōu)于基于DCT 的脈沖星周期估計方法,限制條件過高。因此本文采用DCT 進行脈沖星周期估計。
本文提出了一種基于DCT-CS 的脈沖星周期超快速估計算法,該算法不需要進行嘗試性歷元折疊,而是通過低頻測量矩陣、輪廓畸變字典和超分辨率稀疏恢復(fù)算法來直接估計脈沖星周期。
本算法有三個優(yōu)點:1)計算負(fù)荷低。DCT 對低頻信號有著優(yōu)秀處理能力,該算法舍棄了常用的FFT方法,而是通過DCT 來進行低頻矩陣的提取和構(gòu)建測量矩陣。理論分析表明基于DCT的脈沖星周期估計方法計算負(fù)荷要低于基于FFT的脈沖星周期估計方法。2)計算精度高。文中通過實驗評估了基于DCT的脈沖星周期估計方法的計算精度,實驗結(jié)果表明其計算精度達(dá)到了3.82×10-12s,優(yōu)于傳統(tǒng)脈沖星周期估計方法。3)計算時間短。由于基于DCT的脈沖星周期估計方法的計算量較低,所以相應(yīng)的計算時間也較短。實驗結(jié)果表明基于DCT的脈沖星周期估計方法的計算時間僅為9.31 ms,要少于基于FFT的脈沖星周期估計方法的14.9 ms。
綜上所述,基于DCT-CS的脈沖星周期快速估計算法可在降低系統(tǒng)計算負(fù)擔(dān)的同時提高周期估計精度。