星载GNSS-R海面风场观测载荷关键技术设计与验证

王延光，白照广，朱雪萍，王崇羽，韩琳

1. 空间电子信息技术研究院，西安 710100

2. 航天东方红卫星有限公司，北京 100094

1 引言

GNSS-R技术是一种被动式的遥感手段，由Martin-Neira在20世纪90年代提出[1]。经过多年的发展，该技术逐渐走向成熟并进入应用阶段。该技术已经广泛应用于海面高度测量[1-2]、海面风场反演[3]、海冰反演[4]、土壤湿度测量[5-6]、海面浮油探测等方面[7]。与传统的遥感手段相比，该方法有如下优点：1)采用无源探测方法，不需要发射设备；2)全球覆盖均匀，获取数据量大。把GNSS反射信号处理技术与低轨卫星技术相结合，充分发挥GNSS信号全球分布的特性，并利用低轨卫星全球扫描的特性，覆盖全球的海面风场测量系统已经走向工程应用。ESA和SSTL等机构都在致力于星载GNSS-R接收机的研制[8]，NASA发射成功的用于飓风观测的CYGNSS也是多颗联合星载接收机系统[9]。

与地基信号处理系统相比，卫星系统信号处理资源受限，包括数据传输带宽和星上信号处理系统的存储量及信号实时处理能力均受限于卫星平台能力。时域算法通过时域相关的方式实现反射信号处理[10]，处理结构简单，但是需要进行海量相关运算，复杂度大占用大量硬件资源，在卫星资源约束条件下无法实现多点反射信号实时处理；文献[11][12]直接使用频域相关算法，不考虑GNSS-R信号处理积分时间和频率覆盖范围的特殊性会因为码频率和载波频率失配造成DDM图失真。为了实现对所有落入反射天线波束的GNSS反射信号高精度实时信号处理，需要从系统设计和信号处理算法两个方面考虑，提高信号处理速度和精度。

本文充分考虑了GNSS-R信号处理在积分时间、信号动态在数据处理过程中对相关值有效相位的影响，以及星载实时信号处理系统的资源约束，改进了文献[11][12]的频域处理方法在积分过程中本地码相位偏移的缺陷，设计了基于短时循环相关与动态相位补偿的高速处理算法和信号处理系统。经过地面和在轨测试，系统处理性能和处理损失均满足设计要求，有效地支撑了中国第一颗GNSS-R系统的研制。对推动GNSS-R技术的研究与应用具有一定参考价值。

2 GNSS-R信号模型与风场观测原理

GNSS信号由载波、扩频码、导航电文调制生成，在信号发射端，第i颗GNSS卫星发射的信号可以表示为：

si(t)=aicos(2πfL)ci(t)Di(t)

(1)

式中：ai表示信号幅度，fL是载波频率，ci，Di表示导航电文。GNSS-R接收机接收到的反射信号是反射面所有面元上散射信号的和，可以表示为：

exp(j2π(f-fρ)t)d2ρ+n(t)

(2)

式中:Aρ为散射面元处的信号幅度；c为伪随机码；d为数据码；τρ和fρ分别为散射面元上的时间延迟和多普勒频率；n(t)为高斯白噪声。

GNSS 海面散射信号相关功率模型是时间延迟τρ、多普勒fρ的函数。散射功率在时延-多普勒域的分布被称为延迟多普勒图像(delay Doppler mapping，DDM)，是反射信号处理经典方式，如图1所示。该二维相关函数图像的时延、多普勒的相关值包含着海面波浪变化等重要的物理特性[3，13]。风速通过风应力对海面作用，风生成海表面波，这些表面波改变了海面粗糙度，而GNSS 海面散射信号相关功率后沿斜率与海面粗糙度相关[14]。图1为不同海况下GNSS 海面散射信号功率曲线，其后沿斜率与风速具有直接关系，海面风速越大，波形的后沿变化越趋于平缓，反之陡峭。通过精确测量GNSS 海面散射信号相关功率可以反演海面风速。

图1 时延多普勒二维相关功率和风速的关系Fig.1 Relationship between DDM and wind speed

由于GNSS卫星和接收机的几何分布离散，GNSS-R探测的幅宽不如传统卫星遥感载荷，只能形成若干分散的探测条带。为了获得更好的观测覆盖，信号处理系统需要实时完成对落入反射天线的所有反射信号进行处理。与地面系统高性能硬件不同，受限于卫星平台功率、数传能力、芯片处理能力和可靠性要求等因素，通过高效算法提高有限硬件资源的信号处理能力是GNSS-R载荷设计的必由之路。捕风一号卫星GNSS-R载荷及信号处理算法就是基于以上需求进行的设计。

3 星载GNSS-R载荷系统及信号处理算法设计

3.1 GNSS-R载荷系统组成

GNSS-R海面风场观测载荷主要任务是采集GNSS-R接收天线足印区域内可见GPS/BDS卫星镜面点的散射信号，并实时生成镜面反射点对应区域的DDM图，回传至地面。

星载GNSS-R载荷系统实现框图如图2所示，载荷主要由直达天线、反射天线、LNA、信号处理机组成。载荷接收GNSS卫星的直达和反射信号，经过滤波放大后在信号处理机完成直达和反射信号的处理。其中直达信号用于实现对本星的导航定位和反射信号参数估算，反射信号经处理后实时生成DDM图。为了获得较好的反射信号可见性，系统设计了双路反射天线倾斜安装。

图2 GNSS-R海面风场观测载荷系统框图Fig.2 Block diagram of GNSS-R sea surface wind field observation system

3.2 GNSS-R信号处理算法设计与性能分析

反射信号可被视为反射区域不同时延和不同多普勒信号的和，其生成DDM图的信号处理过程是本地再生信号和接收散射信号在时域和频域的二维相关值，处理过程包含大量乘加运算。

从反射信号处理时延-多普勒二维信号功率分布计算需求出发，最直接的方式是对AD原始采样信号进行不同频率的下变频，再通过时域相关的方式获得不同频点和不同时延的相关值[10]。如图3，该处理方式需要二维交叉遍历所有频点的码相位并进行长时间积分，运算量较大。

图3 DDM基本处理循环Fig.3 DDM basic processing loop

为了提高反射信号处理效率，可将反射信号处理与GNSS捕获算法相结合。GNSS-R反射处理，是对信号进行码相位和载波频率二维搜索。这和GNSS信号捕获的区别有两个方面。首先，由于反射信号功率较弱，为了满足观测数据信噪比要求，需要进行长时间的非相干累加。其次，由于反射信号处理关心的码相位和频率范围分别是以镜面反射点处信号时延和多普勒为中心±16 chip和±5 kHz，所以处理时无需对全码相位进行相关。

为了提高信号处理速度，实现实时5幅DDM图输出，系统设计中考虑了以下几个特点进行展开。

1)信号处理算法设计充分借鉴了GNSS循环相关捕获技术；

2)结合GNSS-R信号处理在时域、频域关注范围小的要求；

3)充分考虑由于信号频率不匹配造成本地码和接收码相位不匹配。

如图4，采用基于FFT的并行码相位搜索方式来提高运算效率。该方法通过FFT运算和IFFT运算获得信号s(n)和本地码c(n)的相关值。单次运算可以获得整个码周期的全部相关值。

图4 并行码相位相关算法Fig.4 Parallel code phase correlation algorithm

由于FFT运算的高效性，并行码相位算法和直接进行卷积相关的计算方法相比，显著降低了计算量。由于反射信号DDM图只关注镜面反射点对应信号时延附件(±16码片)的相关功率，并行码相位的计算方法计算了所有码相位上的相关结果，需要舍弃大部分计算结果，计算效率受到限制。

本文把通常用于GNSS长码捕获或剔除数据调制的DBZP相关算法和旋转变化技术引入到反射信号处理中，并通过相位偏移估计预先补偿由于本地码频率和接收码频率不匹配造成的码相位偏移问题，其基本思路是：

1)利用DBZP通过短时循环相关缩短傅里叶变换的数据长度，避免过度计算不需要的相关值。由于缩短了相干积分时间，等效为提高了相关值的采样率，可以通过旋转变化的方式获得出不同频率维度的相关值序列。

2)通过频率和接收频率误差估计，估计旋转变化相关值索引。

3)对2)获得的索引对应相关值进行旋转变换和累加，获得最终DDM图。

其具体流程如图5所示。

图5 DBZP 相关算法Fig.5 DBZP correlation algorithm

步骤1 数据预处理：

是对原始AD采样数据进行下变频至零中频(下变频中心频率为估算镜面反射点处反射信号的多普勒)，通过低通滤波器滤除和频分量，再对零中频数据进行重采样，目的是降低数据采样率，提高运算效率，在DDM信号处理中选择4倍码速率作为基带信号重采样频率，生成复数基带序列s(n)。

根据反射信号多普勒对应的码多普勒fD_code产生本地码序列c(n)。

(4)

式中:fD为载波多普勒，fc为扩频码频率。

以DDM图关注码相位范围(±16chip)为长度对s(n)和c(n)分块，单块长度为128个采样点。

步骤2 DBZP循环相关：

对s(n)分块后各块表示为blocks(1)～blocks(N+1), 对c(n)分块后获得blockc(1)～blockc(N)，为了DBZP计算数据采集分块总数比扩频码分块总数多1，N=T/(128fsamp)，T为总处理时间。其中fsamp为信号采样频率。通过[blocks(k),blocks(k+1)]和[blockc(k),O]，k=1～N，做循环相关。每次截取前128个相关值可获得相关值矩阵corr(m,n)m=1～N,n=0～127。具体流程如图5所示。

考查循环相关后的相关值矩阵corr(m,n),每行相关值代表以反射点时延相位为中心±64个采样点对应的相关值，每一列代表同一码相位连续采样时刻对应的相关值。

第K列向相关值可表示为[15]：

(5)

式中：sin(πfDTcoh)/(πfDTcoh)表示由于频率不一致造成的处理损耗；ej[2πfD(tn)θen]表示残留载波。

步骤3 生成相关值偏移索引：

考查DBZP获得的相关值矩阵corr(m,n)，如果恢复码频率和恢复载波频率一致，矩阵corr(m,n)中的每一列代表同一码相位在连续时间段内的相关值，对相关值矩阵进行列向相干和非相干累加可以获得对应时延的相关值，如果进行旋转后累加，可获得不同频率处的相关值。但是对corr(1:N,n)累加的时间长度和码多普勒相关。由于反射信号处理的总积分长度大于1 s，在此条件下即使一个很小的码频率不匹配也会造成相关值矩阵中第一行的相关值对应的码相位和最后一行相关值对应的码相位发生变化。如果直接对矩阵中的元素进行列向累加，会造成DDM处理结果与设计不符。

以BDS B1I信号为例，fcode=2.046 MHz，载波频率fL=1 561.098 MHz,fD_code=fD/763。意思是码相位偏移1个chip每1/fD_code秒。考虑到反射信号处理原始数据采样率为4fc，所以每1/4fD_code秒相位偏移1个采样点。

(6)

由于码多普勒的存在，对于corr(m,n)，corr(0,n)对应的相关值码相位和corr(k,n)相关值对应码相位存在偏移，具体偏差量为：

(7)

根据公式(8)，如果多普勒偏差较小或总积分时间足够小，codeshiftsamp(N)<1时码相位偏移可以忽略不计。以反射信号处理常用的积分时间设为1 s为例，当频率动态范围小于192 Hz时，则总相位偏移为1个采样点，可以忽略。

由于码频率失配造成的相位漂移现象如图6所示，此图设置的场景为频偏1 kHz，高亮斜线表示随着时间增加。相位偏移的过程，且随着处理码速率的增加倾斜度也增加。

图6 由码频率失配造成的码相位偏移仿真Fig.6 Demo of code phase shifting due to code Doppler mismatch

根据公式(8)计算结果，对不同积分时间和不同载波覆盖范围的码相位偏差修正补偿。在积分过程中，如果由于频偏造成的偏差小于1个采样点就进行列向累加。随着积分时间的增加，当相位偏移量大于一个采样点时就沿多普勒方向偏移一个采样点，形成倾斜累加位置索引。

步骤4 倾斜旋转累加:

(8)

沿着不同多普勒时信号码相位偏移计算出的相关值索引对corr(m,n)中的数据进行累加。累加过程包括128个连续相关值的相干累加和1 000次非相干累加，获得DDM相关值矩阵corrDDM(m,n)不同相位的相关值。

3.3 载荷运算负担分析

反射信号处理的运算量取决于DDM图需要覆盖的频率范围和时延范围，一般条件下单幅DDM图的时延维相关范围为正负16 chip，频率维覆盖范围为正负5 000 Hz，频率步进为500 Hz，所以单幅DDM图包含64×21个相关值点，每个相关值都由1 ms相干积分加1 000次非相干积分获得。

如果采用传统时域相关算法，单个1 ms(4 092个采样点)相干积分的运算量为4 092次乘加运算，所以单点的运算量为4 096×1 000；单幅DDM图的运算量为4 092×1 000×128×21=1.1×1010，所以总运算负责度为O(1.1×1010)。

使用传统循环相关的方式实现DDM处理，和时域相关的方式相比，1 ms(4 092点)采样数据通过两个4 096点的正FFT和一个IFFT实现，然后对IFFT结果进行非相干累加，一次获得4 096个码相位的相关值，4 092×128点的乘加运算被替换为3个4 096点的FFT运算器运算，由于FFT运算的复杂度为：

Nlog2N=4 096×12×3(N代表FFT长度)，总运算量为4 096×12×3 ×1 000×21=3.1×109，所以总运算负责度为O(3.1×109)。

使用DBZP循环相关配合带码相位偏移实时补偿的DDM算法(以128点分段长度为例)，计算量包含三部分。第一部分是分段补零循环相关，每次分段循环相干运算包含2次128点的FFT变换和1个128点的IFFT变换，对1 s内fsamp/64个分段独立计算。第二部分是对相关值矩阵列向旋转累加，同时进行相位偏移补偿，每个独立频率独立计算。第三部分是由于128点FFT不能覆盖全部关注码相位，需要将第一步和第二步操作更换相位后重复一次。其总的运算量为:

3.45×108

各算法的运算复杂度对比如表1所示。

表1 各算法运算负担表

可见改进算法的运算量与直接时域处理运算或传统循环相关相比，运算量分别减小到传统方法的1/31.9或1/8.99，具有明显优势。

4 GNSS-R载荷系统试验验证

4.1 地面验证

星载GNSS-R载荷系统地面测试验证使用海面风场模拟器进行，如图7所示。

图7 海面风场模拟器设备组成框图Fig.7 Equipment composition block diagram of sea wind field simulator

风场模拟数据根据GNSS-R观测几何数据集和气象样本数据集，结合GNSS-R载荷设计要求，进行海面风场气象样本的重采样形成，以GNSS-R观测几何集合气象样本数据集作为输入，完成GNSS-R等时延多普勒空间几何划分、海面风场海浪谱优选、Fresnel反射系数计算、KA-GO光学近似模型和Z-V电磁散射模型分析和模拟数据生成并存储。测试过程中，GNSS-R载荷直接接收模拟器提供的直达及反射射频信号，实时计算反射点和反射信号相位，并完成DDM图生成。

如图8所示，载荷绘制的DDM图与风场模拟器模拟的DDM图特性一致，风速反演结果与模拟值一致，误差为2.51 m/s。

图8 捕风一号卫星载荷地面测试结果Fig.8 Test results of BF-1 satellite

4.2 在轨验证

2019年6月5日，捕风一号卫星由长征11号运载火箭在距山东省海阳市200km海域圆满发射成功。卫星配套了本文描述的GNSS-R载荷设备。如图9所示，在轨测试期间设备成功完成星载DDM图绘制。

图9 捕风一号卫星在轨观测获取的DDM图像实例Fig.9 DDM acquired by BF-1 satellite on-orbit observation

为了验证捕风一号卫星GNSS-R载荷数据可用性，国家卫星气象中心与航天科技集团共同利用大洋浮标和再分析场完成了捕风一号L1级数据风场反演可用性评估，载荷数据实时数据产品成功率优于90%，评估表明系统时延计算精度为23 ns，镜像点位置计算精度为2 m，DDM数据满足2米/秒或10%的风速反演要求。至今，系统运行稳定，利用捕风一号卫星已经建立了全球海面GNSS-R信息观测能力。

地面和在轨试验表明，快速处理算法在提高处理速度的同时未对DDM处理结果产生明显影响，算法及信号处理系统可行。

5 结论

本文以捕风一号卫星为研究对象，介绍了GNSS-R反射信号处理载荷和改进的实时信号处理算法设计。算法结合了短时循环相关和时域旋转方法并通过实时相位偏移估算解决了由于本地码频率和实际码频率不一致造成的处理失真问题。该处理算法有效地提升了星载接收机信号处理能力,在相同条件下信号处理运算量减小到传统循环相关方法的1/8.99。系统设计有效地缓解了卫星平台在体积、质量、功耗和数传能力等方面的负担，满足实时信号处理的要求。GNSS-R载荷系统在轨测试表明本设计提供的DDM观测数据有效，数据质量满足应用需求。系统设计与验证为星载GNSS-R技术进化和业务化运行打下坚实基础。