CN104678386B - 一种利用gnss海面反射信号相关功率探测目标的方法 - Google Patents

Abstract

一种利用GNSS海面反射信号相关功率探测目标的方法，它有两大步骤：步骤一：计算海面反射信号的时延‑多普勒二维相关功率DDM；包括GNSS直射信号的捕获和跟踪；反射信号时延和多普勒频移估计；探测区域的划分和计算反射信号二维相关功率；步骤二：利用压缩感知理论对退化后的DDM进行恢复运算；包括建立基于压缩感知理论的退化后DDM的复原模型；DDM的稀疏性分解；选取高斯型随机矩阵作为测量矩阵和DDM的复原运算。本发明在遥感技术应用领域里有广阔地应用前景。

Description

一种利用GNSS海面反射信号相关功率探测目标的方法

技术领域

本发明涉及一种基于GNSS海面反射信号的目标探测方法，更特别地说，是一种利用GNSS海面反射信号相关功率的时延-多普勒映射(Delay Doppler Mapping，DDM)探测目标的方法，属于遥感技术应用领域。

背景技术

目前的目标探测主要通过雷达来实现，利用定向天线将调制信号的电磁波能量发射到特定区域进行目标搜索，通过雷达接收机处理目标反射的回波，得到目标的特征信息。随着电子技术的不断发展，传统雷达受到电子干扰、反辐射导弹、各种隐身武器的威胁。现实环境中，存在着大量FM广播信号、通信基站信号、电台信号和GNSS卫星信号等，可为建立非合作式双(多)基地无源雷达提供丰富的外辐射源，以应对传统雷达所面临的威胁。

GNSS为用户提供精确导航、定位和授时服务的同时，全天时、全天候地提供高度稳定、长期可用、安全性较高且覆盖全球的微波信号源。随着GNSS的不断发展，利用海面反射后的GNSS信号实现海面特征参数的提取是近年来GNSS-R(GNSS Reflectometer)技术的研究热点。GNSS信号经海面反射后，信号的波形、极化、功率及频率等参量都受反射面的特征所影响，与反射面直接相关；因此可以利用GNSS卫星作为辐射源，建立非合作式双(多)基地无源雷达，通过对海面反射信号的处理来实现目标的探测和定位。

在利用GNSS海面反射信号探测目标的双(多)基无源雷达系统中，当发射机和接收机速度、高度、天线增益等参数固定时，探测区域的反射信号相关功率主要由海面的雷达散射截面(Radar Cross Section，RCS)决定，而海面有无目标时的RCS差异较大，因此可以通过海面反射信号相关功率的时延-多普勒映射(DDM)来实现目标的探测与定位。

发明内容

本发明的目的是提供一种利用GNSS海面反射信号相关功率探测目标的方法，通过对海面反射区域的GNSS反射信号进行相关处理，得到反射信号的时延-多普勒二维相关功率分布，再利用图像处理的思想，将得到的DDM看作是原图像受到某种干扰后的退化图像，采用压缩感知理论对DDM进行恢复运算，得到退化之前的原图像，进而从原图像中提取出目标是否存在和目标位置的信息。

一种利用GNSS海面反射信号相关功率探测目标的方法，包括以下几个步骤：

步骤一：计算海面反射信号的时延-多普勒二维相关功率(DDM)；

具体内容包括：

(1)GNSS直射信号的捕获和跟踪；

通过对直射信号的捕获跟踪，得到信号的C/A码码相位和载波多普勒频移，并完成导航定位解算，获得卫星和接收机的位置信息。

(2)反射信号时延和多普勒频移估计；

根据步骤(1)中卫星和接收机的位置信息，计算镜面反射点的位置，并估计出反射信号相对于直射信号的时延和多普勒频移。

(3)探测区域的划分；

以步骤(2)中镜面反射点的位置为中心，设定探测区域的范围，并对探测区域进行单元格划分。

(4)计算反射信号二维相关功率；

将步骤(3)中划分后的探测区域从空间域转换到时延-多普勒域，根据步骤(2)中反射信号相对于直射信号的时延和多普勒频移信息，利用步骤(1)中直射信号的码相位和载波多普勒辅助产生本地反射信号，并与接收到的反射信号进行相关运算，得到反射信号的DDM。

步骤二：利用压缩感知理论对退化后的DDM进行恢复运算；

具体内容包括：

(1)建立基于压缩感知理论的退化后DDM的复原模型；

在压缩感知理论框架下建立DDM退化模型，将模糊函数作为退化系统的冲击响应函数，把DDM转换为一个列向量；再建立DDM的复原模型，对原DDM进行估计。

(2)DDM的稀疏性分解；

对退化后的DDM进行小波变换，设定小波变换的层数，选择二维小波基矩阵作为稀疏基，对DDM进行稀疏表示。

(3)选取高斯型随机矩阵作为测量矩阵；

(4)DDM的复原运算；

利用贪婪迭代算法对原始DDM稀疏系数进行估计，将其与二维小波基矩阵相乘，并对所得结果进行小波逆变换，得到原始DDM，再从时延-多普勒域转换到空间域，从反射信号相关功率分布的变化中可以探测到该区域是否有无目标，如果有目标，则可以得到目标所在单元格的位置信息。

本发明的优点为：

1、作为收发分置的双(多)基雷达，具有抗反辐射导弹攻击的能力，利用前向散射信号，可对隐形目标实现有效探测；

2、将图像处理思想运用到时延-多普勒映射中，用于目标探测，较好的契合了反射信号DDM的特点；

3、利用压缩感知理论对DDM进行复原运算，与传统的直接反卷积运算相比，减少了运算量和复杂度。

附图说明

图1是本发明的方法流程图；

图2是本发明中探测区域划分实例图；

图3是本发明中压缩感知理论框架下的DDM复原模型；

图中符号说明如下：

图1中的τ_E和f_E分别表示反射信号相对于直射信号的时延和多普勒频移估计值；DDM_de表示反射信号的二维相关功率分布；Φ为压缩感知理论框架下的测量矩阵；θ为对原始DDM稀疏分解后的稀疏系数矩阵，是一个一维的列向量矩阵；W表示对原始DDM稀疏分解的稀疏基矩阵，是一个二维小波基矩阵。图3中m,n分别表示时延和多普勒划分个数，为退化后的DDM，即DDM_de；G(m,n)为通过测量矩阵Φ对DDM_de在加性噪声的干扰下的观测结果；n表示加性噪声；H^-1[·]表示估计过程；为对原始DDM的估计结果。

具体实施方式

下面结合附图对本发明的具体实施方式进行详细说明。

本发明的一种利用GNSS海面反射信号相关功率探测目标的方法，流程如图1所示，包括以下几个步骤：

步骤一：计算海面反射信号的时延-多普勒二维相关功率(DDM)；

由于海面反射信号比较微弱，要直接捕获跟踪反射信号比较困难，本发明采用直射信号辅助反射信号的模式，首先对GNSS直射信号进行捕获跟踪，通过定位解算，得到卫星和接收处理平台的位置信息，计算出镜面反射点的位置，估计出反射信号相对于直射信号的时延和多普勒，结合直射信号的C/A码和载波信息，产生本地复制的反射信号；然后以镜面反射点为中心，对探测区域进行空间域的划分，并转换到时延和多普勒域，再与接收的反射信号进行相关处理，得到海面反射信号的时延-多普勒二维相关功率。

本发明的反射信号相关功率计算方法具体包括：

(1)GNSS直射信号的捕获跟踪；

根据GNSS空间星座图和目标探测区域的几何关系选择要捕获跟踪的卫星，利用传统的方法，通过在时延和多普勒域并行搜索，完成GNSS直射信号的捕获，然后对其进行精确跟踪。接收到的直射信号S(t)表示为：

S(t)＝A(t)D[t-τ(t)]C[t-τ(t)]cos[φ(t)] (1)

其中A(t)表示接收信号幅度电平，D[·]为导航电文数据，C[·]表示卫星的PRN码，φ(t)表示接收信号的载波相位，τ(t)表示码相位的空间传播延迟。在对GNSS直射信号精确跟踪后，进行导航定位解算，获得卫星和接收平台的位置和姿态信息。

(2)反射信号时延和多普勒频移估计；

根据步骤(1)中GNSS卫星和接收平台的位置和速度信息，可以得到接收平台相对于基准水平面的高度h_R和卫星仰角利用反射事件的几何关系计算出镜面反射点的位置，并且估计出反射信号相对于直射信号的时延τ_E为

其中c为光速；反射信号相对于直射信号的载波多普勒频移f_E为

f_E＝[v_t·m-v_r·n-(v_t-v_r)u_tr]/λ (3)

式中，v_t和v_r分别表示GNSS卫星和接收平台的运行速度；m和n分别表示入射和反射方向矢量；λ为载波波长；u_tr表示卫星和接收平台之间的单位方向矢量。

(3)探测区域的划分；

将接收天线的椭圆照射区域近似为一个以镜面反射点为中心的矩形区域，如图2所示，设定探测区域的范围(例如，120km×120km的区域范围)，然后按照单元格对探测区域进行划分(例如，将探测区域划分为1km×1km大小的单元格)。

(4)计算反射信号二维相关功率

利用步骤(1)中直射信号的C/A码和载波信息和步骤(2)中对反射信号相对于直射信号的时延和多普勒估计，产生本地反射信号，与接收到的反射信号进行相关处理，得到反射信号相关功率分布为

$\⟨{\left|Y(\mathrm{\τ},f_d)\right|}^2\⟩=\frac{P_T{\mathrm{\λ}}^2{T_i}^2}{{\left(4\mathrm{\π}\right)}^3}{\∫}_{\stackrel{\→}{\mathrm{\ρ}}}\frac{D^2\left(\stackrel{\→}{\mathrm{\ρ}}\right){\mathrm{\Λ}}^2\left(\mathrm{\τ}\right){\left|S\left(f_d\right)\right|}^2}{{R_T}^2\left(\stackrel{\→}{\mathrm{\ρ}}\right){R_R}^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)}{\mathrm{\σ}}_0\left(\stackrel{\→}{\mathrm{\ρ}}\right)d^2\stackrel{\→}{\mathrm{\ρ}}---\left(4\right)$

其中，

$S\left(f_d\right)=\frac{1}{T_i}{\∫}_0^{T_i}\exp (-2\mathrm{\πi}{f}_d{T}_i{t}^{\′})d{t}^{\′}=\frac{\sin \left(\mathrm{\π}{f}_d{T}_i\right)}{\mathrm{\π}{f}_d{T}_i}\exp (-\mathrm{\πi}{f}_d{T}_i)---\left(5\right)$

式中P_T为卫星信号发射功率，λ为载波波长，T_i为相干积分时间，R_T和R_R分别表示GNSS卫星和接收平台到海面散射点的距离，为天线增益函数，Λ(τ)表示C/A码自相关函数，为海水表面处的双基地散射界面RCS。

作为一种无源雷达，其模糊函数(Woodward Ambiguity Function，WAF)为

WAF＝χ²(τ,f_d)＝Λ²(τ)|S(f_d)|² (6)

可以将反射信号相关功率分布转换成二维卷积形式为

<|Y(τ,f_d)|²>＝χ²(τ,f_d)**∑(τ,f_d) (7)

其中

$\mathrm{\&Sigma;}(\mathrm{\&tau;},f_d)=\frac{P_T{\mathrm{\&lambda;}}^2{T_i}^2}{{\left(4\mathrm{\&pi;}\right)}^3}{\&Integral;}_{\stackrel{\&RightArrow;}{\mathrm{\&rho;}}}\frac{D^2\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}\right){\mathrm{\&sigma;}}_0\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}\right)}{{R_T}^2\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}\right){R_R}^2\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}\right)}\mathrm{\&delta;}(\mathrm{\&tau;}-\mathrm{\&tau;}\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}\right))\mathrm{\&delta;}(f_d-f_d\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}\right))d^2\stackrel{\&RightArrow;}{\mathrm{\&rho;}}---\left(8\right)$

由于海面探测区域的每个点(x，y)都有与之对应的时延多普勒(τ，f_d)，利用雅克比行列式和δ函数的特性，将∑(τ,f_d)从空间域转换到时延多普勒域有

$\mathrm{\&Sigma;}(\mathrm{\&tau;},f_d)=\frac{P_T{\mathrm{\&lambda;}}^2{T_i}^2}{{\left(4\mathrm{\&pi;}\right)}^3}\&CenterDot;\frac{D^2\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}({\mathrm{\&tau;}}_{\mathrm{xy}},f_{d,\mathrm{xy}})\right){\mathrm{\&sigma;}}_0\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}({\mathrm{\&tau;}}_{\mathrm{xy}},f_{d,\mathrm{xy}})\right)}{{R_T}^2\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}({\mathrm{\&tau;}}_{\mathrm{xy}},f_{d,\mathrm{xy}})\right){R_R}^2\left(\stackrel{\&RightArrow;}{\mathrm{\&rho;}}({\mathrm{\&tau;}}_{\mathrm{xy}},f_{d,\mathrm{xy}})\right)}.\left|J(\mathrm{\&tau;},f_d)\right|---\left(9\right)$

其中

$\left\{\begin{array}{c}{\mathrm{\&tau;}}_{\mathrm{xy}}=\mathrm{\&tau;}(x,y)\\ f_{d,\mathrm{xy}}=f_d(x,y)\end{array}\right.---\left(10\right)$

|J(τ,f_d)|为从空间域转换到时延多普勒域的雅克比行列式。可以看出，模糊函数χ²(τ,f_d)和∑(τ,f_d)均为关于时延-多普勒频移的二维矩阵，反射信号相关功率(DDM)分布<|Y(τ,f_d)|²>也是关于时延-多普勒频移的二维矩阵。

步骤二：利用压缩感知理论对退化后的DDM进行恢复运算；

由式(6)可以看出，在将探测区域进行划分后，反射信号相关功率(DDM)是一个二维矩阵，设时延上对探测区域进行M个划分，在多普勒域上进行N个划分，则DDM为一个M×N的二维矩阵。从图像处理的角度看，将∑(τ,f_d)看作是探测区域反射信号的一种数字成像，记为F(m,n)，将模糊函数χ²(τ,f_d)理解为图像的退化函数，记为H(m,n)，那么接收处理得到的海面探测区域反射信号二维相关功率(DDM)可以看作是原图像F(m,n)受到某种干扰后的模糊图像<|Y(τ,f_d)|²>，记为则可将式(6)转化为

$\tilde{G}(m,n)=\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}F(i,j)\&CircleTimes;H(m-i,n-j)---\left(11\right)$

其中“”表示循环卷积。如果对DDM进行复原，得到退化之前的F(m,n)，便可以从中清晰的看出反射信号相关功率的变化情况，进一步提取出是否有目标出现以及目标出现的位置等信息。

本发明的对退化后DDM进行复原运算具体包括：

(1)建立基于压缩感知理论的退化后DDM的复原模型；

为了便于分析计算，将式(10)中的矩阵按辞书式排列，即从第一行起把二维序列的每个行转置向量一个接一个排成单列向量，使用该方法及矩阵运算对离散卷积进行表述，则式(10)的矩阵表达式为

$\tilde{G}=\mathrm{HF}---\left(12\right)$

其中，和F表示MN×1维的列向量，分别由M×N的矩阵和F(m,n)的各行堆积而成，如下所示：

$\tilde{G}={\left[\begin{array}{c}G\left(\mathrm{0,0}\right)\\ G\left(\mathrm{0,1}\right)\\ .\\ .\\ .\\ G(0,N-1)\\ G\left(\mathrm{1,0}\right)\\ G\left(\mathrm{1,1}\right)\\ .\\ .\\ .\\ G(M-1,N-1)\end{array}\right]}_{\mathrm{MN}\&times;1},F={\left[\begin{array}{c}F\left(\mathrm{0,0}\right)\\ F\left(\mathrm{0,1}\right)\\ .\\ .\\ .\\ F(0,N-1)\\ F\left(\mathrm{1,0}\right)\\ F\left(\mathrm{1,1}\right)\\ .\\ .\\ .\\ F(M-1,N-1)\end{array}\right]}_{\mathrm{MN}\&times;1}---\left(13\right)$

而H为MN×MN维矩阵，具有以下特殊结构：

其中，1≤i≤M。可以看出，如果直接使用H来进行F的求解，运算量将会非常大，因此，本发明利用压缩感知理论来对原图像F(m,n)进行估计。

如图3所示，经过测量矩阵Φ，在加性噪声的干扰下得到G(m,n)，根据式(11)，转换成矩阵表示为

$G=\mathrm{\&Phi;}\tilde{G}+n=\mathrm{\&Phi;HF}+n---\left(15\right)$

其中，Φ∈R^K×MN，K＜＜MN，为测量矩阵；n为K×1维加性噪声；G表示通过测量矩阵Φ对退化图像在加性噪声n的干扰下的观测结果，为K×1维列向量。

根据压缩感知理论，对信号F在变换域(基为W)进行稀疏分解所得系数为θ是稀疏的，如果测量矩阵Φ和稀疏基W不相关，那么可以从观测结果中精确重构原始信号。可以得到

G＝ΦHWθ+n,θ＝W^-1F (16)

其中θ为对信号F稀疏分解的系数矩阵，只要估计出θ，就可以重构出F，进而可以得到原图像F(m,n)的估计

(2)DDM的稀疏性分解；

本发明中将小波基支撑的小波域作为图像重构时的处理域，对信号F稀疏表示的基矩阵W选为二维正交小波逆变换矩阵，满足W^TW＝I，则有F＝Wθ。在对信号F进行重构时，先对进行小波变换，本发明中选用得到广泛承认的性能最优的小波滤波器之一的9/7小波滤波器对进行离散小波变换。设在理想情况下，利用基矩阵W对信号F进行稀疏表示时，稀疏度为L，即稀疏系数θ中包含L个非零元素。

(3)选取高斯型随机矩阵作为测量矩阵；

压缩感知理论要求测量矩阵和稀疏基矩阵不相关，由式(15)可知，实际测量矩阵相当于ΦH或者实际的稀疏基矩阵为HW，因此，本发明中利用高斯型随机矩阵Φ对信号进行观测。

(4)DDM的复原运算

本发明中利用正交匹配追踪算法(Orthogonal Matching Pursuit，OMP)对原图像F(m,n)进行重构，OMP算法是一种贪婪迭代算法，能以较高概率重构信号，计算复杂度为O(MNL²)，相对较低。通过贪婪迭代的方法选择感知算子Θ＝ΦHW中与当前误差向量r相关性最大的列，根据选择的列向量求得当前的近似解及新一次的迭代误差，反复迭代，直到迭代误差满足预设要求。感知算子Θ中各列表示为其中j＝1,2,…,MN；迭代计数为t；索引为λ_t；存储计算出的稀疏系数θ向量中的非零元素的集合为Λ_t；两次迭代所得观测量之差，即余量为r^(t)；第t次选出的列向量与t-1次选出的列向量进行正交化后的列向量为η_t；由前t次迭代得到的η_t组成的矩阵为Ω_t。

OMP算法的迭代步骤如下：

输入：感知算子Θ，观测向量G，稀疏度L；

输出：稀疏系数θ的估计重构误差向量r；

step1初始化余量r⁽⁰⁾＝G，Ω₀为空矩阵，集合迭代计数t＝1；

step2找到索引λ_t，使其满足

step3将与正交化得η_t，且Ω_t＝[Ω_t-1,η_t]；

step4用最小二乘法获得新的观测量估计值为G^(t)＝Ω_t(Ω_t ^TΩ_t)^-1Ω_t ^TG，求解θ中的对应

的非零元素集合 ${\mathrm{\&Lambda;}}_t={\mathrm{\&Lambda;}}_{t-1}\&cup;\left\{{\mathrm{\&theta;}}_{{\mathrm{\&lambda;}}_t}\right\};$

step5计算新的余量r^(t)＝G-G^(t)；

step6判断t＜L是否成立，如果成立，则t＝t+1，返回step2；

step7最后由所有根据索引λ_t排列成列向量，即为对稀疏系数向量θ的的估计

得到稀疏系数向量的估计值之后，可以得到对信号F的估计转换成二维矩阵之后，再经过小波逆变换，便可以得到对原图像F(m,n)的精确估计从而完成对DDM的复原运算。

将恢复的DDM从时延-多普勒域转换到空间域，可以通过反射信号相关功率的分布变化情况判断出探测区域是否有目标，如果有目标，则可以得到目标出现的单元格的位置，再通过镜面反射点的位置便可得到出目标在海面出现的具体位置信息。

本发明中，针对有无目标时海面反射信号相关功率分布的特点，利用图像处理的思想来对反射信号相关功率的分布进行处理；根据DDM的稀疏性，采用压缩感知理论来对DDM进行恢复运算，相比于直接进行反卷积运算，本发明方法有效的减少了运算量。

Claims (1)

1.一种利用GNSS海面反射信号相关功率探测目标的方法，其特征在于：包括以下几个步骤：

步骤一：计算海面反射信号的时延-多普勒二维相关功率；

(1)GNSS直射信号的捕获和跟踪；

通过对直射信号的捕获跟踪，得到信号的C/A码码相位和载波多普勒频移，并完成导航定位解算，获得卫星和接收机的位置信息；

(2)反射信号时延和多普勒频移估计；

根据步骤(1)中卫星和接收机的位置信息，计算镜面反射点的位置，并估计出反射信号相对于直射信号的时延和多普勒频移；

(3)探测区域的划分；

以步骤(2)中镜面反射点的位置为中心，设定探测区域的范围，并对探测区域进行单元格划分；

(4)计算反射信号二维相关功率；

将步骤(3)中划分后的探测区域从空间域转换到时延-多普勒域，根据步骤(2)中反射信号相对于直射信号的时延和多普勒频移信息，利用步骤(1)中直射信号的码相位和载波多普勒频移辅助产生本地反射信号，并与接收到的反射信号进行相关运算，得到反射信号的时延-多普勒二维相关功率；

步骤二：利用压缩感知理论对退化后的时延-多普勒二维相关功率进行恢复运算；

(1)建立基于压缩感知理论的退化后时延-多普勒二维相关功率的复原模型；

在压缩感知理论框架下建立时延-多普勒二维相关功率退化模型，将模糊函数作为退化系统的冲击响应函数，把时延-多普勒二维相关功率转换为一个列向量；再建立时延-多普勒二维相关功率的复原模型，对原时延-多普勒二维相关功率进行估计；

(2)时延-多普勒二维相关功率的稀疏性分解；

对退化后的时延-多普勒二维相关功率进行小波变换，设定小波变换的层数，选择二维小波基矩阵作为稀疏基，对时延-多普勒二维相关功率进行稀疏表示；

(3)选取高斯型随机矩阵作为测量矩阵；

(4)时延-多普勒二维相关功率的复原运算；

利用贪婪迭代算法对原始时延-多普勒二维相关功率稀疏系数进行估计，将其与二维小波基矩阵相乘，并对所得结果进行小波逆变换，得到原始时延-多普勒二维相关功率，再从时延-多普勒域转换到空间域，从反射信号相关功率分布的变化中探测到该区域是否有无目标，如果有目标，则得到目标所在单元格的位置信息。