褚梓樾, 王 雪, 趙 航, 商 鵬
(1. 中國科學院國家授時中心, 陜西 西安 710600; 2. 中國科學院大學, 北京 101408;3. 中國科學院精密導航定位與定時技術重點實驗室, 陜西 西安 710600)
隨著衛(wèi)星通信、遙感、導航等空間應用需求的不斷提升,空間在軌衛(wèi)星數目急劇增加。為了更好地了解和掌握空間軌道環(huán)境,避免衛(wèi)星碰撞等災難性事故的發(fā)生,需要借助輻射源識別技術對非合作衛(wèi)星進行偵測和識別[1-3]。輻射源識別是通過發(fā)射信號表現出的差異,從復雜的電磁環(huán)境中區(qū)分出不同目標的。由于輻射源硬件設備固有的不理想特性,發(fā)射信號的調制特征參數會產生一定的偏差[4-7](如:正交相移鍵控信號的正交調制偏差),可將調制特征參數作為細微特征,用于輻射源識別。在軌衛(wèi)星的正常運行以及與地面站之間的常態(tài)信息交互離不開測控系統[8],測控信號作為日常工作信號,其細微特征可作為衛(wèi)星輻射源識別的穩(wěn)定特征參數。目前在衛(wèi)星測控任務中多采用統一S頻段(unified S band, USB)測控信號[9-10]。調制指數作為USB測控信號的一個重要參數,決定了信號的功率分配關系[11],即使發(fā)射相同信號的輻射源之間也存在功率分配的細微差異。因此,本文將調制指數作為非合作USB測控信號輻射源的一個細微特征,對其進行精確估計。
關于信號調制指數估計的研究,已有的一些方法[12-15]主要針對連續(xù)相位調制信號,該類方法均是在連續(xù)相位的基礎之上建立的,而USB測控信號由于測距碼和副載波調制信息的存在,導致信號相位并不連續(xù)。其中,文獻[13]利用數據輔助的方法來提升算法性能,本文的研究是基于非合作通信,需要對信號進行盲分析,因此這些方法并不適用。文獻[16-17]針對角度調制信號的調制指數提出了一種盲估計方法,主要是通過對信號先進行微分再取希爾伯特變換,得到相應的解析形式,對其取模,再根據幅度差異得到調制指數估計值,其研究對象僅具有一個調制指數且基帶調制信號形式簡單,而USB測控信號具有多個調制指數且基帶調制信號往往是多個信號相加復合的形式,利用相應估計算法無法得到各個調制指數的估計值。
因此,本文針對非合作USB測控信號的調制指數,提出了一種基于頻域譜分析的估計方法,該方法在低信噪比條件下仍具有較高的估計精度。
USB測控體制上/下行信號載頻比為221/240,上行信號的載頻范圍為2 025~2 110 MHz,主要用于傳輸遙控信號和測距信號。下行信號的載頻范圍為2 200~2 290 MHz,主要用于傳輸遙測信號和回傳的測距信號。本文利用高增益天線,在地面接收衛(wèi)星下行信號,因此主要對下行信號進行介紹。
根據文獻[10-11]可知,USB測控下行信號的表達式一般為
s(t)=Acos[ωct+β1x1(t)+β2x2(t)]
(1)
式中,A為信號幅度;ωc為載波角頻率;x1(t)和x2(t)分別為遙測信號和測距信號;β1和β2分別為遙測信號和測距信號的調制指數。
遙測信號的構成方式一般為脈沖編碼調制/相移鍵控,信號表達式如下:
(2)
測距信號的構成方式一般為偽隨機噪聲(pseudo-random noise, PRN)序列或單音信號。
(1) 當測距信號為PRN序列時,信號表達式如下:
(3)
(2) 當測距信號為單音信號時,信號表達式如下:
x2(t)=sin(2πfrt)
(4)
式中,fr為單音信號頻率。
本文所研究的非合作USB測控信號表達式為
s(t)=Acos[ωct+bc0(t)+mc1(t)sin(Ωt)]
(5)
式中,c0(t)為PRN序列構成的測距信號;c1(t)sin(Ωt)為脈沖編碼調制/二進制相移鍵控方式構成的遙測信號,c1(t)為遙測信息碼,Ω為副載波角頻率;b和m分別對應測距信號和遙測信號的調制指數。相應的仿真信號頻譜如圖1所示。通過圖1可以看出,該信號頻譜特性明顯,具有較高的辨識度。其中,標簽“1”“8”表示ωc±2Ω處的副載波分量,標簽“2”“7”表示調制在ωc±Ω處c0(t)c1(t)的頻譜,標簽“3”“6”表示調制在ωc±Ω處c1(t)的頻譜,標簽“4”表示調制在ωc處c0(t)的頻譜,標簽“5”表示ωc處的載波分量。
圖1 仿真信號頻譜Fig.1 Simulation signal spectrum
為了便于理論推導,將式(5)轉化為相應的解析信號形式:
s(t)=Aej[ωct+b c0(t)+m c1(t)sin(Ω t)]
(6)
式中,c0(t)和c1(t)為+1或-1,結合貝塞爾函數與三角函數之間的轉化的關系:
(7)
(8)
同時根據貝塞爾函數的對稱性質J-m(x)=(-1)m·Jm(x),最終s(t)可整理為
s(t)=A[cosb+jc0(t)sinb]·
(9)
對s(t)進行傅里葉變換,可得
F[s(t)]=A[2πδ(ω)cosb+jF[c0(t)]sinb]*
(10)
式中,F[·]表示傅里葉變換;δ(·)表示單位沖激函數;*表示卷積運算。
J2k-1(m)2πδ(ω-(ωc+(2k-1)Ω))*F[c1(t)]]*
(11)
由于圖1所示的信號頻譜存在明顯的沖激譜線特征(標簽“1”“5”“8”),所以考慮將沖激譜線的幅度與調制指數之間建立聯系。分別取k=0和k=1(或k=-1),可得
(12)
(13)
(14)
將m值代入式(12)中的第一個等式,可得
(15)
為了方便表示,記
則有
(16)
若將m值代入式(12)中的第二個等式,可得
(17)
同理有
(18)
式(16)和式(18)是等價的。
上述求解調制指數過程中所涉及到的非線性函數JP(m)與f(b)對應的函數曲線如圖2所示。
圖2 非線性函數曲線Fig.2 Nonlinear function curve
該算法基于對信號頻譜特征的分析結果,主要的計算量花費在對信號作傅里葉變換,因此該算法的復雜度為O(NlogN),N為信號長度。為了便于更好地理解調制指數的估計過程,圖3給出了相應的估計流程。
圖3 USB測控信號調制指數估計算法流程圖Fig.3 Flow chart of modulation index estimation algorithm for USB telemetry tracking and command signal
根據上述理論推導過程以及流程圖可知,調制指數估計需要獲知一些先驗參數,包括:載波角頻率ωc,副載波角頻率Ω,信號幅度A,c0(t)的碼速率和碼相位。這些先驗參數可分別根據文獻[18-21]中的方法得到其精確估計值。
同時,本文算法具有一定的普適性,其核心思想是建立信號頻域解析表達式與信號頻譜特性之間的聯系,構造參數估計方程,通過對方程進行求解,進而得到相應參數的估計值。針對一些其他類型非合作信號的調制指數估計問題(如:文獻[16-17]中涉及到的當基帶調制信號為單音信號或二進制相移鍵控信號時的相位調制信號),可對本文信號模型作相應簡化與修改,實現具體參數的求解。
本文提出的調制指數估計方法主要受頻譜譜線幅度和信號幅度影響,因此根據文獻[22]中的誤差傳遞公式,通過參數估計式(14)對M0和M2及式(16)對M0和A求偏導數以獲取估計值的均方根誤差,具體表達式如下:
(19)
(20)
式中,
由文獻[23]可知:
(21)
(22)
(23)
(24)
(25)
式中,
(26)
(27)
其中,E(·)為均值;D(·)為方差。根據高階矩[24]公式:
(28)
式中,σ為噪聲標準差,則
(29)
故有
(30)
進而可得
(31)
從式(22)和式(23)可以明顯看出,調制指數估計值的均方根誤差與信噪比呈負相關,即隨著信噪比的增大,均方根誤差隨之減小。
根據式(5)利用蒙特卡羅方法在高斯白噪聲條件下進行仿真,信號相關仿真參數設置如下:信號幅度A=1,載波頻率fc=2 210 MHz,副載波頻率f1=4 kHz,調制指數m=0.6,b=0.8,碼速率fc0=2 MHz,fc1=1 kHz,碼元初始相位均為0 s,采樣頻率fs=5 MHz,采樣時間為1 s。信噪比變化范圍為[-10,10]dB,變化步長為2 dB。在每個信噪比條件下仿真Ns=2 000次,若將某信噪比條件下調制指數m的第i次估計值記為mSNR,i,那么在該信噪比條件下估計值的均方根誤差可表示為
(32)
式中,m為真實的調制指數。
同理,對于調制指數b有
(33)
式中,b為真實的調制指數。
根據式(32)和式(33)分別給出調制指數m與調制指數b的仿真估計誤差隨信噪比變化曲線,同時利用在第3.2節(jié)得到理論估計誤差式(22)和式(23),給出理論估計誤差隨信噪比變化曲線,并與仿真曲線進行對比。
為了分析調制指數取值對估計誤差的影響,利用控制變量法的思想在原有仿真參數的基礎上,僅對m和b做相應修改,進行以下兩組對照試驗:① 保持b=0.8不變,取m=0.4;② 保持m=0.6不變,取b=0.6。將以上兩組對照實驗得到的估計誤差曲線與初始設定參數(m=0.6,b=0.8)對應的估計誤差曲線進行對比。以上仿真實驗對比結果如圖4所示。
圖4 調制指數估計誤差隨信噪比變化關系Fig.4 Relationship of modulation index estimation error varies with signal to noise ratio
從圖4中可以看出,隨著信噪比的增加,調制指數的估計誤差也隨之降低,同時理論估計誤差曲線與仿真估計誤差曲線也基本一致。通過對比兩個調制指數誤差曲線可以發(fā)現,相同參數條件下在信噪比為-10 dB時,調制指數m的誤差小于調制指數b,這是由于m的值僅由信號頻譜譜線幅度決定,而b的值除譜線幅度以外還受到信號幅度的影響,信噪比對信號幅度的影響遠大于對頻譜譜線幅度的影響,當信噪比逐漸增大時,b的估計誤差下降趨勢明顯快于m,當信噪比達到10 dB時,b的估計誤差小于m的估計誤差。
為便于更清晰地描述對照組實驗結果,表1給出了在SNR=10 dB條件下的量化數據對比結果。
表1 估計誤差對比結果(SNR=10 dB)
結合圖4中的曲線關系與表1中的量化數值可知,當b取值一定時,隨著m取值的減小,m的估計誤差σm隨之增大,而b的估計誤差σb則基本保持不變;當m取值一定時,隨著b取值的減小,m的估計誤差σm也隨之減小,而b的估計誤差σb則隨之增大。根據式(22)和式(23)可知,σm同時受m和b取值的影響,而σb僅受b取值的影響。同時結合圖2中的JP(m)和f(b)函數曲線可知,m與JP(m)呈反比例關系,b與f(b)呈反比例關系,而σm與JP(m)呈正比關系與f(b)呈反比關系,σb與f(b)呈正比關系。因此,σm與m呈反比例關系,與b呈正比例關系;σb與b呈反比例關系。上述仿真結果與理論分析結果保持一致。
本文采用的實測信號是在中國科學院國家授時中心西安場區(qū),利用13 m高增益天線采集到的某顆低軌衛(wèi)星發(fā)出的USB測控信號,實測信號數據采樣頻率為250 MHz,采樣時長為2 s。取3組等長且不同時段的實測信號數據,按照圖3所示的流程得到調制指數的估計值,同時其他參數的估計結果如表2所示。
表2 實測信號參數估計值
為了驗證本文提出的調制指數估計方法的準確性,利用上述參數估計值生成恢復信號,并與相應的實測信號從頻域和時域兩方面進行對比。觀察以上3組參數估計值,除了由時間段選取所導致的碼相位不同以外,其余參數估計值的差異較小,在此選用第1組估計值進行對比驗證。
3.2.1 頻域對比
根據USB測控信號頻域式(10)可知,信號頻譜幅度受調制指數的影響,因此通過對比實測信號與恢復信號頻譜幅度的差異可以判斷調制指數估計算法的準確性。相應的信號頻譜對比結果如圖5所示。
圖5 實測信號與恢復信號頻譜對比Fig.5 Spectrum comparison between measured signal and recovered signal
為了更好地展示實測信號與恢復信號的頻譜對比結果,下面將標簽a~標簽e進行適當放大,此外還將標簽b和標簽d進行了平滑處理,結果如圖6所示。
圖6 經放大和平滑處理后標簽a~標簽e對比結果Fig.6 Comparison results of labels a~e after amplification and smoothing
根據式(12)可知,標簽a,標簽c,標簽e的頻譜幅度主要由調制指數決定,結合上述對比結果,可以看出恢復信號和實測信號的頻譜幅度基本一致。其中,標簽a、標簽b、標簽d對應的實測信號頻譜幅度大于恢復信號頻譜幅度,而標簽c和標簽e對應的實測信號頻譜幅度小于恢復信號頻譜幅度,并且經過平滑處理后的標簽b和標簽d中恢復信號頻譜與實測信號頻譜的波動性存在些許差別。這些現象產生的原因是,在利用高增益天線采集衛(wèi)星信號的過程中,采集到的信號中不僅包含噪聲,還包含許多干擾信號,而采集設備的采樣率有限,會導致一定的頻譜混疊現象,將一些干擾信號的頻譜混疊進來,這些干擾信號會對實測信號的頻譜幅值產生影響,而恢復信號中不包含干擾信號,因此不會產生頻譜混疊的問題。通過計算相應點幅度差,可以得到標簽c對應的頻譜幅值誤差小于0.01 dB,標簽a、標簽b、標簽d和標簽e對應的頻譜幅值誤差小于0.08 dB,表明利用本文提出的估計算法得到的調制指數較為準確。
3.2.2 時域對比
若求整體時域基帶信號,通過式(6)可知,調制指數m和b包含在基帶信號的相位中,無法直觀地進行對比,利用幅度受調制指數決定的基帶信號進行對比則更為直觀。結合USB測控信號時域式(9)與圖5所示的信號頻譜圖,記圖5中標簽c對應的基帶信號為sbase1,b對應的基帶信號為sbase2,標簽a對應的基帶信號為sbase3。以上3個基帶信號的近似表達式為
(34)
式(34)直觀地體現出調制指數對基帶信號幅度的影響。標簽d對應基帶信號與標簽b的僅相差一個負號;標簽e對應基帶信號與標簽a的相同。此外由于標簽a和標簽e的能量較低,通過解調無法得到較好的基帶波形,且sbase1和sbase3呈倍數關系,因此只需對比標簽b和標簽c對應的實測與恢復基帶信號幅度,即可驗證調制指數估計的準確性。相應的實測與恢復基帶信號對比如圖7所示。
圖7 實測信號與恢復信號的基帶信號對比Fig.7 Comparison of baseband signal between measured signal and recovered signal
根據基帶信號的對比結果以及式(34),sbase1的I路在理想情況下應為常數,但卻出現波動現象,并且Q路在碼片持續(xù)時間內也具有波動性。sbase2的I、Q兩路也存在類似情況,這均是由于信號頻譜相鄰成分的旁瓣導致的。sbase1對應的實測信號與恢復信號一致性較好,sbase2對應的實測信號與恢復信號波動性存在些許差異。這是因為標簽c信號頻譜能量較高,不易受噪聲干擾,且實測信號頻譜處基本不存在干擾信號,而標簽b相對標簽c能量較低,且實測信號頻譜存在干擾信號。總體上,sbase1和sbase2對應的實測信號與恢復信號幅度基本一致,即驗證了調制指數估計的準確性。
為了量化時域對比結果,引入歸一化互相關函數[25]的概念。通過計算實測基帶信號與恢復基帶信號的歸一化互相關函數最大值,若該值越接近于1,說明實測基帶信號與恢復基帶信號的一致性越高。歸一化互相關函數計算公式如下:
(35)
式中,sreal(t)表示實測基帶信號;sreco(t)表示恢復基帶信號;(·)*表示共軛。
利用式(35)分別計算sbase1和sbase2對應實測信號與恢復信號的歸一化互相關函數最大值,在表3中列出了相應的計算結果。
表3 基帶信號歸一化互相關函數最大值
上述結果表明,sbase1和sbase2對應的實測信號與恢復信號的歸一化互相關函數最大值均接近于1,即實測基帶信號與恢復基帶信號一致性較高,能夠說明調制指數估計得較為準確。
本文針對非合作USB測控信號的調制指數,提出了基于頻域譜分析的估計方法。首先根據信號頻譜的譜線特性,由信號頻域表達式推導出調制指數估計方法,并對該方法的估計誤差進行了理論分析。通過蒙特卡羅仿真得到調制指數估計誤差隨信噪比變化關系曲線,同時將理論估計誤差曲線與仿真估計誤差曲線進行對比,驗證了理論分析的正確性。并仿真分析了在不同調制指數取值的條件下,調制指數估計誤差的變化情況。在此基礎上,利用實測低軌衛(wèi)星信號數據對調制指數估計方法進行了驗證,根據估計參數生成恢復信號,從頻域和時域兩個方面與實測信號進行對比。通過信號頻譜幅度、時域基帶信號波形及其歸一化互相關函數最大值的對比結果表明,實測信號與恢復信號的頻譜幅度誤差較小、時域基帶信號波形的一致性較高、歸一化互相關函數最大值接近于1。進一步說明本文提出的調制指數估計方法精度較高,可為非合作衛(wèi)星輻射源識別提供高精度的參數支撐。