一種改進的補償因子區域不變CS算法及Simulink模型建模

宋栩潮，周 鵬，張振華，肖忠源

（1 中國石油大學（華東）海洋與空間信息學院 青島 266580；2 北京遙測技術研究所 北京 100076）

引言

傳統實時SAR（Synthetic Aperture Radar，合成孔徑雷達）成像系統架構大都采用FPGA（Field-Programmable Gate Array，現場可編程門陣列）與DSP（Digital Signal Processor，數字信號處理器）相結合的架構[1?3]。其中，FPGA主要用于實現數據采集和預處理，DSP 主要用于實現SAR 成像算法[4]。近年來，隨著FPGA 器件的快速發展，以及其在星上特殊環境下具有優于DSP 的抗單粒子翻轉能力，采用FPGA 逐步替代DSP 來實現SAR 成像算法，成為星上SAR 實時成像系統的重要研究方向之一[5?7]。

然而，算法設計人員在進行FPGA 程序開發時常會受到以下限制[8?11]：①算法設計人員通常不熟悉硬件描述語言，導致開發周期長；②FPGA的并行程序設計與基于Matlab 或DSP 的串行程序設計有很大的差異。而利用Simulink 工具實現FPGA 快速開發的方式則較好地解決了上述問題[12]。算法設計人員利用Matlab 驗證算法的正確性后，可利用Simulink 搭建算法模型并自動轉換為FPGA 代碼，大大降低了FPGA 開發的難度、提高了開發效率。另外，當算法需要改進時，只需要在Simulink模型對應的模塊上進行修改，通過自動生成硬件語言即可實現算法的改進，靈活性較強。

CS 算法是合成孔徑雷達實時成像處理中常采用的一種算法。在經典CS 算法的處理過程中，涉及到相位補償、距離向補償和方位向補償這三種補償因子的計算。這些補償因子涉及到正余弦和開根號等超越運算，采用FPGA 實現這些計算需要較長的運算時間和很大的存儲空間。若對于數據矩陣中的每條數據都要更新補償因子，則嚴重影響SAR 成像處理的實時性。針對上述問題，文獻[13]提出一種具有補償因子區域不變特點的CS 算法，其特點是區域內的各條數據采用固定的補償因子。該方法在少量損失了三種補償因子精度的情況下大幅減少了計算量，提高了處理效率。但在該方法中，區域內的統一補償因子采用的是區域內第一條數據對應的補償因子，會造成相位補償誤差不均勻的問題。

針對上述問題，本文對文獻[13]中的方法進行了改進，將區域內的統一補償因子改用該區域內所有單元頻率平均值對應的補償因子進行代替，使得采用統一補償因子造成的相位誤差在區域內更加均勻。顯然，改進算法與原有算法的計算量相當。但由于改進后算法對應的相位誤差在區域內更加均勻，使得改進算法的成像質量得到提升。為此，首先開展了Matlab 仿真工作，驗證了改進后算法相比于原有算法在成像指標方面的提升。接下來，搭建了用于生成改進算法中三種補償因子的Simulink模型。之后，將Simulink模型輸出的代碼加載到Vivado 軟件中，通過對比Vivado 軟件輸出的補償因子和Matlab 軟件輸出的精確補償因子，驗證了所搭建Simulink模型的高精度。

1 補償因子區域不變CS 算法的改進及Matlab驗證

1.1 經典CS算法原理

雷達發射信號一般為LFM（Linear Frequency Modulated，線性調頻）信號。對去載頻后的點目標回波信號進行方位向FFT（Fast Fourier Transform，快速傅里葉變換），可得信號在距離多普勒域的表達式為[14]：

式中，τ表示距離向快時間，fη表示方位向頻率，r表示雷達與點目標的距離，A表示調幅項，包括目標后向散射系數、距離向天線方向圖和方位向天線方向圖等因素的影響，Ks表示距離信號多普勒域調頻率因子，c表示光速，λ表示雷達發射脈沖信號波長，v表示平臺飛行速度，Cs表示距離彎曲因子。

為使在不同距離上的移動軌跡都與在參考距離上的移動軌跡一致，需對式(1)乘以如下CS 相位補償因子：

式中，rref表示參考距離。

對乘以相位補償因子后的信號進行距離向FFT，將信號轉換至二維頻域，然后將信號二維頻域表達式乘以距離向補償因子，完成距離壓縮、二次距離壓縮和距離徙動校正。距離向補償因子的表達式為：

式中，fτ表示距離向頻率。

將二維頻域信號與式(5)相乘后進行距離向IFFT，再乘以方位向補償因子后完成方位壓縮和殘余相位補償。方位向補償因子的表達式為：

再將乘以方位向補償因子后所得的信號進行方位向IFFT（Inverse Fast Fourier Transform，快速傅里葉逆變換），即可得到目標圖像。

為方便在FPGA 上實時生成三種相位因子，需要對fη、fτ進行量化處理。將fη等分為Nη份，將fτ等分為Nr份（Nη和Nr分別表示方位向和距離向點數），即

式中，fs表示雷達信號采樣頻率，fPRF表示脈沖重復頻率。

CS算法的成像質量與式(4)、式(5)和式(6)表示的三種補償因子的精度有直接關系，因此在實際應用中，補償因子計算模塊的精度非常重要[15]。

1.2 補償因子區域不變CS算法的原理

CS 算法成像過程中含有三種補償因子的計算，且每種補償因子都需要與對應的數據矩陣進行復乘。同時CS 成像算法中補償因子的計算比較復雜，包括正余弦和開根號等多種超越運算，采用FPGA 實現需要較長的運算時間和很大的存儲空間，對SAR 成像處理系統的實時性產生較大影響。

針對上述問題，文獻[13]中提出一種具有補償因子區域不變特點的CS 算法，即選定一個更新步長?n=2n（n為正整數），每個更新步長內只計算第一個距離向或方位向的補償因子，將后幾個距離向或方位向的補償因子用其替換。該文中指出：?n越大，在區域內采用統一的補償因子帶來的誤差也越大；當?n不超過16 時，三種誤差均不超過6×10-5。該算法使得三種補償因子的計算量大幅降低，從而有效減少補償因子的計算時間及其占用的存儲器資源。該算法的原理如圖1所示。

圖1 補償因子區域不變CS算法的處理流程圖Fig.1 Flow chart of a chirp scaling algorithm featured by invariance of compensation factors in a region

對于補償因子區域不變CS 算法，在一個更新步長內，其統一補償因子滿足如下特點：

1.3 改進的補償因子區域不變CS算法

在1.2 節所述的現有算法中，對于更新步長內各單元的統一補償因子，使用的是第一個單元對應的補償因子。在該算法中，對于更新步長內首行/列的相位誤差為0，但后續的每個距離向和方位向單元對應的誤差會逐漸增加，導致補償后的相位誤差在區域內分布不均勻。

針對上述問題，本文對文獻[13]中的算法進行了改進，將一個更新步長內的統一補償因子，采用區域內所有單元頻率平均值對應的補償因子進行代替。相比于改進前的算法，在更新步長對應的區域內，中間單元對應的誤差最小，第一個和最后一個單元對應的誤差最大，補償結果的相位誤差在區域內更加均勻。

對于改進后的補償因子區域不變CS 算法，在一個更新步長內，其統一補償因子對應的頻率值滿足如下特點：

在FPGA 實現CS 成像算法的過程中，會頻繁地從存儲器中讀取提前存儲好的補償因子，影響處理效率。經典CS 算法在進行補償因子的復乘操作時，每進行一次復乘都需要給存儲器發送地址來讀取該地址下的數據，而改進后的補償因子區域不變CS 算法只需要每隔一個更新步長取出一次數據即可，大大降低了進行復乘操作時的存儲器讀寫頻率，同時節約了存儲器資源，相比于經典CS 算法更易于FPGA 實現。在實際應用時，只要在成像之前完成所有補償因子的計算并存入存儲器中即可，不會占用成像時間。改進算法相比于文獻[13]中的算法，計算量相當。由于最終的補償結果相比于文獻[7]中算法的補償結果更加均勻，成像質量將有一定提升。下面對改進算法中3種補償因子的相位誤差進行量化分析。

1.4 改進算法中3種補償因子的相位誤差分析

式(4)～式(6)中的Ks和Cs僅隨方位向頻率變化，而雷達有效速度v和方位向頻率fη在距離向不發生變化，因此v和方位向頻率fη不受更新步長的影響；而τ和r是隨著距離單元變化的，即：

式中，R0表示最近斜距。在一個更新步長內，相位因子的最大相位誤差為：

下面對三種補償因子的相位誤差進行量化分析。

1.4.1 CS相位補償因子的相位誤差

由式(4)、式(15)可得，第n個距離單元內的CS相位補償因子的相位可以表示為：

在一個更新步長內，對應的最大相位誤差?φ1為：

下面通過仿真，對兩種算法的最大相位誤差進行比較。仿真參數設置為：雷達發射信號頻率為5.3×109Hz，場景中心斜距為2 000 m，雷達有效速度為150 m/s，距離向調頻率為2×1013Hz/s，距離向采樣率為6×107Hz，方位向采樣率為200 Hz，距離向采樣點數為640，方位向采樣點數為2 048。兩種算法最大相位誤差隨更新步長變化的仿真結果見表1?？梢钥闯?，改進算法的最大相位誤差小于原始算法的最大相位誤差。

表1 CS相位因子最大相位誤差隨更新步長的變化Table 1 Variation of maximum phase error of chirp scaling phase factor with the update step size

1.4.2 距離向補償因子的相位誤差

由式(5)可得，第n個方位單元內的距離向補償因子的相位可以表示為：

在一個更新步長內，對應的最大相位誤差?φ2為：

利用與前文相同的參數進行了仿真，兩種算法最大相位誤差隨更新步長變化的仿真結果見表2?？梢钥闯?，改進算法的最大相位誤差小于原始算法的最大相位誤差。

表2 距離向補償因子最大相位誤差隨更新步長的變化Table 2 Variation of the maximum phase error of the range compensation factor with the update step size

1.4.3 方位向補償因子的相位誤差

由式(6)、式(16)可得，第n個距離單元內的方位向補償因子的相位可以表示為：

利用與前文相同的參數進行了仿真，兩種算法最大相位誤差隨更新步長變化的仿真結果見表3?？梢钥闯?，改進算法的最大相位誤差小于原始算法的最大相位誤差。

表3 方位向因子最大相位誤差隨更新步長的變化Table 3 Variation of maximum phase error of the azimuth compensation factor with the update step size

1.5 三種CS算法的Matlab驗證

對經典CS算法、補償因子區域不變CS算法和改進的補償因子區域不變CS 算法，采用Matlab 語言編制程序，分別對點目標仿真回波和實測數據進行了成像處理。其中，對點目標的成像指標采用距離向和方位向峰值旁瓣比、積分旁瓣比和分辨率進行比較，對實測數據的成像結果采用均方根誤差進行比較。由于FPGA 本身進行浮點數運算會存在誤差，更新步長的值不宜偏大，因此本文中假設更新步長等于4。

1.5.1 點目標仿真成像

首先，以更新步長等于4的情況為例，分別采用經典CS算法、補償因子區域不變CS算法以及改進的補償因子區域不變CS 算法，對點目標仿真回波數據進行了成像處理，得到詳細的成像結果及成像指標。在處理原始數據時，在方位向加矩形窗以對天線雙程方向圖的數據范圍進行限制。仿真參數設置如下：發射信號脈寬為2.5×10-6s，雷達發射信號頻率為5.3×109Hz，場景中心斜距為2 000 m，雷達有效速度為150 m/s，距離向調頻率為2×1013Hz/s，距離向采樣率為6×107Hz，方位向采樣率為200 Hz，距離向采樣點數為640，方位向采樣點數為2 048，波束方向為正側視。成像過程中三種算法均未進行加窗處理。圖2、圖3 和圖4 分別給出了三種算法的點目標成像結果、距離向剖面圖和方位向剖面圖。觀察圖2、圖3 和圖4 中的剖面圖與成像結果圖可以看出，三種算法對點目標的成像質量都很高。

圖2 經典CS算法的點目標成像結果Fig.2 Imaging results of a point target by the traditional CS algorithm

圖3 補償因子區域不變CS算法的點目標成像結果Fig.3 Imaging results of a point target by the CS algorithm featured by invariance of compensation factors in a region

圖4 改進的補償因子區域不變CS算法的點目標成像結果Fig.4 Imaging results of a point target by the improved CS algorithm featured by invariance of compensation factors in a region

為便于對三種算法的性能進行定量比較，分別計算了距離向和方位向的峰值旁瓣比、積分旁瓣比和分辨率，列于表4和表5。

表4 點目標仿真結果的距離向成像性能指標Table 4 Imaging performance indexes of a point target's result in range direction

表5 點目標仿真結果的方位向成像性能指標Table 5 Imaging performance indexes of a point target's result in azimuth direction

為了比較改進前、后兩種算法在更新步長較大情況下的成像質量，分別采用改進前、后兩種補償因子區域不變CS 算法，在更新步長等于8、16 的情況下對點目標進行了成像處理，成像結果如圖5和圖6所示，成像指標列于表6和表7。仿真參數設置與更新步長等于4的情況一致。

圖5 改進前后兩種補償因子區域不變CS算法的點目標成像結果（?n=8）Fig.5 Imaging results of point targets with original and improved compensation factor region invariant CS algorithm（?n=8）

圖6 改進前后兩種補償因子區域不變CS算法的點目標成像結果（?n=16）Fig.6 Imaging results of point targets with original and improved compensation factor region invariant CS algorithm（?n=16）

對比圖2、圖3、圖4 和表4、表5 中數據可發現：三種算法的成像質量較為接近。其中，經典CS 算法的成像指標最佳。但需說明的是，盡管經典CS 算法的指標最佳，但與另兩種算法的指標相比提升有限，而其付出的代價是運算復雜度高。相比于補償因子區域不變CS 算法，改進后算法的成像性能指標有提升。

對比表6、表7 中數據可以看出，隨著更新步長的增大，改進后算法的成像性能指標相比改進前的指標，其提升會更加明顯，進一步體現了在數據量大、硬件資源有限進而需要選取更大的更新步長的情況下，改進后算法相比于原始算法有更大優越性。

表6 點目標仿真結果的成像性能指標（?n=8）Table 6 Imaging performance performance of a point target's result（?n=8）

表7 點目標仿真結果的成像性能指標（?n=16）Table 7 Imaging performance of a point target's result（?n=16）

1.5.2 實測數據成像

實測數據采用一組來源于RADARSAT-1 衛星的某海岸線周邊區域的回波數據，令更新步長等于4，得到經典CS 算法、補償因子區域不變CS 算法和改進的補償因子區域不變CS 算法的成像結果如圖7所示。

圖7 3種算法對RADARSAT-1原始數據的成像結果Fig.7 Imaging results of three algorithms on raw data of RADARSAT-1

以經典CS 算法的成像結果為基準值，分別以補償因子區域不變CS 算法和改進的補償因子區域不變CS 算法的成像結果作為對比值，計算出兩種算法的均方根誤差指標。從表8數據可以看出，相比于改進前算法，改進后算法的成像性能指標得到提升。

表8 補償因子區域不變CS算法和改進的補償因子區域不變CS算法的均方根誤差Table 8 Mean root squares of the previous and the improved CS algorithms featured by invariance of compensation factors in a region

2 改進算法中補償因子計算模塊的結構設計和Simulink模型

2.1 補償因子計算模塊的結構設計

本文采用乒乓工作機制對改進算法中用于計算補償因子的模塊進行了設計，其信號實時處理能力是單個補償因子模塊的兩倍。當有數據輸入時，該結構采用了乒乓處理原理，將數據分別輸入通道1 和通道2 中，在輸出端對處理結果進行乒乓接收。圖8 和圖9 分別給出了所設計的補償因子計算模塊的乒乓結構示意圖和時序圖。

圖8 補償因子計算模塊的結構示意圖Fig.8 Structure diagram of the designed module for calculating compensation factor

圖9 補償因子計算模塊的時序圖Fig.9 Timing diagram of the designed module for calculating compensation factor

2.2 補償因子計算模塊的Simulink模型搭建

所搭建的三種補償因子計算模塊的Simulink模型分別如圖10(a)～圖10(c)所示。為了節約硬件資源，Simulink模型中所有輸入、輸出和中間變量的數據格式均采用32 bit 的單精度浮點數，數據格式中有7位有效數字。

圖10 三種補償因子計算模塊的Simulink模型圖Fig.10 Three Simulink models for calculating compensation factors

需注意的是，在補償因子的計算過程中，某幾個中間變量相乘后的值會非常大，使得整數部分占有較多的有效數字個數，從而導致產生較大的計算誤差。因此，為減小上述誤差來源，需對模塊中的部分組成進行專門的設計。下面以方位向補償因子計算模塊的部分組成為例進行說明。

在方位向補償因子計算模塊中，先提出了2π作為指數項的公因子，使得最終補償因子的計算結果僅取決于輸入數據的小數部分，再對計算過程中兩個帶有小數的中間變量相乘的項采用如圖11 所示的子模塊代替。圖11 中，fix 表示向0 靠攏取整，a、b表示輸入的被乘數，c表示乘數。通過該子模塊的專門設計，可使得兩個較大的含小數中間變量相乘后，整數部分占有的有效數字個數減少，從而提高最終計算結果的精度。

圖11 用于計算小數部分的子模塊Fig.11 Sub module for calculating decimal part

3 實驗結果

仿真參數設置如下：脈沖重復頻率為1 400 Hz，信號采樣頻率為3×107Hz，發射信號中心頻率為109Hz，發射信號距離向調頻斜率為2×1012Hz/s，雷達有效運動速度為1 000 m/s，距離向點數為1 024，方位向點數為1 024，方位向信號的多普勒調頻率為-8 357 Hz/s。

根據上述仿真參數，執行如圖8 所示的三個Simulink 模塊并輸出FPGA 代碼。將Simulink 模型輸出的代碼加載到Vivado軟件中進行執行。圖12給出了Vivado軟件輸出的三種補償因子的部分結果。

圖12 Vivado軟件的部分仿真結果Fig.12 Partial simulation results by Vivado software

利用Vivado 軟件將補償因子的計算結果存入txt 文件，然后由Matlab 讀取并將結果轉換為十進制小數存入.mat 文件中，每個補償因子的實部和虛部分別構成一個矢量。編寫改進的補償因子區域不變CS 算法的Matlab 程序，輸出正確的補償因子仿真結果，每個補償因子的實部和虛部分別構成一個矢量。將Vivado 輸出的補償因子與由Matlab 輸出的精確的補償因子進行相減，圖13 給出了三種補償因子實部和虛部的誤差。

圖13 Vivado輸出結果相比于Matlab輸出結果的誤差Fig.13 Error of Vivado output result compared with Matlab output result

從圖13 可以看出，CS 相位補償因子的誤差在±(2×10-3)以內，距離向和方位向補償因子的誤差在±0.01 以內，誤差均較小，驗證了本文所搭建Simulink 模型的準確性。其中，CS 相位補償因子和距離向補償因子的誤差曲線圖具有明顯對稱性，且越靠近成像場景中心的補償因子對應的誤差越小，體現了本文所提出改進算法的特點。

4 結束語

對一種具有補償因子區域不變特點的CS 算法進行了改進，對于更新步長內各單元的統一補償因子，采用該區域內所有單元頻率平均值對應的補償因子代替原來使用的第一個單元對應的補償因子，使得采用統一補償因子造成的相位誤差在區域內更加均勻。原算法與改進后的算法計算量相當。經過Matlab 仿真驗證，改進算法的成像指標得到提升，且更新步長越大，改進算法的性能提升越明顯。為提高FPGA 代碼的生成效率，采用Simulink工具搭建了用于生成改進算法中的三種補償因子模型。將Simulink 模型輸出的代碼加載到Vivado 軟件中，通過對比Vivado 輸出的補償因子與由Matlab 輸出的精確補償因子，驗證了所搭建Simulink 模型的高精度。本文所研究的利用Simu‐link 建模生成FPGA 代碼的方法具有較好的通用性，可推廣應用至SAR 實時成像或其他電子信息領域，提高FPGA代碼的生成效率。