風電機葉片微多普勒頻率求解及其影響因子

唐 波，陳 昊，黃 力，袁發庭，奉 彭，李耀偉

(三峽大學 電氣與新能源學院，湖北 宜昌 443002)

風電機旋轉葉片會對雷達入射電磁波產生微多普勒效應[1-4]，導致空間電磁環境變得復雜，嚴重影響雷達的探測性能[5-6]。從當前研究來看，獲取旋轉葉片產生的微多普勒頻率，進而在雷達側進行濾波是解決風電場對雷達臺站電磁干擾的關鍵步驟[7-8]。

當前，獲取葉片微多普勒頻率的方法有實驗測量[5,9-10]和數學仿真計算[1-4,11]兩類。雖然實驗測量是最直接獲取葉片微多普勒頻率的方法，但是由于實驗條件的制約以及外界環境的干擾，導致實驗測量結果效果不佳。因此，目前關于葉片微多普勒頻率的研究重點集中于采用數學仿真計算的方法。

現有的葉片微多普勒頻率的仿真計算方法可分為高頻近似算法[11]與散射點積分算法[1-4]兩類。其中，高頻近似算法基于準靜態的思想，通過獲取葉片動態雷達散射截面的時間序列，進而求解出葉片的微多普勒頻率。雖然該算法求解結果較為準確，但是需要對葉片每個時刻的運行姿態進行精確建模求解，從而導致計算量和復雜度過大而難以接受。散射點積分算法通過積分的方式，獲取葉片回波表達式，進而求解出葉片的微多普勒頻率。散射點積分算法因其計算速度快的優勢，已成為當前葉片微多普勒頻率求解的通用算法。

文獻[1]在研究風電場對航管一次雷達的干擾時，最早采用散射點積分算法求解了葉片的微多普勒頻率，但是僅考慮了雷達和葉片位于同一平面時的情況，沒有考慮微多普勒頻率的影響因子。后續研究中，文獻[2]將散射點積分算法應用于三維空間，采用短時傅里葉變換求解了雷達和葉片不在同一平面時的微多普勒頻率；文獻[12]同樣基于該方法，分析了方位角和俯仰角兩個影響因子對葉片微多普勒頻率的影響。但由于散射點積分算法忽視了散射點的客觀離散性，從而導致該算法獲取的微多普勒頻率不夠準確，僅適用于定性分析。

為此，本文避開散射點積分算法的假設缺陷，基于葉片微多普勒頻率的產生機理，引入電磁散射中心的概念，推導了葉片散射中心位置參數與微多普勒頻率的數學關系，從而實現了葉片微多普勒頻率準確求解，并結合實際風電機和雷達的運行工況，分析了各種影響因子對葉片微多普勒頻率的影響。

1 葉片微多普勒頻率及其散射點積分算法

1.1 風電機葉片的微多普勒頻率的產生機理

雷達在進行目標探測時，天線會在一定角度范圍內發射電磁波。鄰近雷達臺站的風電機不可避免地會出現在其探測路徑上。當入射電磁波照射到風電機時，葉片的旋轉運動會對入射電磁波進行周期性調制，從而產生了葉片微多普勒頻率。微多普勒頻率會改變雷達入射電磁波的頻率，導致空間電磁環境變得復雜，嚴重影響了雷達的探測性能，由此形成了風電機對雷達臺站的電磁干擾。從目前研究成果來看，解決干擾的關鍵在于獲取準確的葉片微多普勒頻率。

在工程實際中，由于風電機葉片的旋轉速度會根據風向和風速進行調整，雷達也會在不同的方位角、俯仰角下發射不同頻率的電磁波，從而導致旋轉葉片產生的微多普勒頻率存在較大的差別。因此，需要采用控制變量的方法，分別討論葉片轉速、入射電磁波頻率、方位角以及俯仰角對葉片微多普勒頻率的影響。

1.2 微多普勒頻率的散射點積分算法

當前，在葉片微多普勒頻率求解領域，通常采用的是散射點積分算法。該算法首先將葉片等間距分割成若干個薄片，并取每個薄片的中心作為等效散射點；然后根據雷達回波的基本方程，推導葉片上任一散射點的回波表達式，將該表達式沿葉片軸線進行積分，從而得到整個葉片的時域回波表達式；最后采用短時傅里葉變換法，對葉片時域回波進行時頻域分析，即可求解出葉片的微多普勒頻率。

上述散射點積分算法要求散射點沿葉片軸線分布，且默認了相鄰散射點之間具備數學連續性。但是，雷達在實際工作中，受制于工作帶寬和轉角，其分辨率為有限值，從而導致葉片相鄰散射點之間必定存在幾何間距。綜合以上分析可以看出，由于散射點積分算法過于理想化，僅適用于葉片微多普勒頻率定性分析。

2 基于散射中心的葉片微多普勒頻率求解

2.1 散射中心的概念

為實現旋轉葉片微多普勒頻率的準確求解，根據高頻電磁散射理論，在光學區，目標總的電磁散射可以認為是某些局部位置上電磁散射的合成，這些局部散射源被稱為散射中心[13-14]。依據這一機理，通過求解每個散射中心產生的微多普勒頻率，進而采用疊加原理，以實現整個葉片的微多普勒頻率求解。

2.2 葉片散射中心與微多普勒頻率的數學關系

如圖1所示的葉片旋轉模型，以葉片旋轉軸心為坐標原點O，垂直于葉片旋轉面的方向為x軸，平行葉片旋轉面的方向為y軸，豎直向上的方向為z軸，建立三維坐標系。假設葉片旋轉角速度為ω、旋轉方向為逆時針。雷達固定于空間點R處，其距離葉片旋轉軸心O的距離為R0.Pi是葉片上任意一個散射中心，其空間坐標為(0，yi，zi).當雷達發射機在俯仰角為θ、方位角為φ的方向發射頻率為fc的電磁波照射葉片時，由微多普勒的定義可知[15]，散射中心Pi產生的微多普勒頻率為：

(1)

式中：fDi(t)為t時刻葉片散射中心Pi產生的微多普勒頻率；fc為雷達發射機所發射的電磁波頻率；c為電磁波波速；vi(t)為散射中心Pi相對雷達R的瞬時旋轉線速度。

由于散射中心瞬時旋轉線速度vi(t)遠小于電磁波波速c，因此式(1)可近似為：

(2)

從式(2)可以看出，求解葉片微多普勒頻率的關鍵在于獲取瞬時旋轉線速度vi(t).

假設t時刻，葉片處于運行姿態1，此時散射中心Pi與雷達R之間的距離為RPi(t).經過時間間隔Δt后，葉片處于運行姿態2，此時散射中心與雷達之間的距離為RPi(t+Δt)，則該時間段內散射中心與雷達之間距離的變化量為ΔR=RPi(t+Δt)-RPi(t).當時間間隔Δt無限小時，可用該時間段內的平均速度代替散射中心Pi的瞬時旋轉線速度：

(3)

由式(3)可知，散射中心Pi的瞬時旋轉線速度vi(t)即為散射中心Pi與雷達R之間的距離RPi(t)對時間t的導數。因此葉片微多普勒頻率的求解就轉化為對散射中心Pi與雷達R之間的距離RPi(t)的求解。對于△ORPi，根據余弦定理可求得RPi(t)為：

(4)

式中：γi(t)為t時刻散射中心Pi和葉片軸心O的連線OPi與雷達視線之間的夾角。

實際工程中，由于雷達R與葉片通常滿足遠場條件，即R0?OPi.因此，式(4)可近似為：

(5)

由于式(5)中的參量cosγi(t)難以直接得到，因此需要利用連線OPi與連線OR之間的夾角公式進行求解：

cosγi(t)=sinφsinθcosβi(t)+cosθsinβi(t) .

(6)

式中，φ為連線OR在xoy面內的投影與x軸正方向的夾角，即雷達相對于葉片的方位角；θ為連線OR與z軸正方向的夾角，即雷達俯仰角；βi(t)為t時刻連線OPi與y軸正方向的夾角，βi(t)=βi0+ωt，βi0為連線OPi與y軸正方向的初始夾角。

將式(6)和式(5)代入式(3)，即可求解出散射中心Pi所產生的微多普勒頻率為：

(7)

由上述求解過程可知，散射中心Pi所產生的微多普勒頻率與其縱坐標和豎坐標密切相關，并且受葉片轉速、雷達入射電磁波頻率、方位角、俯仰角的影響。基于此，若能獲取葉片散射中心的位置參數，則可求解出整個旋轉葉片的微多普勒頻率，進而可分析各個影響因子與葉片微多普勒頻率的關系。

2.3 葉片散射中心位置參數的求解

為實現葉片散射中心位置參數的求解，首先需要尋求能夠準確描述葉片高頻電磁散射特性的散射中心數學模型。從現有研究來看，幾何繞射理論(geometrical theory of diffraction，GTD)散射中心數學模型被廣泛應用于雷達目標的散射中心求解[14]。根據GTD，當雷達在俯仰角為θ的方向發射頻率為f的電磁波照射葉片時，其后向電磁散射場可用以下散射中心數學模型表示為[14]：

(8)

式中，E(f,θ)為葉片的后向散射場；I為散射中心總個數；yi和zi分別為第i個散射中心的縱坐標和豎坐標，為待求參量；Ai和αi分別為第i個散射中心的散射強度和散射類型參數；f為電磁波的步進頻率，f0為電磁波的起始頻率，f=f0+m·Δf，m=0,1,…,M-1，其中Δf為頻率采樣間隔，M為頻率采樣點數；θ為俯仰角，θ=θ0+n·Δθ，n=0,1,…,N-1，其中θ0為雷達初始俯仰角，Δθ為俯仰角采樣間隔，N為俯仰角采樣點數；W(f,θ)為零均值二維復高斯白噪聲。

在實際應用中，散射中心數目I是未知的，需要利用葉片的后向散射場數據對其進行估算。從現有研究來看，目標散射中心的數目估算通常采用基于特征值取值范圍的蓋氏圓盤(gerschgorin disk estimator,GDE)法[13]。采用GDE法估計散射中心數目時，根據葉片后向散射場數據構成的協方差矩陣所具備的特點，即信號蓋氏圓盤半徑大而噪聲蓋氏圓盤半徑小，從而通過統計半徑大的蓋氏圓盤個數確定散射中心數目。

估算出散射中心數目I之后，即可對式(8)中的散射中心位置參數進行求解。由于式(8)中冪函數和指數函數同時存在，導致參數求解困難，因此需對式(8)進行化簡。當Δf/f0?1時，可用指數函數(1+αi·Δf/f0)m代替冪函數(f/f0)αi，從而式(8)可簡化為僅含指數函數的形式：

(9)

(10)

式(9)是空間譜估計的標準模型，可以采用空間譜估計算法[14]對式(10)中的P1i和P2i進行估計，進而求解出葉片散射中心的位置參數yi和zi.考慮到目標散射中心位置參數的求解技術已經非常成熟，如可用旋轉不變技術信號參數估計法[14]對其進行求解，具體方法在此不再贅述。

3 葉片微多普勒頻率的求解及其影響因子

3.1 葉片微多普勒頻率求解結果

選取金風77/1500型風電機，建立相同尺寸的葉片幾何模型進行微多普勒頻率算例分析，其幾何模型如圖2所示。該型風電機的葉片長度為35 m，并令第1支葉片的初始位置與y軸正方向平行，順時針分別為第2支葉片和第3支葉片。

圖2 風電機葉片幾何模型

根據現有的散射中心提取方法，葉片散射中心位置參數的求解結果如圖3所示。圖3中分布有33個散射中心，其中，第1支葉片上分布5個散射中心，第2支葉片上分布13個散射中心，第3支葉片上分布14個散射中心，以及在葉片軸心處分布有1個散射中心。造成第1支葉片散射中心數目少于其它兩支葉片散射中心數目的原因在于，第1支葉片的初始位置平行于雷達入射電磁波，該支葉片被雷達信號照射的有效寬度較小，導致雷達難以分辨出該支葉片。

圖3 葉片散射中心位置分布結果

設置雷達入射電磁波的頻率fc為3 GHz，方位角φ為90°，俯仰角θ為90°；葉片轉速RPM為15，對應的旋轉角速度ω為0.5π rad/s；觀測時間設置為10 s.將葉片33個散射中心的位置參數代入式(7)，可得到旋轉葉片微多普勒頻率的求解結果如圖4所示。

圖4 葉片微多普勒頻率的求解結果

從圖4所示的仿真結果可以看出，葉片微多普勒頻率包括1條零頻帶和32條正弦曲線。其中，零頻帶與文獻[2]的仿真結果相一致，是由于葉片軸心部位相對雷達沒有發生位置偏移而引起的。但是，相比文獻[2]中僅有1條正弦包絡曲線，圖4中出現了微多普勒頻率正弦曲線分層現象，共有33條微多普勒頻率正弦曲線。其中，第1支葉片有5條微多普勒頻率曲線，第2支葉片有13條微多普勒頻率曲線，第3支葉片有14條微多普勒頻率曲線。葉片微多普勒頻率具有的這一特點與圖3中葉片散射中心的數目相一致。此外，圖4中出現了正負微多普勒頻率，當葉片遠離雷達運動時，出現正微多普勒頻率，當葉片靠近雷達運動時，出現負微多普勒頻率。結合葉片的旋轉速度可知，當葉片與入射電磁波垂直時，微多普勒頻率正弦曲線取峰值。根據圖2所示的葉片初始位置可知，最先與入射電磁波垂直的是第2支葉片，然后是第1支，最后是第3支，因此圖4中微多普勒頻率正弦曲線最先出現峰值的是第2支葉片，隨后是第1支，最后是第3支。

3.2 葉片轉速對微多普勒頻率的影響

保持雷達入射電磁波的頻率為3 GHz、方位角為90°以及俯仰角為90°不變，分析葉片轉速對微多普勒頻率的影響。根據金風77/1500型風電機的技術參數可知，其RPM為9～17.3，對應的旋轉角速度為0.3πrad/s～0.58πrad/s.本文選取葉片RPM為9、12、17.3三種情況進行分析，得到旋轉葉片微多普勒頻率的結果如圖5所示。

從圖5(a)、5(b)以及5(c)中可以看出，其他參數一定的條件下，葉片轉速對微多普勒頻率有兩方面的影響。一是葉片轉速越快，微多普勒頻率的變化周期越短。其中，當葉片RPM為9時，微多普勒頻率變化周期約為6.67 s；當葉片RPM為12時，微多普勒頻率變化周期為5 s；當葉片RPM為17.3時，微多普勒頻率變化周期為3.47 s.二是葉片轉速越快，每個散射中心產生的微多普勒頻率越大。為了更有效地說明轉速對葉片微多普勒頻率的影響，選取最大微多普勒頻率進行分析。當葉片轉速RPM為9時，其最大微多普勒頻率約為662 Hz；當葉片RPM為12時，其最大微多普勒頻率約為883 Hz；當葉片RPM為17.3時，其最大微多普勒頻率約為1 273 Hz.從圖5(d)中可以看出，葉片最大微多普勒頻率與轉速呈現出嚴格的正比關系。

圖5 葉片轉速對微多普勒頻率的影響

3.3 電磁波頻率對微多普勒頻率的影響

控制葉片RPM為15、雷達入射電磁波方位角為90°以及俯仰角為90°不變，分析入射電磁波頻率對葉片微多普勒頻率的影響。選取入射電磁波頻率為1 GHz，2 GHz以及4 GHz三種情況進行分析，得到葉片微多普勒頻率的結果如圖6所示。

從圖6(a)、6(b)以及6(c)中可以看出，其他參數一定的條件下，雷達入射電磁波頻率越高，每個散射中心產生的微多普勒頻率越大。同樣地，選取最大微多普勒頻率分析雷達入射電磁波頻率對微多普勒頻率的影響。其中，當入射電磁波頻率為1 GHz時，其最大微多普勒頻率約為368 Hz；當入射電磁波頻率為2 GHz時，其最大微多普勒頻率約為736 Hz；當入射電磁波頻率為4 GHz時，其最大微多普勒頻率約為1 472 Hz.如圖6(d)所示，葉片最大微多普勒頻率與入射電磁波頻率同樣呈現出嚴格的正比關系。

圖6 入射電磁波頻率對微多普勒頻率的影響

3.4 方位角對微多普勒頻率的影響

控制葉片RPM為15、雷達入射電磁波頻率為3 GHz以及俯仰角為90°不變，分析入射電磁波方位角對葉片微多普勒頻率的影響。選取入射電磁波方位角為0°、30°以及60°三種情況進行分析，得到葉片微多普勒頻率的結果如圖7所示。

圖7 入射電磁波方位角對微多普勒頻率的影響

從圖7(a)、7(b)以及7(c)中可以看出，其他參數一定的條件下，在0°～90°范圍內，雷達入射電磁波方位角越大，每個散射中心產生的微多普勒頻率越大。同樣地，選取最大微多普勒頻率分析雷達入射電磁波對微多普勒頻率的影響。其中，當入射電磁波方位角為0°時，其微多普勒頻率為0 Hz；當入射電磁波方位角為30°時，其最大微多普勒頻率約為552 Hz；當入射電磁波方位角為60°時，其最大微多普勒頻率約為956 Hz.值得注意的是，微多普勒頻率與方位角之間并不是正比關系，如圖7(d)所示，葉片最大微多普勒頻率與入射電磁波方位角呈現出正弦變化關系。具體表現為，當方位角為0°和180°時，微多普勒頻率為0，即入射電磁波垂直照射葉片平面時，不存在微多普勒效應。當方位角為90°和270°時，微多普勒頻率達到最大值，即入射電磁波平行葉片平面時，微多普勒效應最嚴重。

3.5 俯仰角對微多普勒頻率的影響

控制葉片RPM為15、雷達入射電磁波頻率為3 GHz以及方位角為90°不變，分析入射電磁波俯仰角對葉片微多普勒頻率的影響。選取入射電磁波俯仰角為0°、45°、135°以及180°四種情況進行分析，得到葉片微多普勒頻率的結果如圖8所示。

圖8 入射電磁波俯仰角對微多普勒頻率的影響

從圖8(a)、8(b)、8(c)以及8(d)中可以看出，其他參數一定的條件下，雷達入射電磁波俯仰角不影響葉片微多普勒頻率的大小，但是改變了每支葉片微多普勒頻率出現峰值的時刻。其原因在于，改變入射電磁波俯仰角時，三支葉片與入射電磁波垂直的先后順序放生了變化，從而影響了每支葉片微多普勒頻率出現峰值的時刻。

4 結論

1) 為準確獲取風電機旋轉葉片產生的微多普勒頻率，基于高頻電磁散射理論，推導了葉片散射中心位置參數與微多普勒頻率的數學關系，從而提出了基于散射中心的葉片微多普勒頻率求解方法，為葉片微多普勒頻率特性分析提供了新的方法。

2) 獲取了葉片微多普勒頻率隨葉片轉速、雷達入射電磁波頻率、方位角以及俯仰角的變化規律，為后續研究風電機雜波分量的剔除技術以及風電機陣列條件下的微多普勒頻率特性奠定了理論基礎。