岸基多星多参数北斗GEO 反射信号海面风速反演

王峰，杨东凯，李杰，张波

（北京航空航天大学 电子信息工程学院，北京 100191）

海面风速作为反映海洋表面状态的物理参数之一，是影响海水运动以及气候变化的重要因素. 近岸海域与人类关系密切，存在大量海上交通运输、渔业、养殖等人类活动. 有效监测近岸海域风速对近岸海域的经济发展、海上安全作业具有重要的现实意义[1].常用的海面风速探测手段主要分现场和遥感监测两种类型. 现场监测主要有浮标、海洋站等，但因其建设成本高，不利于在近岸海域大范围部署. 高度计、散射计以及辐射计等卫星遥感技术能够提供全球尺度的风速，但由于时间和空间分辨率的限制，无法有效监测近岸海域的风速. 机载遥感技术灵活度高，但持续性差，无法提供长时间的连续观测. 岸基遥感平台观测区域固定、且可长期稳定观测，可作为星载和机载遥感技术的补充手段，提供近岸海域风速. 欧美国家已在沿海区域部署了各自的岸基海洋雷达监测网络，用于海洋安全、溢油响应、海啸警报等[2].

全球导航定 位系 统(global navigation satellite system, GNSS)不仅可为用户提供实时的位置、速度以及时间信息，而且可作为一种新型的双基地雷达信号源. 该技术被称作GNSS-Reflectometry (GNSS-R),通过接收、处理直达和地球表面反射的GNSS 信号探测地球表面物理参数. 被动无源的特点使其功耗和成本低，利于在近海岸，尤其野外供电难的区域进行大规模部署. 自1993 年Martin-Neria 提出Passive Reflectometry and Interferometry System (PARIS)[3]概念后，大量理论和试验论证该技术已可用于海面风速[4]、有 效 波 高[5-7]、海 面 高 度[8]、海 冰[9]、土 壤 湿度[10-11]等参数的探测. 地球同步轨道（GEO）卫星相对于地球表面固定，在岸基场景可提供稳定的观测几何，能够实现长期稳定的观测. 利用同步卫星发射的电磁信号作为外辐射源的探测技术成为一种研究热点[12-15]. 随着我国北斗系统的建成，将其应用于探测地球表面物理参数成为一种趋势. 同其他同步卫星类型，北斗系统的地球同步轨道（GEO）卫星在岸基场景能够实现长期稳定的观测，为我国近岸海域的 风 速 监 测 提 供 有 效 手 段[16-19]. LI 等[16]首 次 在 阳 江开展了利用北斗B1I 海面反射信号探测风速的验证试验，结果表明时延波形特征参数可反演高达35 m/s的风速，在“尤特”和“飞燕”台风期间的风速反演精度为3.2 m/s. 杨东凯等[19]研制了北斗B1I 反射信号海面风速实时反演系统，外场测试结果表面在低风速场景，通过多星平均和平滑处理风速反演精度可达1.28 m/s. 随着我国北斗-3 的建成，高码率的B3I 信号已在北斗GEO 卫星上成功播发. 目前尚未针对B3I 信号进行风速反演的研究. 在岸基场景，直射对反射信号的干扰是最重要的误差源之一. 提高观测高度使反射信号相对于直射信号的时延增大是降低直射对反射信号干扰的方法之一. 由于B3I 信号相对于B1I 信号有更宽的信号带宽，距离分辨率高，因此B3I 信号有更好的抗直射信号干扰的特性.当高度角为30°时，B1I 信号在观测高度高于146.63 m时，反射信号不受直射信号干扰；而B3I 信号只需观测高度高于29.32 m. 同一时延波形存在多个特征参数对风速敏感，且多颗GEO 卫星信号稳定在同一反射天线视场内，为多星多参数融合的最优反演提供了必要条件. 本文针对北斗GEO 卫星的B3I 反射信号开展岸基海面风速反演论证，重点研究多星多参数融合的海面风速反演算法，提高风速反演精度和稳定性.

1 GNSS 反射信号模型

1.1 几何模型

岸基场景建立如图1 所示的本地坐标系，其中坐标系原点位于镜面反射点处，Z轴指向镜面反射的法向方向，Y轴在反射面内指向接收机方向，X轴和Y、Z轴构成笛卡尔坐标系. 在本地坐标系内，接收机和GNSS 卫星的坐标为

图1 GNSS 反射信号几何关系Fig. 1 Geometry model of GNSS-R

1.2 时延波形模型

岸基GNSS-R 的基本观测量是GNSS 反射信号一维时延波形，表示了反射信号相关功率在时延上的分布情况. 常用的时延波形模型是基于几何光学和基尔霍夫近似推导的Z-V 模型[20]：

式中：Pt和Gt分别为GNSS 卫星的发射功率和发射天线增益； λ为信号波长；Tcoh为相干积分时间；Gr(r)为接收天线在散射单元r处的天线增益；Rt(r)和Rr(r)分别为GNSS 卫星和接收机到散射单元r的距离；f0为信号载波频率；f(r)为 散射单元r处的多普勒频率；c为电磁波传播速度； σpq(r)为 散射单元r处的散射系数. 从上式可知，GNSS 在粗糙海面经历了3 种空域滤波. 第一种是天线滤波，抑制了天线主瓣外的GNSS 反射信号；第二种是时延滤波，滤除了相对于本地信号的时延大于伪随机码码长的散射信号；第三种是多普勒滤波，滤除了相对于本地信号的多普勒绝对值大于 1/Tcoh的散射信号. 在岸基场景下，海面反射的GNSS 信号既包含了相干成分，也包括了非相干成分. 相干和非相干成分的占比与海态、高度角 相 关[21-22]. 为 了 抑 制 或 消 除 时 延 波 形 中 直 射 信 号和相干成分的干扰，相干累计和非相干累加结合的方 法 已 被 提 出[16-23]. 本 文 采 用LI 等[16]提 出 的 方 法 抑制相干成分的干扰，Nincoh次非相干累计的时延波形表示为

2 风速敏感观测量

针对岸基场景，已有多种海面风速反演模型被提出，可大致分为两类. 第一类利用GNSS 反射信号的相干特征，定义反射信号的相关性度量，主要有相关时间[5]、有效非相干次数[21]以及相干非相干比[22].另一类通过定义时延波形对海态敏感的特征参数进行海面风速的反演，主要有波形面积[6]、后延相对幅度[16]、功率分布比[24]本文采用相关时间和波形面积作为基础特征参数进行海面风速的反演.

2.1 相干时间

时间序列的相干时间定义序列自相关函数的积分宽度为[25]

式中：Ir和Qr为GNSS 反射信号的同相和正交支路的相关值；Id和Qd为GNSS 直达信号的同相和正交支路的相关值，用于消除导航电文、信号功率等对反射信号的影响. 当海面均方高度服从分布，GNSS 反射信号复数相关值的自相关函数近似为高斯函数[5]：

1) 利用北斗B3I 反射和直射信号复数相关值得到式(8)所示的时间序列；

2) 对时间序列进行自相关处理得到反射信号自相关函数；

2.2 归一化时延波形面积

归一化时延波形面积定义为归一化时延波形超出阈值的积分面积：

图2 波形面积与风速的关系Fig. 2 Relationship between waveform area and wind speed

3 多星多参数模型

基于相关时间和波形面积均对风速敏感以及多颗GEO 卫星的反射信号在反射天线视场内的事实，本文利用主成分分析（principal component analysis, PCA）和 线 性 无 偏 最 小 方 差（linear unbiased minimum variance, LUMV）估计进行多星多参数的风速反演模型的研究.

3.1 主成分分析

3.2 线性无偏最小方差估计

4 北斗GEO B3I 信号试验

4.1 试验场景

为了验证北斗B3I 反射信号反演海面风速，在山东省东营市青东五验潮站（北纬37°26′51″；东经119°0′36″）开展了外场试验. 试验时间为2021 年6 月16 日16 时至6 月19 日18 时. 试验场景如图3 所示，直射和反射天线距海面约12 m. 直射天线采用通用的右旋圆极化全频段导航天线. 反射天线采用增益为94 dB，波束宽度为±60°的方向性左旋圆极化天线. 反射天线向下45°指向海面，其方位角为250°. 北斗GEO 卫星PRN01，PRN02以及PRN03 位于反射天线视场内被长期稳定接收. 验潮站装配的风杯式风速测量仪提供了试验期间的同比风速. 试验期间采用的4 通道中频数据信号采集器的量化位数为4 bit，采样率为32.738 MHz，中频为0.1271 MHz.为了减小数据存储量，B3I 直反射信号每隔10 min 采集一组，每组采集时长为2 min.

图3 北斗B3I 试验场景Fig. 3 Experiment scenario of BeiDou B3I signals

4.2 数据处理

GNSS 反射信号处理通常采用先跟踪直射跟踪，然后通过跟踪的码和载波相位信息辅助反射信号开环跟踪的方法. 由于GEO 卫星在岸基场景下相对接收平台固定，且信号接收稳定，因此一方面无需动态捕获直射信号确认信号有无，另一方面直反射信号的载波多普勒相对稳定，无需利用闭环跟踪动态估计直射信号多普勒频率. 本文采用如图4 所示的简化处理结构，首先利用采集系统的中频产生数字载波对直反射信号进行载波剥离；然后利用快速傅里叶变换对直反射信号进行频域码并行相关运算[28].相干积分时间设置为北斗B3I 信号伪码周期1 ms.1.2×105个连续的直反射信号的复数相关值得到式(8)所示的时间序列，并对时间序列进行自相关求解反射信号相干时间. 为了得到归一化时延波形，1.2×105个连续的复数相关波形被非相干累加，总处理时长为2 min.

图4 北斗GEO 直反射信号处理结构Fig. 4 Architecture of processing direct and reflected signals from BeiDou GEO satellites

图5(a)为北斗GEO 卫星PRN01 的直射和反射信号时延波形. 反射信号相对于直射信号的时延为15.30 m，对应的接收平台高度为12.05 m，和实际高度相一致. 由于岸基反射信号中包含大量的相干成分，因此利用式(6)进行相干成分抑制后反射信号的时延波形峰值下降. 反射信号相对于直射信号的时延小于一个B3I 码片，直射信号对反射信号形成干扰. 相比于未经抑制相干成分的时延波形，式(6)有效抑制了直射信号干扰，反射信号时延波形的前沿更陡峭. 图5(b)为北斗GEO 卫星PRN01 的反射信号复数相关值的自相关函数. 由图可知，反射信号的自相关函数符合式(9)所示的高斯分布，相应的相关时间为193.9 ms.

图5 时延波形和自相关函数Fig. 5 Delay waveform and Autocorrelation function

4.3 结果分析

图6 给出了由北斗GEO 卫星PRN02 的反射信号估计的相干时间和波形面积与风速的散点关系.由图可知，相关时间和波形面积分别与风速呈正比例和反比例关系，且与式(12)和(14)所示的关系相一致. 相干时间的拟合优度为0.52，未消除直射干扰和消除直射干扰的归一化时延波形面积和风速的拟合优度分别为0.19 和0.35. 相干时间的拟合优度明显高于归一化时延波形面积的，能够提供更好的风速反演. 消除直射信号干扰后，归一化时延波形面积的拟合优度明显得到提高. 这说明直射信号干扰是风速反演的重要误差源之一. 由于直射信号的干扰使得波形面积并未与风速呈现反比例关系，抑制相干成分后，波形面积和风速呈现了较好的反比例关系.

图6 散点图Fig. 6 Scatter diagram

表1 给出了各反演算法反演风速的均方根误差.从结果可知，由于式(6)并不能完全抑制相干成分，因此相比于相干时间，波形面积获得较差反演结果.利用3.1 节所述流程对各特征参数组成的矢量序列进行主成分分析. 图7 给出了第一和二主成分与风速的关系. 从图可知，第一主成分和风速呈现了反比例关系，可表示为指数衰减函数.

图7 第一和二主成分与风速的关系Fig. 7 Relationship of first and second components versus wind speed

其中，aws和bws为拟合参数. 第二主成分基本和风速无明显关系. 第一和第二主成分和风速的拟合优度分别为0.68 和0.03. 第一成分和风速的拟合优度明显高于图6 的相干时间和归一化时延波形面积的拟合优度，能够提供更好的风速反演结果. 第二成分和风速无明显相关性，表现了极差的拟合优度. 这说明风速信息主要包含在第一主成分内. 图8 给出了同比风速和GEO 01 号卫星相干时间、归一化时延波形面积、PCA 及LUMV 反演的风速时间序列. 从图可知，相比于PCA 和LUMV 联合多颗卫星的多参数反演的风速较相干时间和归一化时延波形面积具有更光滑的风速反演结果. 由表1 可知利用第一主成分反演的风速的均方根误差为1.21 m/s，优于单星单参数反演结果. 利用3.2 节所述的LUMV 估计器对各特征参数独立反演的风速进行最优加权组合. 组合后的风速的均方根误差为1.15 m/s，获得了最佳的反演结果.

图8 反演风速和同比风速的时间序列Fig. 8 Time series of retrieved wind speed and in-situ wind speed

表1 反演风速的均方根误差Tab. 1 RMSE of retrieved wind speed m/s

5 结 论

相比于GPS，Galileo 和GLONASS 系统，北斗系统的GEO 卫星具备独有的优势，在岸基场景下可构成稳定观测几何，为长期连续的观测近岸海域风速提供了必要条件. 随着北斗三代的建成，B3I 信号已经在GEO 卫星上播放. 目前研究均集中在利用北斗GEO 卫星的B1I 信号反演地球表面的物理参数，而针对北斗B3I 信号的研究较少. 此外，由于北斗系统从东经58°～160°部署了5 颗GEO，使得直反射天线可同时接收多颗GEO 卫星信号，为多星多参数融合提供了条件. 本文针对北斗GEO 卫星B3I 信号开展了反演海面风速的研究，提出了PCA 和LUMV 的多星多参数海面风速反演算法. 试验结果表面相比于波形面积，在岸基场景反射信号的相干时间能够提供更好的反演性能. 相比于单星单参数反演，多星多参数方法的风速反演精度更高，而LUMV 获得了精度为1.15 m/s 的最佳风速反演结果. 值得注意的是北斗GEO 卫星同时播放B1I 和B3I 信号，联合双频点的反演将可能进一步提高测量精度.