A FOURTH-ORDER CUMULANT BASED RECURSIVE APPROACH TO DOA ESTIMATION WITH SPARSITY CONSTRAINT

Wan Qun; Yang Wanlin

Volume 23 Issue 11

Nov. 2001

Turn off MathJax

Article Contents

Article Navigation > Journal of Electronics & Information Technology > 2001 > 23(11): 1123-1128

WANG Feng, LI Jianqiang, ZHANG Guodong, ZHANG Qi, YANG Dongkai. Maximum Likelihood Estimation of Ocean Wind Vector Using Subsatellite-Observation Spaceborne Global Navigation Satellite System-Reflectometry[J]. Journal of Electronics & Information Technology, 2024, 46(4): 1418-1427. doi: 10.11999/JEIT230464

Citation:

Wan Qun, Yang Wanlin . A FOURTH-ORDER CUMULANT BASED RECURSIVE APPROACH TO DOA ESTIMATION WITH SPARSITY CONSTRAINT[J]. Journal of Electronics & Information Technology, 2001, 23(11): 1123-1128.

Citation:

Wan Qun, Yang Wanlin . A FOURTH-ORDER CUMULANT BASED RECURSIVE APPROACH TO DOA ESTIMATION WITH SPARSITY CONSTRAINT[J]. Journal of Electronics & Information Technology, 2001, 23(11): 1123-1128.

PDF( 1048 KB)

A FOURTH-ORDER CUMULANT BASED RECURSIVE APPROACH TO DOA ESTIMATION WITH SPARSITY CONSTRAINT

Received Date: 1999-09-29
Rev Recd Date: 2000-02-11
Publish Date: 2001-11-19

Abstract

Abstract

One disadvantage with subspace-based DOA estimate method is that it is difficult to incorporate prior knowledge of the signal correlation into the eigen-decomposition. Herein, an estimator that combines ideas from subspace and four-order cumulant matching method is proposed, which makes it possible to incorporate the prior knowledge of signal correlation and sparsity of the DOA estimation into the estimator. Theoretical analysis and numerical results illustrate the robust and high resolution character of the estimator.

FullText(HTML)

1. 引言

海洋深刻影响着全球气候和生态，对人类社会的生产生活至关重要。此外，海洋动力灾害是对全球沿海各国危害最大的自然灾害^[1]。据统计，21世纪以来海洋动力灾害造成我国人员死亡4134人，直接经济损失约2.56×10¹¹元^[2,3]。高密集测量海洋动力参数，提高对海洋动力灾害的预报预警能力对海洋渔业、国际航运、海上交通、海上能源开发和近海居民生产安全具有重要战略意义和迫切现实需求。海面风矢量表征了大气与海洋间的交互作用，是海洋动力参数的重要组成，也是海洋气象预报的基本观测要素^[4,5]。传统的海面风矢量探测手段包括星载辐射计、高度计和散射计等。由于设备较复杂、成本高而不适合大规模卫星组网，使其空间覆盖性和时间分辨率较低^[6,7]。

全球导航卫星系统反射计(Global Navigation Satellite System-Reflectometry, GNSS-R)技术是一种以导航卫星信号作为机会信号源的被动双/多基遥感探测手段 ^[8]。由于其仅需接收机，因此其设备简单、重量轻、功耗低、成本低，较高度计和散射计更便于微卫星搭载实现星座组网。目前该技术已在海面风场^[9,10]、海面高度^[11]、海冰探测^[12,13]和土壤湿度等参数的探测^[14]上得到应用。海面风矢量包含海面风速和风向。随着美国全球导航卫星飓风计划(CYclone GNSS, CYGNSS) ^[15]、捕风A/B^[16]、风云3-E^[17]等卫星的在轨运行，海面风速逐步由数据产品开始向业务化运营转变。但目前的星载GNSS-R主要采用镜像反射信号进行观测和反演，由于其对风向敏感性弱而难以直接用于风向反演。大多工作均集中在通过海浪谱模型建立海面散射系数与海面风向的关系，研究海面风向对反射信号的影响，进而探索GNSS-R反演海面风向的可行性。机载场景的GNSS-R时延波形后沿对海面风向敏感^[18]，可通过匹配理论模型和实测时延波形反演海面风向^[19]。文献[20]提出了机载GNSS-R延迟多普勒图(Delay-Doppler Maps, DDM)对称性偏斜角的概念，并利用仿真和实测数据验证了其对风向的敏感性。文献[21]将DDM非对称性用于星载场景，初步验证了利用DDM反演海面风向的可行性。在镜向观测模式时，由于需较大区域的非镜向散射信号使DDM对风向敏感，因此该方法的空间分辨率较差。上述研究均表明远离镜面反射点的散射信号对海面风向敏感。基于此，非镜像观测模式被提出，利用非镜像散射信号对风向的敏感性反演海面风向^[22,23]。文献[24]利用英国技术验证1号卫星(United Kingdom-TechDemoSat-1, UK TDS-1)采集的全球定位系统(Global Positioning System, GPS)海面后向散射信号，首次论证了星载场景接收后向散射的卫星导航信号的可能性，为非镜像观测模式提供了试验基础。目前就作者所知尚未见星载GNSS-R海面风矢量直接反演的相关研究论述。

本文针对星载GNSS-R海面风矢量反演难点，首先分析非镜向海面散射信号对海面风矢量的敏感性，构建星下点非镜向观测模式；然后定义该模式下对海面风矢量敏感的特征观测量，研究其与海面风矢量的关系，提出星下点观测模式的海面风矢量极大似然估计(Maximum Likelihood Estimation, MLE)反演算法，并提出数值搜索算法提高反演效率；最后搭建非镜像观测模式的星载GNSS-R仿真平台，验证所提算法的可行性，并评估算法性能。

2. GNSS-R非镜像观测模式

GNSS信号在海面发生漫散射，其散射信号中既包含镜向散射成分，也包含非镜向散射成分。如图1所示，镜向观测模式指镜面反射点位于反射天线的有效覆盖区内，且星载接收设备的时延和多普勒窗以镜面反射信号的到达时间和多普勒频率为参考点。接收设备接收到的信号主要为镜面反射点附近的镜向散射成分。非镜向观测模式指天线指向镜面反射点以外区域，接收设备的时延和多普勒窗以天线相位中心指向点的散射信号的到达时间和多普勒频率为参考点。接收的信号为远离镜面反射点的信号。

图 1 镜向模式与非镜向观测模式示意图

下载: 全尺寸图片幻灯片

由于镜向散射信号对海面风向不敏感，为寻找可行的海面风矢量反演模型，建立如所示的坐标系。坐标系的原点位于镜面反射点 ${{\boldsymbol{S}}_{O}}$ ；坐标系的Z轴指向镜面反射点切平面法向；导航卫星与GNSS-R接收卫星位于YOZ平面，且均位于Z轴正半轴侧，但分别位于Y轴正半轴和负半轴侧。值得注意的是本文中的海面风向定义并非为海面风矢量与北向夹角，而是风矢量与X轴正向的夹角 ${\varphi _{\text{w}}}$ ，范围为[0°, 360°]。

图 2 星载GNSS-R本地坐标系示意图

下载: 全尺寸图片幻灯片

在建立的坐标系内，导航卫星 ${\boldsymbol{{\rm T}}}$ 、镜面反射点 ${{\boldsymbol{S}}_{O}}$ 和接收卫星 ${\rm{R}}$ 的位置坐标分别为

$\qquad\quad {\boldsymbol{T}} = \left( {0,{R_{\text{t}}}\sin {\theta _{{\text{sp}}}},{R_{\text{t}}}\cos {\theta _{{\text{sp}}}}} \right)$

(1)

$\qquad\quad {{\boldsymbol{S}}_O} = \left( {0,0,0} \right)$

(2)

${\boldsymbol{R}} = \left( {0, - {R_{\text{r}}}\sin {\theta _{{\text{sp}}}},{R_{\text{r}}}\cos {\theta _{{\text{sp}}}}} \right)$

(3)

其中， ${\theta _{{\text{sp}}}}$ 为镜面反射点处导航信号的入射角或散射角。散射点 ${\boldsymbol{S}} = ({S_x},{S_y},{S_z})$ 处散射向量 ${\boldsymbol{q}}$ 为

${\boldsymbol{q}} = {\boldsymbol{n}} - {\boldsymbol{m}} = {\boldsymbol{\hat x}}{q_x} + {\boldsymbol{\hat y}}{q_y} + {\boldsymbol{\hat z}}{q_z}$

(4)

其中， ${\boldsymbol{\hat x}}$ , ${\boldsymbol{\hat y}}$ 和 ${\boldsymbol{\hat z}}$ 分别为X轴、Y轴和Z轴的单位向量；入射信号和散射信号的单位向量 ${\boldsymbol{m}}$ 和 ${\boldsymbol{n}}$ 分别为

${\boldsymbol{m}} = {\boldsymbol{\hat x}}\sin {\theta _{\text{i}}}\sin {\varphi _{\text{i}}} + {\boldsymbol{\hat y}}\sin {\theta _{\text{i}}}\cos {\varphi _{\text{i}}} - {\boldsymbol{\hat z}}\cos {\theta _{\text{i}}}$

(5)

${\boldsymbol{n}} = {\boldsymbol{\hat x}}\sin {\theta _{\text{s}}}\sin {\varphi _{\text{s}}} + {\boldsymbol{\hat y}}\sin {\theta _{\text{s}}}\cos {\varphi _{\text{s}}} + {\boldsymbol{\hat z}}\cos {\theta _{\text{s}}}$

(6)

其中， ${\varphi _{\text{i}}}$ 和 ${\varphi _{\text{s}}}$ 分别为入射信号和散射信号相对于镜面反射点入射面的方位角； ${\theta _{\text{i}}}$ 和 ${\theta _{\text{s}}}$ 分别为导航信号的入射角和散射角，范围为[–90°,90°]，满足

${\theta _{\text{i}}} = \left\{ \begin{aligned} { \ge 0},&\;\;{{S_y} \le {R_{\text{t}}}\sin {\theta _{{\text{sp}}}}} \\ { < 0},&\;\;{{S_y} > {R_{\text{t}}}\sin {\theta _{{\text{sp}}}}} \end{aligned} \right. \quad$

(7)

${\theta _{\text{s}}} = \left\{ \begin{aligned} { \ge 0},& \;\;{{S_y} \ge - {R_{\text{r}}}\sin {\theta _{{\text{sp}}}}} \\ { < 0},& \;\;{{S_y} < - {R_{\text{r}}}\sin {\theta _{{\text{sp}}}}} \end{aligned} \right.$

(8)

星下点观测模式的观测区域位于卫星星下点。由于星下点的几何坐标易于求解，因此接收处理时的时延和多普勒窗的位置易于确定。

3. 双基电磁散射

利用双基散射系数衡量GNSS信号在散射面上产生的回波强度。在镜像观测模式下，由于对GNSS信号产生散射的主要是大尺度粗糙海面，因此通常采用基尔霍夫近似几何光学模型(Kirchoff Approximation-Geometric Optics, KO-GA)进行散射系数建模。由于海面大小尺度粗糙海面均对非镜像GNSS散射信号产生作用，因此需考虑大小尺度粗糙海面的散射强度。双尺度模型(Two-Scale Model, TSM)将散射面的粗糙度视为大尺度和小尺度粗糙度的叠加。本文采用TSM计算GNSS信号在海面的散射系数。TSM模型的散射系数 ${\sigma _{{\text{TSM}}}}$ 表示为

${\sigma _{{\text{TSM}}}}{\text{=}}{\sigma _{{\text{KA-GO}}}}{\text{ + }}{\sigma _{{\text{SPM}}}}$

(9)

其中， ${\sigma _{{\text{KA-GO}}}}$ 和 ${\sigma _{{\text{SPM}}}}$ 分别为KA-GO散射系数和微扰动法(Small Perturbation Method, SPM)散射系数

${\sigma _{{\text{KA-GO}}}}{{ = \pi }}{\left| \Re \right|^2}\frac{{{{\left| {\boldsymbol{q}} \right|}^4}}}{{q_z^4}}P\left( { - \frac{{{{{q}}_\bot }}}{{{q_z}}}} \right)$

(10)

${\sigma _{{\text{SPM}}}}{\text{ = }}8{\left| {{k^2}\cos {\theta _{\text{i}}}\cos {\theta _{\text{s}}}\Re } \right|^2}S\left( {2k\sin {\theta _{\text{i}}},{\varphi _{{\text{i,w}}}},{u_{10}}} \right)$

(11)

其中， $\Re$ 为反射系数； ${\boldsymbol{q}}$ 为散射向量； ${{{q}}_ \bot }$ 和 ${q_{{z}}}$ 分别为 ${\boldsymbol{q}}$ 的水平分量和垂直分量， ${\varphi _{{\text{i,w}}}}$ 为入射信号相对于海面风向的方位角。如所示为基于所示的坐标系，当入射角 ${\theta _{{\text{sp}}}}$ 为20°时的海面散射系数随散射角的变化。由图3可知，镜向方向具有最强散射，随着不断偏离镜向方向，散射强度逐渐降低。

图 3 海面散射系数随散射角的变化

下载: 全尺寸图片幻灯片

图4为当海面风速为10 m/s，入射角为20°时，不同散射角的海面散射系数随海面风向的变化。当散射角为10°和30°时，即靠近镜向观测模式下，海面散射系数随海面风向微弱周期性变化。这说明镜向信号难以反演海面风向。当散射角偏离镜面反射时，即散射角为–30°, –10°和50°时，海面散射系数随海面风向呈显著性的周期性波动，即非镜向散射信号对海面方向敏感，能够反演海面风向。

图 4 不同散射角的海面散射系数随海面风向的变化

下载: 全尺寸图片幻灯片

海面散射系数同时受海面风速和风向影响。图5为散射角分别为20°(镜向观测模式)和0°(星下点观测模式)时不同海面风速下的海面散射系数与海面风向的关系。由图5可知，海面风速对镜向散射系数的影响远大于海面风向，且散射系数随风速的增加而降低；非镜向散射系数同时受海面风速和风向影响，随海面风速增加而增加，随海面风向呈周期性波动，且不同海面风速下随风向的波动幅度不同。

图 5 镜向和非镜向散射系数与海面风速和海面风向的关系

下载: 全尺寸图片幻灯片

4. 星下点观测模式信号链路

星载GNSS-R功率链路通过双基雷达方程进行计算^[25]

${P_{\text{r}}} = \frac{{{P_{\text{t}}}{G_{\text{t}}}{\lambda ^2}{G_{\text{r}}}\sigma A}}{{{{(4{\pi})}^3}R_{\text{t}}^{\text{2}}R_{\text{r}}^{\text{2}}}}$

(12)

其中， ${P_{\text{r}}}$ 和 ${P_{\text{t}}}$ 分别为接收信号功率和发射信号功率； ${G_{\text{r}}}$ 和 ${G_{\text{t}}}$ 分别为接收天线增益和发射天线增益。GNSS散射信号在接收处理时首先通过相干积分，然后对连续复数相干值进行非相干累加。相干积分和非相干累加处理后的信噪比为

${\text{SNR = }}10\lg \left(\frac{{{P_{\text{r}}}}}{{{P_{{\text{noise}}}}}}\right) + {G_{{\text{coh}}}}({T_{{\text{coh}}}}) + {G_{{\text{incoh}}}}({N_{{\text{incoh}}}})$

(13)

其中， ${G_{{\text{coh}}}}$ 和 ${G_{{\text{incoh}}}}$ 分别为相干积分增益和非相干累加增益； ${T_{{\text{coh}}}}$ 为相干积分时间； ${N_{{\mathrm{incoh}}}}$ 为非相干累加次数。相干积分增益可表示为

${G_{{\text{coh}}}}({T_{{\text{coh}}}}) = 10\lg \left(\frac{{{T_{{\text{coh}}}}}}{{{T_{{\text{chip}}}}}}\right)$

(14)

其中， ${T_{{\text{chip}}}}$ 为伪随机码的宽度。非相干累加增益可表示为^[25]

$\begin{split} {G_{{\text{incoh}}}}({N_{{\text{incoh}}}}) =\,& 10\lg ({N_{{\text{incoh}}}}) \\ & - 10\lg \left(\frac{{1 + \sqrt {1 + 9.2{N_{{\text{incoh}}}}/{D_{\text{c}}}} }}{{1 + \sqrt {1 + 9.2/{D_{\text{c}}}} }}\right) \end{split}$

(15)

其中， ${D_{\text{c}}}$ 为检测因子，是虚警概率和检测概率的函数。当虚警概率和检测概率分别为10^–7和98%时， ${D_{\text{c}}}$ 为26.3。假设散射信号相关功率中仅存在热噪声，利用表1中的参数分析星下点非镜向观测模式散射信号信噪比与入射角的关系。如图6所示，由于入射角越大，星下点区域距镜面反射点越远，散射信号越弱，因此信噪比随入射角的增大而降低；当入射角大于18°时信噪比随海面风速增加而增大，变化趋势与镜向观测模式相反。这是由于海面风速越大，非镜向散射分量强度越大。

表 1 星下点非镜向配置参数表

符号	参数	值
${P_{\text{t}}}$	发射信号功率	26.8 W
${G_{\text{t}}}$	发射天线增益	12.1 dB
${h_{\text{t}}}$	发射机高度	20 200 km
${h_{\text{r}}}$	接收机高度	510 km
${G_{\text{r}}}$	接收天线增益	12.1 dB
${T_{{\text{coh}}}}$	相干积分时间	1 ms
${N_{{\text{incoh}}}}$	非相干累加次数	1 000
${D_{\text{c}}}$	检测因子	26.3
${f_{\text{B}}}$	接收机带宽	2.5 MHz
${T_{{\text{eff}}}}$	等效温度	25°C
${\theta _{\text{i}}}$	入射角	[0,90°]
${\varphi _{\text{w}}}$	风向	90°
${u_{10}}$	风速	5～20 m/s

下载: 导出CSV

| 显示表格

图 6 星下点非镜向观测模式散射信号信噪比与入射角的关系

下载: 全尺寸图片幻灯片

通常当信噪比低于0 dB时，GNSS-R接收机接收的散射信号相关功率不可信^[26]。由表1所示参数计算的功率链路可知，入射角低于35°可保证信噪比高于0 dB，星载GNSS-R接收机可成功接收星下点非镜向散射信号。提高反射信号天线增益可进一步扩大入射角范围和接收信号的信噪比。尽管在星下点非镜向观测配置下星下点散射的信号弱，但仍可接收、处理特定入射角范围内的散射信号用于反演风向。

5. 海面风矢量MLE反演模型

星载GNSS-R的基本观测量是DDM，其描述了信号功率在时延-多普勒域内的分布。由于空间域和时延-多普勒域之间存在映射关系，因此DDM也描述了信号功率在空间域的分布。本文以星下点附近特定时延-多普勒窗内的DDM均值(Delay-Doppler Map Average, DDMA)作为对风速和风向敏感的观测量。仿真得到不同风速、风向和入射角的星下点DDM数据集，并计算得到对应的DDMA。图7为星下点DDMA与信号入射角、海面风速和风向的关系。由图7可知，DDMA是海面风向、海面风速和入射角的3元函数，即为反演海面风速和风向，需建立一个3参数的经验地理模式函数(Geophysical Model Function, GMF)。

图 7 星下点非镜向DDMA与海面风向的关系

下载: 全尺寸图片幻灯片

采用如式(16)的三角函数表示DDMA与风向、风速和入射角的经验GMF

${\text{DDMA}} = a({u_{10}},{\theta _{\text{i}}}) + b({u_{10}},{\theta _{\text{i}}}) \cdot \cos (w{\varphi _{\text{w}}}{\text{ + }}\Delta \varphi )$

(16)

其中， $a(\cdot)$ , $b(\cdot)$ , $w$ 和 $\Delta \varphi$ 为拟合参数。由图7可知，上述的三角函数能很好的拟合DDMA与海面风速、风向及入射角的关系。在经验GMF中DDMA是星载GNSS-R的已知特征观测量数据。入射角根据收发卫星位置计算得到，而海面风矢量为待求解未知变量。当每一个风矢量观测区域内有两个独立的DDMA观测值，可根据经验GMF构成两个方程求解得到风矢量。在如图8所示的多星星下点观测配置下，星载GNSS-R卫星配置一个反射信号天线，位于卫星正下方，垂直向星下点照射，在同一观测海域同时接收到2颗及以上导航卫星的星下点散射信号，利用多颗方位角不同的导航卫星散射信号实现海面风矢量反演。

图 8 多星非镜向观测示意图

下载: 全尺寸图片幻灯片

5.1 极大似然估计

经验GMF与海面风矢量之间是复杂的非线性关系，若采用直接求逆的方法求解海面风矢量，难度大、准确度低。为解决这一“模型已定，参数未知”问题，本文采用MLE直接处理多颗导航卫星星下点散射信号实现海面风矢量反演，将先前概率密度函数的风矢量作为变量，寻找使似然函数最大的风矢量。在噪声条件下DDMA可表示为

${\text{DDM}}{{\text{A}}_i} = {\text{DDMA}}_i^0({u_{10}},{\varphi _{\text{w}}},{\theta _{{\text{i}}i}}) + {\varepsilon _i}({u_{10}},{\varphi _{\text{w}}},{\theta _{{\text{i}}i}})$

(17)

其中， ${\text{DDM}}{{\text{A}}_i}$ 为第i个DDMA测量值； ${\text{DDMA}}_i^0$ 为第i个无噪声DDMA，在极大似然估计过程中为模型DDMA； ${u_{10}}$ 和 ${\varphi _{\text{w}}}$ 分别为海面风速和风向； ${\theta _{{\text{i}}i}}$ 为第i个测量值DDMA的入射角； ${\varepsilon _i}$ 为各种随机噪声引起的随机误差，假设满足均值为0，方差为 ${V_{{\varepsilon _i}}}$ 的高斯分布，即 ${\varepsilon _i}{\text{～}}N(0,{V_{{\varepsilon _i}}})$ 。对给定的海面风矢量， ${{\mathrm{DDMA}}_i}$ 和模型预测值 ${\text{DDMA}}_i^0({u_{10}},{\varphi _{\text{w}}},{\theta _{{\text{i}}i}})$ 间的残差 ${R_i}$ 定义为

$\begin{split} {R_i}({u_{10}},{\varphi _{\text{w}}},{\theta _{{\text{i}}i}})& {\text{ = DDM}}{{\text{A}}_i}-{\text{ DDMA}}_i^0({u_{10}},{\varphi _{\text{w}}},{\theta _{{\text{i}}i}})\\ & {\text{ = }}{\varepsilon _i}({u_{10}},{\varphi _{\text{w}}},{\theta _{{\text{i}}i}})\\[-1pt] \end{split}$

(18)

残差 ${R_i}$ 为均值为0，方差为 ${V_{{R_i}}}$ 的高斯分布。假设在同一观测海域内，海面风矢量为 ${\boldsymbol{W}}{\text{ = (}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}$ 。星载GNSS-R同时接收N颗导航卫星的星下点散射信号，并得到相应的DDMA测量值。由于各DDMA测量值相互独立，即残差 ${R_i}$ 相互独立，因此残差的联合条件概率密度函数为

$p({R_1},{R_2},\cdots,{R_N}|{\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}) = \prod\limits_{i = 1}^N {p({R_i}|{\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}})}$

(19)

当测量值 ${\text{DDM}}{{\text{A}}_1},{\text{DDM}}{{\text{A}}_2},\cdots,{\text{DDM}}{{\text{A}}_N}$ 给定时， $p({R_1},{R_2},\cdots,{R_N}|{\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}})$ 是参数 ${\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}$ 的函数，记为

$\begin{split} & L({\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}|{\text{DDM}}{{\text{A}}_1},{\text{DDM}}{{\text{A}}_2},\cdots,{\text{DDM}}{{\text{A}}_N}) \\ & \quad = p({R_1},{R_2},\cdots,{R_N}|{\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}) = \prod\limits_{i = 1}^N {p({R_i}|{\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}})} \end{split}$

(20)

函数 $L({\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}|{\text{DDM}}{{\text{A}}_1},{\text{DDM}}{{\text{A}}_2},\cdots,{\text{DDM}}{{\text{A}}_N})$ 即为似然函数，求 ${\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}$ 使似然函数取最大值，即满足式(21)

$\begin{split} & L({\text{(}}{\hat u_{10}}{\text{,}}{\hat \varphi _{\text{w}}}{\text{)}}|{\text{DDM}}{{\text{A}}_1},{\text{DDM}}{{\text{A}}_2},\cdots,{\text{DDM}}{{\text{A}}_N})\\ & \quad{{ = }}\mathop {\sup }\limits_{{\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}} \in ((0,{u_{10\max }}],[0,360])} \{ L({\text{(}}{u_{10}}{\text{,}}{\varphi _{\text{w}}}{\text{)}}\\ & \qquad\; |{\text{DDM}}{{\text{A}}_1},{\text{DDM}}{{\text{A}}_2},\cdots,{\text{DDM}}{{\text{A}}_N})\} \end{split}$

(21)

对应的 ${\text{(}}{\hat u_{10}}{\text{,}}{\hat \varphi _{\text{w}}}{\text{)}}$ 即为求解的海面风矢量。

5.2 数值搜索算法

由于似然函数形式复杂，难以直接求得海面风矢量解，因此本文采用数值搜索方法得到海面风矢量解。式(16)所示的经验GMF和式(20)所示的似然函数均为非线性的，在海面风矢量数值搜索时难以预测似然函数局部最大值的具体位置和数量。根据DDMA与海面风向呈简谐函数特点，参考文献[27]提出的散射计极大似然估计海面风场反演算法，本文提出如下搜索算法：

(1)取海面风向为0°，以给定的起始风速6 m/s为搜索起点，在风速区间按照设定的风速搜索间隔寻找使似然函数取最大值的风速，并将对应的似然函数值和海面风向、海面风速记录下来。风速搜索的具体步骤为：分别计算风速起始点和右边相邻点的似然函数值，比较两者大小，如果起始点的似然函数值小于相邻点的似然函数值，则继续向右搜索，反之向左搜索，直到找到使似然函数取最大值的点，并记录似然函数值、海面风速和风向。

(2)令海面风向值增加一个搜索间隔，风速搜索起点为上一个风向搜索到的风速，按照与第1步相同的搜索步骤寻找似然函数取最大值的点，并记录似然函数值、海面风速和风向。

(3)重复第2步的操作，将风向区间0°～360°搜索完，共得到361组似然函数值、海面风速和风向。

(4)根据第3步得到的结果在风向区间0°～360°寻找似然函数的局部极大值，并记录似然函数值、海面风速和风向。

(5)将第4步搜索得到的似然函数局部极大值从大到小排序，海面风矢量的模糊解为前四个局部极大值对应的海面风速和风向。其中，第1模糊解为排名第1的海面风矢量模糊解。

上述数值风矢量搜索算法利用了似然函数在风矢量 ${\text{(}}{u_{10}}{\text{,}}{\varphi _{\mathrm{w}}}{\text{)}}$ 2维空间的分布特征，避免了在整个海面风速区间维度逐点搜索似然函数最大值，搜索效率高。需要注意的是，在实际应用中要综合考虑反演精度和搜索效率来设定合适的海面风速、风向搜索间隔和海面风速起始点。

6. 算法验证

利用星下点观测模式的GNSS-R仿真平台得到不同方位角、入射角、海面风速、海面风向下的星下点DDM数据集。数据集的风速范围为限制在[2 m/s, 25 m/s]，风向范围为[0°, 360°]。

6.1 双星观测

双星观测中两颗不同入射角和方位角的导航卫星散射信号被用于反演海面风矢量。利用海面风矢量极大似然估计反演算法处理双星观测验证数据(式(19)中N=2)。为第1模糊解对应的海面风速和风向。由可知，随着风速增加，海面风速反演结果离散度增大，即海面风速越高，风速反演精度越低。海面风向存在4个海面风向模糊解 ${\varphi _{{\text{w}}1}}$ , ${\varphi _{{\text{w}}2}}$ , ${\varphi _{{\text{w}}3}}$ 和 ${\varphi _{{\text{w4}}}}$ ，且4个模糊解(见图10)满足

图 9 双星观测海面风矢量反演结果

下载: 全尺寸图片幻灯片

图 10 海面风向反演多模糊解示意图

下载: 全尺寸图片幻灯片

$\begin{split} {\varphi _{{\text{w}}2}}& {\text{ = }}{360^ \circ }{{ - }}{\varphi _{{\text{w1}}}},{\varphi _{{\text{w3}}}}{\text{ = }}{180^ \circ }{\text{ + }}{\varphi _{{\text{w1}}}},{\varphi _{{\text{w4}}}}\\ & {\text{ = }}{180^ \circ }{{ - }}{\varphi _{{\text{w1}}}} \end{split}$

(22)

6.2 3星观测

图10海面风向反演多模糊解示意图。3星观测中每个时刻接受3颗入射角和方位角不同的导航卫星在同一星下点海域的散射信号。利用海面风矢量MLE反演算法处理3星观测验证数据(式(19)中N=3)。如图11所示，相比双星观测，海面风速的反演精度得到大幅提升，且3星观测消除了由观测几何对称性引起的模糊解，使海面风向对应的模糊解由4个降为2个。这两个模糊度由海浪谱对称性导致，无法通过增加导航卫星数量方式消除。一种有效消除该模糊度的方法是多卫星遥感数据融合，即借助其他气象辅助方法确定正确风向。

图 11 3星观测海面风矢量反演结果

下载: 全尺寸图片幻灯片

6.3 算法性能影响因素分析

图12为双星观测的海面风速和风向的均方根误差随信噪比的变化，其中海面风向取4个模糊解中最接近真实风向的解，海面风速则为该模糊解对应的解。由图12可知，风速和风向的均方根误差均随信噪比增加而减小，即信噪比越高，海面风矢量反演精度越高。

图 12 双星观测海面风矢量反演精度与信噪比的关系

下载: 全尺寸图片幻灯片

图13为3星观测的海面风速和风向的均方根误差随信噪比的变化。由图13可知，和双星观测变化规律一致。当信噪比高于11 dB时，风速的反演精度优于2 m/s，风向反演精度优于15°。由图6可知，当导航卫星入射角小于17°时，信噪比大于11 dB。当限制入射较以满足信噪比要求时，可见卫星数下降导致空间采样率下降。因此增大反射信号接收天线的增益是同时兼顾反演精度和空间采样率的最有效解决方案。值得注意的是当信噪比低于7 dB时，3星观测的海面风矢量反演精度低于双星观测，尤其是海面风向反演精度。这是因为当信噪比较低时，由于信号质量较差，导航卫星信号越多，似然函数的不确定性越大，导致反演精度降低。

图 13 3星观测海面风矢量反演精度与信噪比的关系

下载: 全尺寸图片幻灯片

7. 结论

星载GNSS-R技术已被扩展到诸多应用领域，尤其海面风场反演逐渐趋于业务化。目前针对星载GNSS-R技术的研究主要集中在镜向观测模式。由于海面镜向散射信号对海面风向敏感性较弱，难以反演海面风向。本文提出星下点观测模式，提出同时反演海面风速和风向的MLE反演方法。本文首先研究了非镜向海面散射系数与海面风速、风向及入射角的关系；然后构建了多卫星的星下点观测模式，提出了海面风矢量MLE反演算法，通过利用多颗导航卫星的星下点散射信号反演得到了海面风矢量；最后利用搭建的星载星下点GNSS-R仿真平台验证了所提算法的可行性，并评估了其性能。结果表明：(1)所提算法可成功反演海面风矢量，且信噪比越高反演精度越高；(2)观测几何关系和海浪谱的对称性导致了海面风向有四个模糊解，通过增加导航卫星的数量可消除观测几何对称性导致的风向模糊解，而无法消除由海浪谱对称性引起的180°风向模糊。由于星下点非镜向散射信号功率低于镜向散射信号，使用高增益反射信号天线、研究GNSS-R高灵敏度接收机、完善星载数据定标和校正方法，可在一定程度上提升星下点非镜向海面风向反演算法的精度。

References(1)

References

魏平,肖先赐,李乐民,基于四阶累积量特征分解的空间谱估计测向方法,电子科学学刊,1995,17(3),243-249[2]M. Jansson, B. Goransson, B. Ottersten, A subspace method for direction of arrival estimationof uncorrelated emitter signal, IEEE Trans. on SP, 1999, SP-47(3), 945-956.[3]魏平,肖先赐,基于四阶累积量的阵列扩展,电子科学学刊,1997,19(6),745-750[4]I.F. Gorodnitsky, B. D. Rao, A recursive weighted norm minimization algorithm and its appli-cation, in Proc. Int. Conf., ASSP, Minneapolis, Apr., 1993, 456-459.[5]B.D. Rao, L. K. K. Delgado, An affine scaling methodology for best basis selection, IEEE Trans.on SP, 1999, SP-47(1), 187-200.[6]P.C. Hansen, Analysis of discrete ill-posed problems by means of L-curve, SIAM Rev., 1992, 34(4), 561-580.[7]I.F. Gorodnitsky, B. D. Rao, Sparse signal reconstruction from limited data using FOCUSS(FOCual Underdetermined System Solver), A re-weighted minimum norm algorithm, IEEE Trans.on SP, 1997, SP-45(3), 600-615.[8]D.A. Linebarger, Redundance averaging with large array, IEEE Trans. on SP, 1993, SP-41(4), 1707-1710.

Relative Articles

Supplements(0)

Cited By

Cited by

Periodical cited type(3)

1.	贺靖，谭鸽伟. 基于ALPFT参数估计的运动目标高分辨率成像. 航空学报. 2025(01): 260-274 .
2.	稂时楠，许奔，崔祥斌. 基于快速后向投影的超宽带冰雷达数据成像算法. 电子与信息学报. 2022(04): 1249-1256 . 本站查看
3.	陈溅来，熊毅，徐刚，张俊超，杨德贵，梁步阁. 基于子图像变标的非线性轨迹SAR成像及其自聚焦方法. 雷达学报. 2022(06): 1098-1109 .

Other cited types(1)

Proportional views

Proportional views

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Get Citation

PDF

XML

Article Metrics

Article views (2040) PDF downloads(479)

A FOURTH-ORDER CUMULANT BASED RECURSIVE APPROACH TO DOA ESTIMATION WITH SPARSITY CONSTRAINT