星載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技術進化和業務化運行打下堅實基礎。