基于DCT-CS的脈沖星周期超快速估計方法

賈運澤,劉 勁,王奕迪,潘 超

(1.武漢科技大學信息科學與工程學院,武漢 430081;2.國防科技大學空天科學學院,長沙 410003;3.湖北經濟大學信息與通信工程學院,武漢 430205)

0 引言

脈沖星是一顆高穩定度的中子星,而脈沖星導航[1-3]是一項不依賴于地面站的自主導航技術。同時,脈沖星能夠為近地空間、太陽系乃至星際旅行的航天器提供定軌、守時等功能,滿足航天任務在不同軌道下的持續高精度導航需求。脈沖星導航的基本原理為[4]:在脈沖星觀測周期內,航天器上安裝的X射線探測器會持續收集到脈沖星發出的X射線信號,并將這些信號在太陽系質心(solar system barycenter,SSB)處折疊成一個固定的脈沖星周期[5-6]。這個過程被稱為歷元折疊(epoch folding,EF),經過歷元折疊后會形成一個累積脈沖輪廓,將之與脈沖標準輪廓進行比較可以得到一個相位差,進而可以將相位差轉換為脈沖到達時間(time-of-arrival,TOA),將脈沖TOA 作為測量量即可求得航天器任意時刻的速度與位置信息。由此可以看出,脈沖星TOA 是脈沖星導航中的重要參數,但是因為航天器的運動加上脈沖傳播過程中的噪聲干擾[7-8],會導致接收到的脈沖星信號周期發生變化,使積聚在固有周期內的脈沖星累積輪廓發生畸變,從而引起脈沖星TOA 漂移[9]。因此,脈沖星周期的觀測對高精度導航有著重要的意義,當前如何快速且高精度地獲取脈沖星周期已經成為了研究的熱點。

脈沖星周期估計的方法可分為兩類:一類是利用脈沖星TOA 的漂移來估計周期誤差,另一類是利用脈沖星輪廓畸變來反演脈沖星周期,后者為當前的主流。其基本原理為:嘗試按不同周期折疊脈沖星信號,得到一系列脈沖畸變輪廓,再找出畸變度最小的脈沖累積輪廓,其對應的周期就是固有周期。目前最典型的就是時域χ2統計法[10]和傅立葉頻譜法[11]。除此之外,還有文獻[12]提出的一種基于頻率細分和連續倫姆周期圖(CLP)法,文獻[13]提出的基于TOA 信息的脈沖星周期估計法等等。以上方法都有一個共同的問題,那就是多次EF 帶來的高計算量,而對于深空計算設備來說,高計算量帶來的高負載是必須要避免的問題之一。

壓縮感知(compressive sensing,CS)[14-15]是一種新興的信號處理方法,它可以通過極少的觀測值重構脈沖輪廓,從而解決脈沖星數據量過大這一問題。CS對于稀疏信號有較好的恢復能力,而脈沖信號就是一個稀疏信號,因此,將CS應用到脈沖信號處理上已經成為一個熱點話題。在CS過程中,測量矩陣的構造和字典的設計是關鍵,如孫海峰等[16]將哈達瑪矩陣作為測量矩陣重構脈沖星信號,吳春艷等[17]提出了基于粗估計與精估計兩級字典的時延估計法。劉勁等[18]將CS應用在脈沖星周期估計中,提出一種基于FFT-CS的脈沖星周期快速估計方法,并將其與蝶形算法結合,減小了計算量并提高了精度。除此之外,CS 在其他方面也得到了應用,如TOA 估計[19]、航天器定位[20]等。

為了進一步減少計算量并提高計算精度,本文提出了一種基于DCT-CS的脈沖星周期快速估計方法。通過構造不同程度的畸變度字典來直接估計脈沖星周期。脈沖星信號的能量主要集中在低頻部分,將信號通過DCT變換為測量矩陣,然后經過超分辨率稀疏恢復直接估計出X射線脈沖星周期。

1 輪廓畸變對DCT的影響

本文提出的基于DCT-CS 的脈沖星周期超快速估計方法需要對脈沖星畸變輪廓進行DCT 計算來得到脈沖星周期。

在脈沖星導航的實際應用中,由于航天器運動會造成多普勒效應,導致觀測信號周期改變。此時如果仍然使用提前預估的脈沖星周期而不是脈沖星的固有周期進行歷元折疊,就會導致累積脈沖輪廓產生輪廓畸變[21]。

離散余弦變換(discrete cosine transform,DCT)是從離散傅里葉變換(discrete Fourier transform,DFT)推導而來,是輸入函數為偶函數情況下的一種特殊的DFT,相當于一個長度為其2倍的實偶序列的DFT。相比于DFT,DCT 主要有以下三個優點:1)DCT 中不存在虛部,計算主要為實數變換,具有更好的計算效率。2)在不引入間斷的同時對信號施加了周期性。在DFT 中,當時間信號被截斷并假定其周期性時,會在時域中引入不連續性,而DCT 即使假設信號存在對稱性,也不會在信號中引入不連續性和相關操作。3)DCT 比DFT 具有更好的能量聚集度,能夠將能量聚集在低頻部分后進行處理。一維DCT 的表達式如式(1)所示。

式中,f(x)為輸入信號;F(k)為輸出信號;k為輸出信號點位數;Y為信號長度。

圖1為標準脈沖輪廓和畸變輪廓進行DCT 的結果對比。圖1(a)為標準輪廓進行DCT 的結果,圖1(b)為使用最大輪廓畸變度M=21的累積輪廓進行DCT 的結果,其中輪廓畸變度為畸變寬度與一個脈沖周期內的周期間隔數N之比,表示輪廓畸變的程度。DCT 系數為脈沖輪廓信號經DCT 后得到的數值,其基本特性為直流和低頻系數數值大,高頻系數數值小,幅值表示信號強度。由圖1可以看出,隨著輪廓的累積,信號強度會逐漸增大,同時信號能量會向周邊分散。因此,相比于直接對標準輪廓進行DCT 計算,對累積輪廓進行DCT 能更好地利用脈沖信號的低頻部分。

圖1 標準輪廓與畸變輪廓DCT結果Fig.1 DCT results of standard contour and distorted contour

2 基于DCT-CS的脈沖星周期超快速估計

提出了一種基于DCT-CS 的脈沖星周期快速估計算法。在傳統脈沖星周期估計方法使用的歷元折疊方法中,脈沖信號將會以不同的周期進行折疊并生成一系列的畸變輪廓,通過求解最小失真度來求解脈沖星周期。而本算法不需要進行嘗試性歷元折疊,該算法的核心思想為通過檢測脈沖星輪廓的畸變度來直接估計脈沖星周期誤差。

本算法包括以下四個部分:1)構建畸變輪廓字典。2)設計低頻離散測量矩陣。3)獲取頻率偏移累積輪廓。4)使用最大元素超分辨率稀疏恢復。算法示意圖如圖2所示。其中,測量矩陣和字典的構建是提前在地面完成的;脈沖星輪廓的積累和超分辨率稀疏恢復在航天器上進行。脈沖星輪廓的積累與脈沖星信號的采集是同步的。在收集完X射線脈沖星光子后,才會開始計算匹配矩陣和使用最大元素超分辨率稀疏恢復。該算法的目標是在脈沖星信號采集后快速獲得脈沖星周期。因此,減少算法第四部分的計算負荷是主要任務。

圖2 基于DCT-CS的脈沖星周期超快速估計算法Fig.2 Ultra fast estimation algorithm for pulsar period based on DCT-CS

在本算法中也用到了歷元折疊,為了與傳統脈沖星周期估計方法中的歷元折疊區分,將本算法中的歷元折疊命名為基本歷元折疊,將傳統算法中的歷元折疊命名為嘗試性歷元折疊?；練v元折疊相比于嘗試性歷元折疊有以下區別:1)計算量不同,基本歷元折疊在整個脈沖星周期估計算法中只進行了一次,而傳統脈沖星周期估計方法中需要進行多次嘗試性歷元折疊。2)計算時間不同,基本歷元折疊可以在脈沖星信號的采集過程中同步進行,而嘗試性歷元折疊需要在脈沖星信號的采集之后進行,這大大節省了脈沖星周期估計的計算時間。

下面將對算法流程的四個部分進行逐個介紹。

(1)構造畸變輪廓字典

不同頻率偏移引起的累積脈沖星輪廓的畸變程度不同,所以構造畸變輪廓字典是必要的。

設h(φ)為歸一化的標準脈沖輪廓,其滿足以下條件

式中,φ為脈沖相位;h(φ)為一個周期為1的函數,則有以下公式

畸變脈沖輪廓可視為標準脈沖輪廓與門函數的循環互相關,因此,構造了M個不同畸變度的脈沖輪廓φm,表達式為

式中,?循環相關;m是畸變脈沖輪廓序號;^h為反向標準脈沖輪廓,表達式如式(5)所示

Gm(N×1)為傳播矢量,表達式如式(6)所示

從式(4)和式(6)可看出,Gm相當于M個寬度不同的矩形窗。脈沖星輪廓畸變度定義為畸變脈沖輪廓序號數m與一個脈沖周期內的周期間隔數N之比,即m/N。隨著m增大,矩形窗變寬,脈沖累積輪廓的畸變度變大。

將M個φm合成畸變輪廓合成一個畸變字典ψ(N×M),其可以表示為

式中,m={1,2,3…M-1};M為脈沖星最大輪廓畸變度,文中為21。

(2)設計低頻離散測量矩陣

因為脈沖星累積輪廓的相位是未知的,所以脈沖星周期估計方法必須不受相位影響。因為DCT系數與相位無關,所以本文選用DCT 來求得測量矩陣,由于脈沖星輪廓的能量集中在低頻部分,所以取用低頻部分提取測量矩陣。

文中的低頻測量矩陣由字典中的一個元素φm經DCT 得出,則低頻離散測量矩陣Φ的表達式為

假設對應字典因素φm的頻率偏移為Δfm,可以表示為

式中,N為脈沖輪廓間隔數。

可以構造感知矩陣Θ,即測量矩陣與輪廓字典的乘積,可表示為

(3)獲取頻率偏移累積輪廓

累積脈沖輪廓是根據脈沖固有周期將脈沖星輻射光子進行折疊而形成的,這一過程稱為歷元折疊。但是,多普勒效應以及脈沖星周期躍變會使航天器接收到的脈沖周期存在誤差ΔT,則歷元折疊以T0+ΔT為周期累積脈沖輪廓,T0為預先設置的固有周期。脈沖星周期誤差與頻率偏移的關系為

式中,Δf是頻率偏移;f0是固有頻率。

如果不考慮相位,不能僅通過累積輪廓來區分Δf和-Δf,為了解決該問題,在此處引入脈沖星周期偏移量Toff?ΔT,-Toff+ΔT和-Toff-ΔT的符號是相同的,但是其振幅不同,由此可以區分出頻率的正負符號。在基本歷元折疊中,采用T0-Toff+ΔT為周期進行輪廓累積并進行如下操作:

收集一個觀測周期內的所有光子信號,并將其按照[0,T0-Toff+ΔT]的脈沖星周期進行折疊,這就意味著其在一個觀測周期內的計算模量為

然后,將周期持續時間分為等長的幾塊,計算每塊中的光子量為

最后,對得到的光子進行歸一化處理,得到累積脈沖輪廓。

(4)使用最大元素超分辨率稀疏恢復

設觀測矢量y的表達式為

匹配矢量的第m項zm,(m={0,1,2,…,M-1})可表示為觀測矢量和感知矢量的乘積

式中,am為感知矩陣Θ的第m列。

設匹配矩陣為S,其第m列Sm可表示為

式中,L為感知矩陣測量次數,Θ(Μ,1∶L)代表感知矩陣第M行的1到L列。

然后得到S最大值的行下標和列下標

最后,脈沖星周期可估計為

3 計算復雜度分析

本章中將計算整個算法的計算復雜度,算法流程如圖2所示。其中,脈沖星累積過程可在脈沖星信號的觀測時間內進行,因此此處只計算提取低頻部分、匹配字典和超分辨率匹配估計的計算量,即乘-加操作次數(multiply-accumulate operations,MAC)。

1)提取低頻部分:在此部分對字典中的單個元素φm進行了DCT來提取信號低頻部分,所以計算量集中在DCT 部分,對于一個長度為Y的信號,一維DCT的計算量為Y2。設脈沖星累積輪廓間隔數為N,因此在此部分中計算量為N2MAC。

2)匹配字典:在此部分中進行了一次DCT 計算,如公式(10)所示,所以計算量依然為N2MAC。

3)超分辨率匹配估計:該部分的計算量主要由匹配矢量和匹配矩陣決定,其中匹配矢量如式(16)所示,共需要(M×L)MAC。匹配矩陣包含式(17)和式(18)兩方面,設匹配矩陣大小為L,故其中式(15)包含L次計算,式(18)包含L次計算,要計算出max(Sm)需要計算一次式(17)和M次式(18),共需(2L+M×L)MAC。

則算法整體需2N2+2L+2M×L次計算。在預先設定N=33 000,L=5 000,M=21的情況下,共需約2.2×109MAC。

為了對比體現本算法的優勢,下面計算同樣情況下使用FFT 算法進行脈沖星周期估計所需的計算量。同樣計算提取信號低頻部分,匹配字典和超分辨率匹配估計三部分的計算量。在提取信號低頻部分中,計算量可近似看作FFT 部分的計算量,其計算量為N×log2NMAC。而在匹配字典部分,需要先進行FFT 再進行IFFT,所以需要進行2N×log2NMAC。最后的超分辨率匹配估計,其中匹配矢量共計算M×L次,由于FFT 需要進行二維計算,所以式(18)的計算量變為2L次,則計算出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 的計算量與脈沖輪廓間隔數的平方,即N2成正比。當已知航天器的速度和位置信息后,可減少搜索點的數量,屆時可大幅減少計算量。而在其后的脈沖星周期估計算法部分里,基于DCT-CS的脈沖星周期估計方法的計算量相較于FFT-CS較少。其原因在于DCT 沒有虛數部分且算法流程更為簡單,所以計算量較少。從結果來看,基于DCT-CS的脈沖星周期估計方法從計算量角度與基于FFT-CS的周期估計方法各有優劣,需要通過實驗驗證分析。

4 仿真結果與討論

本次仿真實驗中使用了歐洲脈沖星網絡(European pulsar network,EPN)的數據,此數據庫收集了超過一千顆脈沖星的標準脈沖輪廓數據,而且,數據都是采用ASCII字符存儲,使用者可以直接免費下載其純文本文件,非常方便。本實驗選取的脈沖星為蟹狀星云脈沖星PSR B0531+21,其標準脈沖輪廓圖如圖3所示,之所以選這顆脈沖星,是因為其具有較大的光子流量密度,便于觀測和實驗。

圖3 脈沖星PSR B0531+21標準輪廓Fig.3 Standard profile of pulsar PSR B0531+21

設時間分辨率為1μs。脈沖輪廓間隔數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。實驗設備為一臺處理器為AMD7840@3.80 GHz,RAM 16 G的電腦。

4.1 觀測時長和探測器面積的影響

本節研究觀測時長和探測器面積對脈沖星周期估計誤差的影響,同時通過對比分析不同觀測時長和探測器面積對結果的影響程度。

圖4 給出了不同觀測時長和探測器面積下的仿真結果,圖中橫縱坐標軸為對數坐標值。由圖4可以看出脈沖星周期估計誤差隨觀測時長和探測器面積增大而減少,且觀測時長對脈沖星周期估計誤差的影響強于探測器面積。由此可以得出,無論是面積的擴大還是觀測時間的延長,都有利于脈沖星周期精度的提高。而且,觀測時間的延長比觀測面積的增大更能提高周期估計的精度。

圖4 時間面積對周期估計誤差的影響Fig.4 The influence of time and area on period estimation error

其原因為根據脈沖星周期估計誤差的克拉美羅下界(Cramer-Rao lower bound,CRLB)推導得出,脈沖星周期估計誤差與探測器面積和觀測時長的立方成反比[18],其性質可以表示為

式中,Te為脈沖星周期估計誤差;Tobs和A分別為觀測時長和探測器面積。由式(23)和式(24)可以看出提高觀測時長和增大探測器面積均可以提高脈沖星周期估計的精度,且提高觀測時長的影響比增大探測面積更大,而實驗結果與該性質相符。

4.2 測量次數的影響

本節研究了測量次數對脈沖星周期估計誤差的影響,由于脈沖星輪廓具有較低的信噪比,且脈沖星輪廓的能量集中在低頻部分,因此本算法通過DCT 提取矩陣中選擇低頻部分,而不是傳統的隨機提取。

圖5表明了測量次數與脈沖星周期估計誤差的關系,圖中縱坐標軸為對數值?？梢钥闯?隨著測量次數的增加,脈沖星周期估計誤差逐漸減小,直至某一點到達最小值,對于不同的探測器面積,達到最小值的測量次數均為約5 000,大于5 000時誤差值趨近于定值。其原因為脈沖星輪廓的能量集中在低頻部分,因此小尺寸的測量矩陣就包含了絕大部分能量,后續再增加測量次數能得到的能量很小且會引入噪聲。因此本次實驗中采用L=5 000作為測量矩陣的測量次數。

圖5 測量次數對周期估計誤差的影響Fig.5 The influence of observation number on period estimation error

4.3 光通量與背景噪聲的影響

光通量和背景噪聲是脈沖星周期估計中的重要參數,本節研究了光通量和背景噪聲對脈沖星周期估計誤差的影響。

圖6給出了3 000次蒙特卡羅實驗中不同光通量和背景噪聲下脈沖星周期估計的誤差結果,圖中橫、縱坐標軸為對數坐標值。由圖6可以看出,脈沖星周期估計誤差隨著脈沖星光子通量的增加而減小,而隨著背景噪聲的增加而增大,脈沖星周期估計誤差隨著脈沖星光子通量的增加而減小,而隨著噪聲級的增加而增大。因此,無論是增大脈沖星光子通量還是減小背景噪聲,都有利于提高脈沖星周期估計精度。

圖6 光通量和噪聲對周期估計誤差的影響Fig.6 Impacts of flux and noise on period estimation error

4.4 DCT與FFT的對比

本算法中使用DCT 替換了傳統的FFT 來進行脈沖星周期估計,現對比這兩種算法的優劣。

表1給出了不同觀測時長和探測器面積的情況下,采用FFT 和DCT 進行脈沖星周期估計的誤差結果和進行歷元折疊需要的計算時間。由表1可以看出,DCT 的計算時間少于FFT 且與參數無關,這是因為本算法的計算時間與光子量無關。而在周期估計誤差方面,在觀測時長較短和探測器面積較小的情況下,采用FFT 的誤差結果會遠大于采用DCT 得到的結果,但隨著觀測時長的增加和探測器面積的增大,兩者的估計誤差結果會逐漸接近,最終采用FFT 得到的結果會優于采用DCT 方法得到的結果。其原因在于FFT 較DCT 多出來一個虛數部分,計算精度損失低于DCT,具有較高的魯棒性,在高分辨率的情況下,達到的結果會好于DCT。

表1 兩種周期估計方法的對比Tab.1 Comparison of two cycle estimation methods

從理論角度出發,DCT 對比FFT 有以下幾個優勢:1)DCT 相比于FFT 有更好的能量聚集度,能更好地處理脈沖信號。2)DCT 沒有虛部,相比于FFT 擁有更小的計算復雜度。從實際角度出發,考慮到航天器高速飛行過程中產生的多普勒效應的影響,需盡快完成脈沖星周期估計的計算過程,才能保證算法的實時性。根據表1的實驗結果表明,基于DCT 的脈沖星周期估計方法需要的計算時間更短。而基于FFT 的周期估計方法只有在高分辨率的情況下的誤差估計結果才會優于基于DCT 的脈沖星周期估計方法,限制條件過高。因此本文采用DCT 進行脈沖星周期估計。

5 結論

本文提出了一種基于DCT-CS 的脈沖星周期超快速估計算法,該算法不需要進行嘗試性歷元折疊,而是通過低頻測量矩陣、輪廓畸變字典和超分辨率稀疏恢復算法來直接估計脈沖星周期。

本算法有三個優點:1)計算負荷低。DCT 對低頻信號有著優秀處理能力,該算法舍棄了常用的FFT方法,而是通過DCT 來進行低頻矩陣的提取和構建測量矩陣。理論分析表明基于DCT的脈沖星周期估計方法計算負荷要低于基于FFT的脈沖星周期估計方法。2)計算精度高。文中通過實驗評估了基于DCT的脈沖星周期估計方法的計算精度,實驗結果表明其計算精度達到了3.82×10-12s,優于傳統脈沖星周期估計方法。3)計算時間短。由于基于DCT的脈沖星周期估計方法的計算量較低,所以相應的計算時間也較短。實驗結果表明基于DCT的脈沖星周期估計方法的計算時間僅為9.31 ms,要少于基于FFT的脈沖星周期估計方法的14.9 ms。

綜上所述,基于DCT-CS的脈沖星周期快速估計算法可在降低系統計算負擔的同時提高周期估計精度。