一種聯(lián)合估計信號參數(shù)和位置的高精度被動定位方法
【專利摘要】本發(fā)明提供一種聯(lián)合估計信號參數(shù)和位置的高精度被動定位方法，在DPD?unknown的基礎(chǔ)上，通過聯(lián)合估計發(fā)射信號的參數(shù)和位置來提高定位性能。基站對接收到的觀測信號進行短時傅里葉變換估計出接收信號的瞬時頻率，再利用網(wǎng)格搜索的方法，聯(lián)合估計發(fā)射信號的波形參數(shù)和目標(biāo)的位置。本發(fā)明充分挖掘發(fā)射信號的波形信息，針對線性調(diào)頻信號的時頻特征，結(jié)合短時傅里葉變換聯(lián)合估計信號參數(shù)和位置，使定位性能在信號參數(shù)的幫助下有極大提高。由于采用了時頻分析工具短時傅里葉變換，使該方法的計算量處于一個合理水平，計算代價小。
【專利說明】
一種聯(lián)合估計信號參數(shù)和位置的高精度被動定位方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明屬于雷達信號處理技術(shù)，特別涉及多基地雷達被動定位和低信噪比下的參 數(shù)估計技術(shù)。
【背景技術(shù)】
[0002] 日益激烈的電子戰(zhàn)要求在復(fù)雜環(huán)境下仍能對微弱目標(biāo)進行高精度的定位，因此如 何提高定位精度受到了越來越多的關(guān)注。同時，隨著芯片技術(shù)的飛速發(fā)展，處理器的數(shù)據(jù)處 理能力得到了極大改善，大數(shù)據(jù)的處理成為了可能。
[0003] 被動定位技術(shù)由其諸多優(yōu)點如隱蔽性、低功耗等，受到了越來越多的關(guān)注。被動定 位通常采用幾個分布較開的雷達基站截獲目標(biāo)發(fā)射機的發(fā)射信號，相比主動雷達，由于雷 達自身不需要發(fā)射源發(fā)射電磁波，因此被動定位也被稱為無源定位。目前，正處在研制當(dāng)中 的無源雷達有兩種:一種是靠目標(biāo)輻射器探測和定位目標(biāo)，而另一種則利用其它輻射器發(fā) 射的電磁波。對前一種而言，目標(biāo)本身就是輻射源，它們可能是機載雷達、干擾機和通信設(shè) 備;對于后者來說，目標(biāo)本身不直接發(fā)射電磁波，而是當(dāng)有非合作輻射器發(fā)射的電磁波通過 目標(biāo)時，利用從目標(biāo)反射的電磁波來獲取目標(biāo)數(shù)據(jù)。
[0004] 對于第一種情況，也可稱之為發(fā)射機定位，較為常用的方法是時差定位TD0A、到達 角定位Α0Α。它們的基本思路是，首先在各基站接收機中提取出有關(guān)目標(biāo)發(fā)射機的到達時間 或到達角度，再將這些參數(shù)傳輸給處理中心，通過求解方程式的方法估計出目標(biāo)位置。這種 方法雖然傳輸?shù)臄?shù)據(jù)量小，計算簡單，但是其定位精度低，無法滿足高精度定位的需要。
[0005] 文南犬"Direct Position Determination of Narrowband Radio Frequency Transmitters，IEEE Signal Process·Lett·，vol · 11，no·5，pp·513-516，May 2004"提出了 一種聯(lián)合處理各基站接收機觀測數(shù)據(jù)的直接定位方法DPD，該方法沒有傳統(tǒng)定位方法的參 數(shù)提取過程，盡可能地保留了目標(biāo)信息，仿真顯示有明顯的優(yōu)勢。這種方法考慮了兩種情 況，一種是目標(biāo)發(fā)射機發(fā)射的信號已知，如訓(xùn)練信號或同步信號，我們稱這種情況的DPD定 位算法為DPD-known算法。另一種情況是目標(biāo)發(fā)射的信號完全未知，主要針對非合作的發(fā)射 機，在電子戰(zhàn)中最為常見。由于不知道目標(biāo)信號的形式，DPD定位方法通過使目標(biāo)函數(shù)特征 值最大的方法實現(xiàn)了發(fā)射機定位，仿真顯示其定位性能仍然優(yōu)于傳統(tǒng)的定位方法，稱這種 方法為DPD-unknown算法。但是，該方法忽略了發(fā)射信號的波形信息，其定位精度受限，不能 適應(yīng)微弱目標(biāo)的高精度定位。

【發(fā)明內(nèi)容】

[0006] 本發(fā)明所要解決的技術(shù)問題是，針對微弱非合作目標(biāo)發(fā)射機提供一種更高精度的 定位方法。
[0007] 本發(fā)明為解決上述技術(shù)問題所采用的技術(shù)方案是，一種聯(lián)合估計信號參數(shù)和位置 的高精度被動定位方法，包括以下步驟：
[0008] 步驟1、初始化系統(tǒng)參數(shù):初始化接收機個數(shù)、各接收機位置、采樣間隔、各接收機 通道上的噪聲協(xié)方差矩陣，短時傅里葉變換的窗長度和步長；
[0009] 步驟2、從雷達各接收機中讀取量測，對量測進行采樣，得到L個離散的量測向量 ri,l = l,2,...,L；
[0010] 步驟3、根據(jù)預(yù)先設(shè)置的短時傅里葉變換窗長度和步長，對各接收機量測向量^作 短時傅里葉變換得到短時傅里葉變換的功率譜％ ;
[0011] 步驟4、根據(jù)短時傅里葉變換功率譜&，在每個頻率上尋找最大功率譜以及最大功 率譜對應(yīng)的頻率作為估計的接收信號的瞬時頻率/；和最大功率
[0012] 步驟5、劃分目標(biāo)位置網(wǎng)格搜索區(qū)間（Xgrid，ygrid)，信號發(fā)射時間的網(wǎng)格搜索區(qū)間 tgrid，以及發(fā)射的信號長度的網(wǎng)格搜索區(qū)間Tgrid;
[0013] 步驟6、對每個網(wǎng)格點rugrid=(叉811〇1^11〇1山11〇1，1 811(〇計算相應(yīng)的信號起始時間1：1,3 (nigrid)和終止時間tl,e(nigrid)，其中，Xgrid,ygrid表示網(wǎng)格點的起始坐標(biāo)，tgrid表示信號發(fā)射 時間的網(wǎng)格搜索區(qū)間，Tgrld表示信號長度的網(wǎng)格搜索區(qū)間；
[0014] 步驟7、根據(jù)起始時間和終止時間對應(yīng)的時間區(qū)間，從瞬時頻率/；中取出對應(yīng)區(qū)間 的頻率值，作為線性調(diào)頻信號的瞬時頻率的估計^"^；)，從最大功率中取出對應(yīng)區(qū) 間的最大功率值);
[0015] 步驟8、對從L個量測中估計出的L個線性調(diào)頻信號的瞬時頻率的估計取 平均，得到線性調(diào)頻信號的頻率估計值對L個量測中取出的最大功率 取平均再求和，得到區(qū)間功率和ω (rUgrid);
[0016] 步驟9、通過最小二乘參數(shù)擬合從線性調(diào)頻信號的頻率估計值估計線 性調(diào)頻信號的初始頻率/〇 ^丨咖）和斜率挪?作為估計的線性調(diào)頻信號參數(shù)
[0017]步驟10、計算最小二乘擬合誤差的區(qū)間平均值l(nlgrid)，將最小二乘擬合誤差的區(qū) 間平均值C(TUgrid)與區(qū)間功率和ω (TUgrid)相除得到質(zhì)量因子e(rugrid);
[0018] 步驟11、根據(jù)網(wǎng)格點nigrid = (xgrid，ygrid，tgrid，Tgrid)和估計的線性調(diào)頻信號參數(shù) 計算網(wǎng)格點nigrid的加權(quán)似然函數(shù)值G，卜;
[0019] 其中，計算網(wǎng)格點nigrid的加權(quán)似然函數(shù)值：
[0020]
[0021 ]其中，Η表不對矩陣求共輒轉(zhuǎn)置，e(rugrid)為質(zhì)量因子，ri為量測向量，Ri代表第1個 接收機通道上目標(biāo)附近噪聲及干擾的協(xié)方差矩陣Λ為根據(jù)信號的起始時間t^JrUgnd)、信 號長度Tgrid、和計估計的線性調(diào)頻信號參數(shù)承構(gòu)建的匹配函數(shù)在時間上離散化后，得 到的離散化后的匹配函數(shù)；
[0022]步驟12、重復(fù)步驟6至11，直到遍歷數(shù)據(jù)平面上所有網(wǎng)格點，得到全數(shù)據(jù)平面的加 權(quán)似然函數(shù)值，并找出取得最大似然函數(shù)值對應(yīng)的網(wǎng)格點的位置作為估計的發(fā)射機所在位 置。
[0023] 本發(fā)明在DPD-unknown的基礎(chǔ)上，通過聯(lián)合估計發(fā)射信號的參數(shù)和位置來提高定 位性能。針對線性調(diào)頻發(fā)射機(指該發(fā)射機發(fā)射一個線性調(diào)頻信號），首先對基站接收到的 觀測信號進行短時傅里葉變換（一種與傅里葉變換相關(guān)的數(shù)學(xué)變換，估計方法見文獻 "Accurate parameters estimation of chirp signal in low SNR,Audio,Language and Image Process. ,2014 International Conference,pp.551-555,July 2014·''）估計出接 收信號的瞬時頻率，再利用網(wǎng)格搜索的方法，聯(lián)合估計發(fā)射信號的波形參數(shù)和目標(biāo)的位置。 本發(fā)明相對于現(xiàn)有DPD-unknown算法，充分挖掘發(fā)射信號的波形信息，針對線性調(diào)頻信號的 時頻特征，結(jié)合短時傅里葉變換聯(lián)合估計信號參數(shù)和位置，使定位性能在信號參數(shù)的幫助 下有極大提高。由于采用了時頻分析工具短時傅里葉變換，使該方法的計算量處于一個合 理水平，在實際工程中具有一定的實用性。
[0024] 本發(fā)明的有益效果是，有效現(xiàn)了對微弱發(fā)射機的高精度定位，且計算代價小，可以 直接應(yīng)用到現(xiàn)有被動雷達的高精度定位中。
【附圖說明】
[0025] 圖1為本發(fā)明提供方法的流程框圖；
[0026] 圖2為對某一接收機信號的短時傅里葉變換功率譜平面；
[0027] 圖3為從某一接收機信號估計的瞬時頻率；
[0028] 圖4為從某一接收機信號估計的最大功率譜；
[0029 ]圖5為本發(fā)明方法與現(xiàn)有DPD-unknown方法、Dro-known方法(假設(shè)發(fā)射的信號參數(shù) 完全已知)定位性能的對比效果圖（其中DPD-STFT為本發(fā)明方法）。
【具體實施方式】
[0030]本發(fā)明主要采用計算機仿真的方法進行驗證，所有步驟、結(jié)論都在MATLAB-R2014a 上驗證正確。具體實施步驟如圖1所示：
[0031]步驟1、初始化系統(tǒng)參數(shù)：
[0032] 初始化系統(tǒng)參數(shù)包括：目標(biāo)發(fā)射機的位置，雷達基站接收機的個數(shù)L = 4,雷達基站 接收機位置。初始化采樣周期Ts= l(T7s，短時傅里葉變換的窗長度400TS，步長60TS，噪聲協(xié) 方差矩陣Ri。
[0033] 步驟2、產(chǎn)生用于仿真的雷達接收機量測向 =1，2,. . .，L。其中Ns表示量測樣本個數(shù)。這里用于仿真的調(diào)頻信號初始頻率為10MHz，調(diào)頻 斜率為0.1 MHz/ys。
[0034] 步驟3、根據(jù)設(shè)置的窗長度和步長，對Γ1作短時傅里葉變換，得到其功率譜 P,.; = |Fr, |2。這里&為Μ X N的矩陣，通過MATLAB-R2014a中命令spectrogram可直接得到，如圖 2〇
[0035] 步驟4、對Pr,，在每一列η上尋找對應(yīng)最大值的行標(biāo)m，對應(yīng)于第η列的時間1和第m行 的頻率f(U)，估計出接收信號的瞬時頻率/；，如圖3,其中/；為行向量，第η列 /(/") = argmax巧'同時計算每一列!!上的最大值) = Pr(，得到向量，如圖4。
[0036] 步驟5、劃分目標(biāo)位置，信號發(fā)射時間以及信號長度的網(wǎng)格搜索區(qū)間。
[0037] 步驟6、對某個網(wǎng)格點TUgrid - ( Xgrid，ygrid，tgrid，Tgrid ) , Xgrid，ygrid表不網(wǎng)格點的目標(biāo)
位置坐標(biāo)，tgrid表不信號發(fā)射的時間的網(wǎng)格搜索區(qū)間坐標(biāo)，Tgrid表不發(fā)射的信號的長度的網(wǎng) 格搜索區(qū)間坐標(biāo);計算相應(yīng)的信號起始時間 和終止 時間 tl, e(rUgrid ) - tl, s (Hlgrid )+Tgrid。其中，C為;/^巡，XI，y 1乃國込·)"安Η乂奶生個不。
[0038] 步驟7、根據(jù)起始時間和終止時間對應(yīng)的時間區(qū)間，從瞬時頻率/；·,中取出對應(yīng)時間 區(qū)間的頻率值，作為線性調(diào)頻信號的瞬時頻率的估計同時，根據(jù)時間區(qū)間在最大 功率1中取出對應(yīng)時間區(qū)間的最大功率值Plftn, 。
[0039] 步驟8、對從L個量測中估計出的L個線性調(diào)頻信號的瞬時頻率的估計 平均得到估計。同時對L個量測中取出的最大功率1^^取平均再求和得到區(qū)間功 率和 ω (rugrid)。
[0040] 步驟9、采用最小二乘參數(shù)擬合，通過參數(shù)擬合從又,；)中估計線性調(diào)頻信號的 初始頻率又)和斜率。在MATLAB-R2014a中可直接調(diào)用擬合函數(shù)ρο 1 yf i t，估計 出初始頻率和調(diào)頻斜率，同時得到擬合誤差Gmin(nigrid)。
[0041] 步驟10、計算最小二乘擬合誤差的區(qū)間平均值(MMCF)
得其與 區(qū)間功率和ω (qlgrid)相除得到質(zhì)量因虧
-\-ngna i
[0042] 步驟11、根據(jù)網(wǎng)格點nigrid = (xgrid，ygrid，tgrid，Tgrid)和估計的線性調(diào)頻信號參數(shù)
，計算加權(quán)似然函數(shù)值G ?冰D。其中，計算加權(quán)似 然函數(shù)值的步驟為：
[0043] 1:根據(jù)信號的起始時間t^sU^nd)，信號長度Tgrid，以及估計出的信號參數(shù) 來構(gòu)
[0044] 建匹配函數(shù)夏,:
[0045]
其中， t表示當(dāng)前時間變量，1 = 1，2,...，L，T Τ表示接收機的觀測時間長度這里觀測時間長度要 求盡量長，以保證接收到完整的發(fā)射信號為根據(jù)信號參數(shù)陜復(fù)出 的發(fā)射信號，分號后的4 ^ ~ 表不該發(fā)射信號的信號參數(shù)；
[0046]
[0047] 2 :將匹配函數(shù)在時間上離散化，得到離散化后的匹配函數(shù) 二[元[0]忑[：1]，,..4[%-1]]。
[0048] 3:計算對應(yīng)網(wǎng)格點nigrid的加權(quán)似然函數(shù)值：
[0049]
[0050] 其中，Η表示對矩陣求共輒轉(zhuǎn)置，h代表第1個接收機通道上目標(biāo)附近噪聲及干擾 的協(xié)方差矩陣。
[0051] 步驟12、重復(fù)步驟6至11，直到遍歷數(shù)據(jù)平面上所有網(wǎng)格點，得到全數(shù)據(jù)平面的加 權(quán)似然函數(shù)值，并找出取得最大似然函數(shù)值對應(yīng)的網(wǎng)格點作為估計值奪，同時得 到信號參數(shù)的估計值⑷)和Α承），從而估計出發(fā)射機位置(又4。
[0052]在給定仿真參數(shù)下，對每一個信噪比（SNR)進行500次蒙特卡洛仿真實驗，以均方 根誤差(RMSE)為定位性能依據(jù)，對比DPD-unknown算法，發(fā)現(xiàn)其定位性能有明顯的優(yōu)勢。同 時，以DPD-known算法作為定位誤差的下界對比本發(fā)明方法在估計信號參數(shù)上的性能損失。 如圖5,在高信噪比情況下（SNR大于10dB)，三個定位算法的定位誤差都較低且相同；當(dāng)信噪 比低于1 OdB時，本發(fā)明提出的方法(DFO-STFT)明顯優(yōu)于DPD-unknown算法，且定位性能在一 定范圍內(nèi)（SNR大于-8dB小于10dB)與理論最優(yōu)的定位算法(DPD-known，假設(shè)發(fā)射的信號參 數(shù)完全已知)相接近。因此，在發(fā)射信號參數(shù)未知的情況下，本發(fā)明提出的算法相比原有的 定位算法(DFO-unknown)在低信噪比下有明顯的優(yōu)勢。
【主權(quán)項】
1. 一種聯(lián)合估計信號參數(shù)和位置的高精度被動定位方法，其特征在于，包括W下步驟： 步驟1、初始化系統(tǒng)參數(shù):初始化接收機個數(shù)L、各接收機位置、采樣間隔、各接收機通道 上的噪聲協(xié)方差矩陣，短時傅里葉變換的窗長度和步長； 步驟2、從雷達各接收機中讀取量測，對量測進行采樣，得到L個離散的量測向量ri，1 = 1,2,...,L; 步驟3、根據(jù)預(yù)先設(shè)置的短時傅里葉變換窗長度和步長，對各接收機量測向量ri作短時 傅里葉變換得到短時傅里葉變換的功率譜Ρι,; 步驟4、根據(jù)短時傅里葉變換功率譜Ρ。，在每個頻率上尋找最大功率譜W及最大功率譜 對應(yīng)的頻率作為估計的接收信號的瞬時頻率/；,和最大功率吟; 步驟5、劃分目標(biāo)位置網(wǎng)格捜索區(qū)間（Xgrid，ygrid)，信號發(fā)射時間的網(wǎng)格捜索區(qū)間tgrid， W及發(fā)射的信號長度的網(wǎng)格捜索區(qū)間Tgrid; 步驟6、對每個網(wǎng)格點rUgrid= (Xgrid，ygrid, tgrid，Tgrid)計算相應(yīng)的信號起始時間tl,s (nigrid )和終止時間11, e ( rilgrid )，其中，Xgrid，ygrid表示網(wǎng)格點的發(fā)射機的位置坐標(biāo)，tgrid表示 信號發(fā)射的時間的網(wǎng)格捜索區(qū)間坐標(biāo)，Tgrid表示發(fā)射的信號長度的網(wǎng)格捜索區(qū)間坐標(biāo)； 步驟7、根據(jù)起始時間和終止時間對應(yīng)的時間區(qū)間，從瞬時頻率Λ中取出對應(yīng)區(qū)間的頻 率值，作為線性調(diào)頻信號的瞬時頻率的估計i：置,?巧，從最大功率，中取出對應(yīng)區(qū)間的 最大功率值; 步驟8、對從L個量測中估計出的L個線性調(diào)頻信號的瞬時頻率的估計i：志,?巧1,.,,,；)取平均， 得到線性調(diào)頻信號的頻率估計值撫對L個量測中取出的最大功率％,,"" (恥y)取平 均再求和，得到區(qū)間功率和ω (rugrid); 步驟9、通過最小二乘參數(shù)擬合從線性調(diào)頻信號的頻率估計值如中估計線性調(diào) 頻信號的初始頻率Λ(恥府)和斜率《:(巧如,）作為估計的線性調(diào)頻信號參數(shù)步驟10、計算最小二乘擬合誤差的區(qū)間平均值Chlgrid)，將最小二乘擬合誤差的區(qū)間平 均值ξ仙grid)與區(qū)間功率和ω (riigrid)相除得到質(zhì)量因子e(nigrid); 步驟11、根據(jù)網(wǎng)格點nigrid二（Xgrid，ygrid，tgrid，Tgrid)和估計的線性調(diào)頻信號參數(shù) 衣?巧喻)，計算網(wǎng)格點nigrid的加權(quán)似然函數(shù)值；運里分號后的rUgrid和 衣(?…)表示該似然函數(shù)的取值郵Igrid和4 (巧?？鞐l件。）； 步驟12、重復(fù)步驟6至11，直到遍歷數(shù)據(jù)平面上所有網(wǎng)格點，得到全數(shù)據(jù)平面的加權(quán)似 然函數(shù)值，并找出取得最大似然函數(shù)值對應(yīng)的網(wǎng)格點的位置作為估計的發(fā)射機所在位置。2. 如權(quán)利要求1所述一種聯(lián)合估計信號參數(shù)和位置的高精度被動定位方法，其特征在 于，網(wǎng)格點化grid的加權(quán)似然函數(shù)值的計算方法是：其中，Η表示對矩陣求共輛轉(zhuǎn)置，e(rugrid)為質(zhì)量因子，r功量測向量，Ri代表第1個接收 機通道上的噪聲協(xié)方差矩陣，馬為根據(jù)信號的起始時間tl,shlgrid)、發(fā)射的信號的長度Tgrid、 和估計的線性調(diào)頻信號參數(shù)衣(巧。W)構(gòu)建的匹配函數(shù)3?[巧iwi,，東(巧在時間上離散化 后，得到的離散化后的匹配函數(shù)；其中，t表示當(dāng)前時間變量，1 = 1，2, . . .，L，T表示接收機的觀測時間長度(運里觀測時 間長度要求盡量長，W保證接收到完整的發(fā)射信號衣(％4]為根據(jù)信號參數(shù) 走?巧g,.w>陜復(fù)出的發(fā)射信號，exp為指數(shù) 函數(shù)。3.如權(quán)利要求1所述一種聯(lián)合估計信號參數(shù)和位置的高精度被動定位方法，其特征在 于，每個網(wǎng)格點nigrid = ( Xgrid，y肝id，tgrid，Tgrid )的信號起始時間11, S (m肝id )和終止時間tl, e (化grid)的計算方法是：tl, e (nigrid) _ tl, s (Hlgrid) +Tgrid ； 其中，C為光速，XI，y功雷達接收站坐標(biāo)。
【文檔編號】G01S5/04GK105974362SQ201610281883
【公開日】2016年9月28日
【申請日】2016年4月29日
【發(fā)明人】易偉, 諶振華, 李帥, 方梓成, 孔令講, 楊曉波
【申請人】電子科技大學(xué)