李一越，胡姝玲

(北京航空航天大學 儀器科學與光電工程學院，北京 100191)

1 引 言

大氣湍流是大氣運動的一種形式，不同于由規(guī)則的大氣分子運動形成的層流，大氣湍流屬于不規(guī)則運動，且運動的劇烈程度通常大于層流。在地表附近，由于受地形起伏、熱量交換等的影響，近地表的大氣運動往往存在強度不等的湍流。這些湍流會導致光信號在傳播過程中質量下降，具體表現(xiàn)為光束漂移、光束擴展、光強閃爍及到達角起伏等。

在風力發(fā)電領域中，湍流作為一種常見的大氣效應，對風電系統(tǒng)的危害主要包括風能利用率下降、發(fā)電量降低、風機軸疲勞及載荷增大等[1]。如新疆大型陸上風電廠中[2]，由于湍流特別是風力發(fā)電機尾流的影響范圍達到3倍葉片直徑以上，即300 m以上，由此導致的最大風速損失率為50%。由于傳統(tǒng)激光測風雷達無法將風機尾流同當?shù)刂行娘L速區(qū)分開，使后方風機在受到干擾時無法做出正確調整，因此年平均風能損失率約為10%。若不能將尾流及時分離出來，則需要根據(jù)當?shù)仫L場情況適當增加風力發(fā)電機間距離，這同樣會造成額外的風能浪費。同時，如果能夠及時觀測尾流等湍流對風場的影響，則可以及時調節(jié)不同風機的參數(shù)，保證系統(tǒng)整體的工作效率最高。

在風力發(fā)電機機組中，由于尾流問題導致不同位置風機的實際可利用風速的大小和方向均有所區(qū)別。若在風速分析中將湍流與中心風速進行混合分析，則所獲得的風速穩(wěn)定性低，同時需要部署多臺測風雷達用于整個風機組陣列的風場測量，以保證綜合利用率最大。但是這種做法會明顯提高系統(tǒng)的成本，同時降低系統(tǒng)的穩(wěn)定性(由于需要多臺測風雷達同時保持正常工作狀態(tài))。因此，需要對湍流進行研究并設計一種能在風機組陣列的任意或部分位置，在獲得先驗的湍流信息后，能夠根據(jù)即時風速信息，對整個陣列風場進行估計的方法。

在20世紀70年代至90年代，激光雷達成為大氣湍流探測領域的一種重要手段，以美國海洋大氣中心和俄羅斯大氣環(huán)境中心等一部分組織和機構對大氣湍流的性質和特點開展了一系列的相關研究[3]。在湍流測量分析領域的相關研究中，2008年Smalikho等[4]利用徑向速度譜來確定湍流的耗散率。2016年Matsuo等[5]利用神經(jīng)網(wǎng)絡對湍流信號進行分析，克服了傳統(tǒng)VPP、VAD算法在湍流分析領域的缺陷。徐春鳳等[6]針對湍流大氣中激光光強的閃爍特性進行了研究。2017年，Zhai 等[7]利用激光雷達對由于地面高度起伏和摩擦造成的近地表湍流進行了研究。

本文基于風力發(fā)電機組中常見的近地表湍流和渦流進行研究，在結合現(xiàn)有non-Kolmogorov湍流模型對湍流的脈動風速特點進行分析后，提出一種利用獨立分量分析(Independent Component Analysis，ICA)對常見的湍流大氣中回波信號的相位噪聲譜進行分離的方法。

2 激光測風雷達的設計參數(shù)與湍流

頻譜特點

2.1 激光雷達的設計參數(shù)

在風力發(fā)電機組中，激光測風雷達往往針對風機前方50～80 m處的近場風速信息進行探測。在這一探測范圍內，采用連續(xù)激光雷達可以獲得更高的探測精度和處理速度。本系統(tǒng)中采用的是由窄線寬半導體激光器作為激光光源，經(jīng)放大器后由收發(fā)一體式透鏡組雷達天線發(fā)射并在目標內獲得粒子的后向散射信號，再經(jīng)過Mach-Zehnder干涉系統(tǒng)將其多普勒頻移信息提取并進行處理，最終獲得相關的風速解算信號，其原理圖及設計參數(shù)如圖1所示。

圖1 激光測風雷達原理Fig.1 Block diagram of wind lidar system

圖中，LD為半導體激光器，ISO為光隔離器，EDFA為激光放大器，AOM為聲光移頻器，circulator為環(huán)形器，antenna為雷達天線，PD為光電探測器，ESA為電頻譜儀。在測量精度方面，雷達系統(tǒng)采用定時轉盤速度校準方法及回光信噪比監(jiān)測方法，確保測量精度在0.5 m/s，系統(tǒng)穩(wěn)定時最優(yōu)精度為0.1 m/s。

表1 激光測風雷達的設計參數(shù)Tab.1 Parameter designed for windlidar

2.2 風機附近風場湍流模型

通常情況下在風力發(fā)電機組之間，以尾流為主導的混合湍流中往往包含一些較為復雜的氣流擾動因素，如地形湍流等。大多數(shù)情況下，隨著氣象條件的穩(wěn)定，尾流在湍流中的影響會愈發(fā)明顯，因此在以尾流為主的大氣湍流的研究中，常用基于大氣凍結湍流假設條件下的non-Kolmogorov模型對湍流特性進行建模。其中，大氣凍結湍流假設在時間尺度不超過100 ms的前提條件下可以認為是有效的[8]，其過程表示為：

φ(z，t)=φ(z-ωtd，t-td)，

(1)

其中：td為系統(tǒng)延遲，ω為橫向風。這一過程表明，在特定的系統(tǒng)時延下，湍流風場中的某一點空間特性可以認為保持不變。在這一假設的基礎上，可以對湍流風場采用non-Kolmogorov模型來描述它的空間分布特性。這一模型具有可變的大氣指數(shù)α值，且其功率譜函數(shù)包括non-Kolmogorov、廣義von Karman譜等[9]。

其功率譜密度函數(shù)如下：

(2)

由式(2)建立湍流信號仿真模型，并得到如圖 所示的仿真結果，通過調節(jié)折射率結構常數(shù)及湍流尺寸可以獲得不同的湍流信號分布情況?？梢钥闯?，隨著湍流頻率分布范圍的擴大以及湍流整體功率的提高，信號的分布范圍逐漸擴大，同時信噪比下降，頻譜峰值穩(wěn)定程度下降，最終被噪聲淹沒。

(a)低湍流強度頻譜 (a)Moderate turbulence spectrum

(b)中湍流強度頻譜 (b)Medium turbulence spectrum

(c)高湍流強度頻譜 (c)Severe turbulence spectrum圖2 不同湍流強度的仿真頻譜Fig.2 Simulated power spectra in different turbulence levels

圖2表明當湍流增強時，系統(tǒng)的回光信號由于被湍流信道調制而出現(xiàn)整體質量下降。但同時，目標區(qū)域附近的湍流本身也是粒子運動所產(chǎn)生的真實后向散射信號，湍流的實時強度及其變化規(guī)律也是風場所需要的重要數(shù)據(jù)。因此，針對湍流信號的功率譜提取方法也顯得更加重要。

3 獨立分量分析模型的構建

獨立分量分析的概念可以追溯到20世紀90年代，由法國學者Jutten等首次提出了獨立分量分析(Independent Component Analysis，ICA)概念，該方法在90年代中期獲得國際信號處理界的廣泛關注。

假設利用3個探測器記錄信號xi(t)(i=1，2，3)為3個信號來源si(t)(i=1，2，3)的一個線性加權[11]，則它可用矩陣形式表示為：

(3)

其中hij(i，j=1，2，3)為信道權系數(shù)，其矩陣組合為信道模型。當信號來源增多時，可以利用更多的探測器對信號進行探測與辨識，同樣用矩陣形式表示，則有：

(4)

其展開式為：

(5)

這種多源信號混合分析問題被稱為“雞尾酒會”問題，如圖3所示。在忽略時延、非線性效應等因素的影響后，系統(tǒng)可視為最簡單的混合系統(tǒng)-線性瞬時混合系統(tǒng)。而構建獨立分量分析模型的目的就是僅利用觀測數(shù)據(jù)估計出或恢復出未知源信號的大部分或全部特征。

圖3 多源信號的“雞尾酒會”模型Fig.3 Cocktail party effect of multisource signal

為了獲得信號源s的一個良好估計y=Wx，引入目標函數(shù)-損失函數(shù)ρ(y，W)，并設其數(shù)學期望為：

R(W)=E{ρ(y，W)}，

(7)

其中R(W)為風險函數(shù)。由于相互獨立性的測度可以有多個衡量標注，在信號中最常用的如KL散度。其定義為：

(8)

在ICA問題中，q(y)是y的邊緣概率密度的乘積，可以表示為：

(9)

其中pi(yi)即為所求的邊緣概率密度，即：

(10)

(11)

其中ρ(y，W)=log(py(y)/q(y))。因此，可以獲得一個便于求解的獨立性測度R(W)。

為使代價函數(shù)(如R(W))最小化，可以利用對隨機梯度搜尋最優(yōu)解的思想獲得獨立分量分析的自然梯度模型：

(12)

(13)

其中激活函數(shù)向量f(y)=[f1(y1)，f2(y2)，f3(y3)，…，fn(yn)]T，為列向量，其單個分量為：

(14)

其中qi(yi)為源信號{si}的近似概率密度函數(shù)。

(15)

為了配合算法的實時性，結合式(13)后，式(15)可改寫為實時自然梯度算法如式(16)所示：

(16)

利用式(16)便可以從混合信號中通過迭代逐步分離出不同類型的信號，進而還原出整個信號的特點。

4 湍流頻譜估計與中心風速分離

在利用獨立分量分析對湍流進行頻譜分解前，首先要對獲得的信號進行空域解處理。在這一步預處理中，分子的熱運動和粒子的布朗運動被去相關化[13]。由于在最終的頻譜分析中，這兩者并不是本文關注的對象，因此這一處理會大幅降低二者對后續(xù)頻譜分離與風速估計的影響，從而提高分離與估計的準確度。

(17)

由輸入信號變?yōu)榘谆幚硇盘柕倪^程則可由式(18)表示：

(18)

其中Q矩陣稱為白化矩陣。當信號在信道傳輸過程中出現(xiàn)混疊時，白化信號矩陣Q由n×n變?yōu)閚×m(m

圖4 白化處理矩陣結構Fig.4 Framework of whiten matrix

因此，經(jīng)由白化矩陣Q處理之后的信號矩陣如下：

(19)

顯然，對這一白化矩陣Q左乘任何一個正交矩陣，式(19)結果不變，因此矩陣Q對信號處理而言并非唯一解。

由于零均值信號的相關矩陣Rxx=E{xxT}通常為對稱正定矩陣，因此，根據(jù)矩陣理論可分解為：

(20)

其中Vx為正交矩陣，Λx=diag{λ1，λ2，λ3，…，λn}為對角矩陣，對應特征值由大到小排列，且由于其為正定矩陣，故所有特征值均大于等于0。

因此，白化矩陣Q的一種求法可以表示為：

(21)

也可以表示為：

(22)

假設兩個信號s1和s2是相互獨立的，如圖5(a)所示，將兩信號合成后，其結果如圖5(b)所示。此時兩信號x1和x2含有相關成分，若直接帶入后續(xù)處理中，則會產(chǎn)生混疊影響估計質量，而通過預白化處理后可以達到如圖5(c)所示的效果，圖中y1與y2分量與s1和s2分量在能量分布上均可以視為獨立的。與s1和s2的區(qū)別僅為分布位置的變化，而與信號自身的獨立性無關。

綜上，白化預處理可以使具有相關性的混合信號解相關，并削弱非相關信號對混合信號的影響，在此基礎上進行的頻譜分離則具有更高的穩(wěn)定性和準確度。

(a)獨立信號成分(a)Independent signals(b)獨立成分混合(b)Mixed signals(c)混合成分白化(c)Whiten signals圖5 混合獨立分量與白化處理示意圖Fig.5 Schematic diagram of mixed and whiten signals

5 測量實驗與結果

基于之前的理論分析，對ICA在風速處理中的實際表現(xiàn)進行了測試，在26 ℃恒溫地下實驗室搭建了激光測風雷達光路并結合電路部分，如圖6所示。圖7所示為激光測風雷達系統(tǒng)樣機，該樣機主要包括基于伽利略望遠鏡結構設計的雙目激光雷達天線，基于Mach-Zehnder干涉系統(tǒng)設計的延時自外差光纖測量光路，經(jīng)光電轉換后轉變?yōu)镸Hz級電信號，再經(jīng)數(shù)字電路簡單處理后可得原始信號(Raw Signal)，將原始信號傳遞至上位機，即可利用ICA方法對原始信號進行湍流分離與風場估計。

圖6 測風雷達光路部分Fig.6 Optical system of wind lidar

圖7 測風雷達系統(tǒng)樣機Fig.7 Prototype of lidar

圖8 混合風速測量結果Fig.8 Raw signal in wind measurement

首先利用轉盤進行速度校準，通過控制轉盤速度及測量位置在0～3 m/s，每隔0.5 m/s采樣10組對回光信號進行校準，其結果表2所示。

表2 轉盤轉速與回光頻率關系Tab.2 Relationship between disk speed andfrequency of backscatter light

在校準中發(fā)現(xiàn)，1 MHz以下部分波動劇烈，頻率校準方差較大且與實際速度關聯(lián)性弱，因此考慮它為光路內反饋或電磁干擾等形成的噪聲。因此利用高通濾波器對噪聲進行濾除，隨后利用擋板障礙物與風扇組合構建復雜風場環(huán)境。每獲得一組1 000×1 024個有效風速數(shù)據(jù)點(約10 min)，之后對測風雷達進行一次定速轉盤速度校準，在校準結果誤差小于5%時繼續(xù)測量，否則停機進行進一步校驗并放棄前一組數(shù)據(jù)?；旌巷L速測量結果如圖8所示。在獲得有效數(shù)據(jù)后對復雜風場回光信號進行預白化處理，得到的白化信號如圖9所示。在白化過程中，雙目信號中具有相關性的部分被逐步分離，通過解相關運算，強相關部分(可認為是湍流或風速的“中心”)被保留，而弱相關部分(可認為是湍流或風速的“邊緣”)則轉化為本底噪聲。這樣的信號在ICA處理中速度更快且收斂性更好。

圖9 原始信號與白化及濾波處理后的信號Fig.9 Raw and whiten signals

利用梯度下降的ICA方法可以將雙目信號分解為相互獨立的信號，在頻譜中表現(xiàn)為兩個信號峰，如圖9所示，其中ICA Processed曲線即為ICA處理結果，其頻譜中的獨立峰參數(shù)如表3中所示。

表3 雙目信號ICA方法在風速分離處理中頻譜獨立峰的平均參數(shù)

Tab.3 Main parameters of two peaks measured by ICA processing

如表3所示，雖然在實際的風速信號處理中，ICA方法的表現(xiàn)不如模擬信號，但是同樣可以獲得較為穩(wěn)定且信噪比滿足要求的頻譜峰值。結合風速的實際情況與包含湍流的風場特性，以高斯函數(shù)和正弦函數(shù)作為先驗信息中的參考進行獨立分量的分離，得到ICA Processed 曲線對應的分離結果，如圖10所示。將一組中所有數(shù)據(jù)分為50個小組進行50次ICA分離，其結果穩(wěn)定性如表4所示。由此可以看出，在可用數(shù)據(jù)量減少時，ICA方法處理中雖然會出現(xiàn)少量錯誤結果，但正確辨識率依舊能夠保證在90%以上。

圖10 原始信號、白化信號與ICA方法處理后的信號Fig.10 Raw， whiten and ICA processed signals

表4 50次ICA分離結果統(tǒng)計
Tab.4 ICA separation results for 50 tests

分離結果次 數(shù)正確識別率/%完全分離46中心風速與渦流混合192.0未解混2無法辨識的混合結果1

圖11 ICA方法所獲得風速頻譜峰值的1 s穩(wěn)定性曲線Fig.11 Two peaks values of wind frequency spectra measured by ICA processing in 1 s

對一組中所有數(shù)據(jù)進行統(tǒng)計可以得到頻域穩(wěn)定性情況，如圖11所示。低頻區(qū)域曲線穩(wěn)定性較差，而上方曲線穩(wěn)定性較好，因此可以初步判斷上方曲線為中心風速，而下方曲線為湍流中心的可能性較大。若考慮到低頻區(qū)域會受到更大的噪聲干擾時，可以結合主成份分析或自適應濾波等方法作為輔助信號處理方法進一步判斷中心風速的頻率。

6 結 論

本文應用ICA方法，對雙目激光測風雷達回光信號的頻譜進行分解，可以從中分離出兩個獨立信號，得到信噪比分別為25.93和31.01 dB的湍流信號和中心風速信號估計，通過觀測其頻域穩(wěn)定性，認為其中(2.59±0.05)MHz為中心風速，能夠實現(xiàn)在湍流大氣中捕捉穩(wěn)定中心風速信號的功能。