基于激光雷達回波的動態尾渦特征參數計算

趙麗雅，谷潤平，魏志強

(中國民航大學空中交通管理學院，天津，300300)

飛行安全間隔很大程度上受限于飛機翼尖產生的尾渦，對尾渦特征參數的分析有助于提高機場跑道利用率、減少飛機延誤，而通過激光雷達探測特定區域內的尾渦并解算其特征參數是一種重要的研究手段。

Rahm等[1]利用機載多普勒雷達測量大型飛機在自由大氣中的尾渦，通過對探測數據的處理得到尾渦運動軌跡。徐世龍等[2]從飛機尾渦激光探測的總體要求出發，設計了尾渦相干激光探測系統的結構與參數。牛鳳梁等[3]對降雨條件下飛機尾流的雷達特性進行研究，得到尾流中雨滴數密度分布特性，基于某雷達參數通過信號處理得到降雨條件下飛機尾流的雷達回波，并發現了其獨特的多普勒特性。Wu等[4]在北京首都國際機場利用激光雷達探測飛機尾渦，對雷達探測系統及探測方式進行了詳細討論。Penkin等[5]在俄羅斯機場利用多普勒激光雷達探測到B737-800飛機的尾渦并進行了動力學分析。Hon等[6]介紹了香港國際機場的飛機尾渦觀測計劃及初步數據分析結果。上述研究主要利用雷達探測尾渦軌跡，并定性分析了尾渦的位置與強度變化，沒有對尾渦的特征參數進行定量反演計算。

利用雷達回波對尾渦特征參數進行定量研究的成果也有不少。K?pp等[7]在Tarbes機場對大型飛機產生的尾渦進行探測，并利用速度包絡法設計了距離-高度顯示器(range-height indicator，RHI)掃描方式下的尾渦特征參數反演算法。吳永華等[8]研究了飛機尾渦激光探測技術，給出激光探測方式和探測系統的參量。徐世龍等[9]利用速度包絡法設計了平面-位置顯示器(plan-position indicator，PPI)掃描方式下的尾渦特征參數反演算法，并以某型號飛機的激光雷達探測實驗數據為例進行了計算。屈龍海[10]圍繞飛機尾流雷達探測需求，研究了晴空、云霧以及降雨條件下飛機尾流的雷達散射特性。Smalikho等[11-12]利用“Stream-Line”脈沖相干多普勒激光雷達對飛機尾渦進行探測，提出了根據雷達數據估計尾渦軌跡的方法，并采用徑向速度法估算飛機尾渦的特征參數。Loaec等[13]用激光雷達探測不同機型的尾渦強度，分析了背景風場的影響。Darracq等[14]鑒于機載雷達的位置特殊性，設計了利用軸向速度推演尾渦特征參數的算法。

上述研究中提出的算法多是針對靜態尾渦的特征參數反演，沒有考慮在實際大氣中尾渦的強度消散和渦核下沉，影響到特征參數的計算精度。本文從飛機尾渦的物理模型出發，分析尾渦附近的徑向速度場與特征參數的關系，并結合尾渦的運動及消散模型，建立動態尾渦特征參數反演算法。最后通過仿真實驗來驗證該算法的有效性。

1 飛機尾渦的物理模型

1.1 尾渦環量模型

飛機尾渦的強度一般用渦旋環量Γ來表示。在理想無黏氣體中，流場內的環量處處相等且保持不變。但通常情況下，在渦旋中心處受到空氣黏性影響較大，其環量為零，而且渦旋環量隨著時間的變化也會逐漸減小。同時，由于相互之間的誘導作用和重力因素，渦旋形成之后會一邊旋轉擴散一邊下降。

在尾渦剛剛形成還沒有擴散衰減時，初始環量Γ0取決于飛機的質量、速度、翼展等參量[15]：

(1)

1.2 尾渦速度模型

尾渦旋轉速度直接關系著尾渦的形狀、大小和演化進程。尾渦速度模型很多，本文選用Hallock-Burnham模型描述尾渦切向速度分布。該模型與尾流的擬合度好、形式簡單、速度變化相對比較平滑。與渦核中心距離為r的某一點的切向速度VT為：

(2)

式中:rc為渦核半徑。渦核半徑為各個方向上切向速度最大值的點與渦核的距離，通常為(3%～5%)b0。

1.3 尾渦消散模型

尾渦消散分為兩個階段：擴散階段和快速消散階段。在擴散階段，尾渦消散的速度非常慢，其強度近似認為不變；在快速消散階段，尾渦強度迅速減小。設尾渦擴散階段的持續時間為tc，同時為了計算方便，通常還定義相對于基準參數的無量綱參數——參考時間t0與參考下沉速度ω0：

(3)

在尾渦擴散階段其強度消散約為10%，且與時間倒數近似為負指數關系，擴散階段的尾渦消散模型為[16]：

(4)

經過時間tc后，尾渦消散會加快，快速消散階段的尾渦消散模型為[16]：

(5)

式中：N為浮力頻率。

1.4 尾渦下沉模型

尾渦產生以后，不僅隨著時間逐漸消散，并且由于相互的誘導作用而向下運動。尾渦下沉速度ω與渦核間距及渦核半徑有關[17]：

(6)

尾渦經過t秒后下降的高度h按下式計算：

ωdt

(7)

式中：h0為初始尾渦高度。

2 靜態尾渦特征參數計算方法

2.1 徑向速度

如圖1所示，激光雷達對飛行中的飛機左右機翼產生的一對尾渦進行扇形掃描。分別以激光雷達中心和兩個渦核的中點為坐標原點建立坐標系YOZ和Y′O′Z′，RO為兩個坐標原點之間的距離。設渦核中心分別為O1和O2，其與激光雷達的徑向距離為R1和R2。選取尾渦流場內任意一點P作為研究對象，P與激光雷達的距離為R，激光雷達掃描到P點時的探測角度為φP，αP為尾渦徑向速度與Y軸的夾角。

圖1 徑向速度計算示意圖

Fig.1Schematicdiagramofthecalculationofradialvelocity

激光雷達測量尾渦流場的速度實質上是測量大氣分子或氣溶膠粒子的運動速度，所以只有與激光波束傳播方向一致的速度分量才能引起后向散射光信號的多普勒頻移。本文使用RHI掃描方式，掃描結果只能獲取尾渦剖面上各點速度沿激光波束的徑向分量，即切向速度在掃描徑向上的投影值，稱為徑向速度VR。通過圖1所示的坐標關系，利用式(2)計算出尾渦附近某一點P受到單個渦影響時的徑向速度為：

αP-φP)

(8)

同時由式(8)可推導出尾渦流場內任意一點受到兩個渦影響時的徑向速度為[7]：

(9)

其中，

Fl(R,φ)=

rl=ROsin(φ-φO)-(-1)l(b0/2)sinφ，

Y(R,φ)=Rcosφ-ROcosφO，

Z(R,φ)=Rsinφ-ROsinφO。

式中：l為尾渦編號，l=1代表右渦，l=2代表左渦。

通過分析可知rc與其他數值相比較小，可以忽略不計，為方便計算同時又不失精度，將上式中的Fl(R,φ)簡化為：

(10)

2.2 渦核位置

讓激光雷達在預定義區域內掃描，將探測到的回波數據進行降噪、插值等處理后，得到有較高速度分辨率的多普勒譜，利用各個探測單元內的多普勒譜計算出尾渦徑向速度。在每個探測單元的多普勒譜中選取適當的門限值，如圖2所示，該門限值與多普勒譜的交點即為尾渦流場中的速度正包絡與負包絡，分別記為Vmax與Vmin，代表各探測單元內徑向速度的最大值和最小值。

圖2 利用多普勒譜得到的速度包絡值

兩個尾渦的位置可利用上述速度包絡來確定。分別計算不同探測單元的正、負速度包絡，如圖3所示。選取正包絡的最大值點A、B，其坐標

圖3 正、負速度包絡值

(11)

2.3 尾渦環量

速度包絡值代表的是每個探測單元內徑向速度的正向最大值和負向最小值，所對應的點為探測單元內的徑向速度極值點，其導數為零，即

(12)

通過求解式(12)可得出：

(13)

計算出正、負包絡上的速度理論值，將式(13)帶入式(9)，得出：

(14)

式中：d(φ)=b0cosφ。式(14)所示即為每個單元內徑向速度正向最大點或負向最小點的速度。

假設在尾渦流場中的某些位置只受到單個渦的影響，則可利用這些點的速度計算單個尾渦的環量大小。如圖4所示，選取與激光雷達距離為R1、R2處的點i和k，探測到i、k點的激光雷達仰角分別為φi、φk。

圖4 尾渦環量計算示意圖

Fig.4Schematicdiagramofvortexcirculationcalculation

由式(14)推導出：

(15)

式中：Γ1ik與Γ2ik為利用i點和k點估算得到的兩個尾渦的環量值。

選取與激光雷達的距離為R1的n個數據點以及與激光雷達的距離為R2的m個數據點進行計算，通過式(15)得到環量估計值Γ1ik和Γ2ik。將這些估算值求和后取平均值，得出根據徑向速度計算尾渦環量的公式如下：

(16)

2.4 尾渦特征參數計算流程

綜上所述，基于激光雷達回波的靜態尾渦特征參數計算流程為：①激光雷達發射激光波束對尾渦周圍進行掃描，并接受后向散射信號；②將接收到的信號進行降噪、傅里葉變換等信號處理得到多普勒譜；③選取門限值確定正、負速度包絡，進而得到探測區域的徑向速度分布；④利用徑向速度分布，通過式(11)反演得到渦核位置；⑤將渦核坐標值帶入式(16)計算尾渦環量。

2.5 實例分析

2.5.1 靜態尾渦流場

下面利用A320飛機的尾渦流場進行計算分析，飛機的相關參數如表1所示。

表1 A320飛機相關參數

假定尾渦形成后強度及高度不變，仿真出靜態尾渦的激光雷達回波。激光雷達與機場跑道中軸線水平距離為700 m，探測開始時飛機距離地面約60 m，激光雷達掃描參數如表2所示。

表2 激光雷達掃描參數

經過對回波信號的處理和徑向速度的解算，得到如圖5所示的掃描區域徑向速度分布。

圖5 掃描區域徑向速度分布

根據徑向速度反演出不消散、不下沉的靜態尾渦左右兩渦的環量分別為283.9731、282.9731 m2/s，與表1中尾渦初始環量290m2/s比較接近。

2.5.2 動態尾渦流場

上述算法在計算靜態尾渦環量時的精度較高，然而大氣中尾渦產生后會不斷下沉且逐漸消散，因此以下分析該算法在計算動態尾渦特征參數時的性能。

利用靜態尾渦特征參數計算方法，結合尾渦下沉運動模型反演出兩個尾渦的渦核高度，如圖6中虛線所示，圖中實線為利用式(7)計算得到的渦核高度隨時間的變化情況，顯然二者存在較大偏差，估算出的渦核位置較高。

圖6 靜態尾渦特征參數算法計算出的動態尾渦高度

Fig.6Heightofdynamicwakevortexcalculatedbythestaticvortexcharacteristicparameteralgorithm

尾渦下沉和消散都可能導致上述算法對環量估算精度的降低。本文根據尾渦下沉與消散的變化情況，仿真動態尾渦的回波信號，并利用該回波信號計算尾渦的特征參數，與尾渦流場內利用機型參數的計算結果進行對比。由于渦核附近流場速度會發生突變，可導致計算值存在一定偏差，因此利用探測到左渦渦核位置之前10 s的數據點以及探測到右渦渦核位置之后10 s的數據點進行計算，不同時刻下尾渦環量的計算結果如圖7所示，其中實線為根據文獻[15]中的研究結果得出的A320飛機尾渦環量消散情況。由圖7可見，估算結果與尾渦流場的環量變化擬合度較差。

圖7 靜態尾渦特征參數算法計算出的動態尾渦環量

Fig.7Circulationofdynamicwakevortexcalculatedbythestaticvortexcharacteristicparameteralgorithm

將根據不同時刻的激光雷達探測數據計算得到的環量值加和取平均得到左、右尾渦環量的最終計算結果分別為365.63、363.67 m2/s，這與利用機型參數計算得到的環量值(290 m2/s)的誤差在25%左右。

綜上所述，對于靜態尾渦，即在不考慮尾渦的消散和運動的情況下，式(11)～式(16)所代表的算法能夠較為準確地計算出尾渦的特征參數。但實際中尾渦的位置與強度不斷變化，通過該方法估算出的尾渦特征參數與尾渦流場內的實際情況相差較多，因此需要考慮尾渦的消散與下沉特性對算法進行優化，以提高特征參數計算精度。

3 尾渦特征參數計算方法的改進

3.1 渦核位置的計算

由于尾渦的下沉運動，其渦核位置會發生變化，故以下對式(11)中的靜態尾渦渦核位置確定方法進行改進。

首先，通過式(11)找到渦核位置，即為當前探測時刻的渦核位置，解算出該時刻渦核與激光雷達之間的水平距離x及垂直距離y：

x=Rlcosφl,y=Rlsinφl

(17)

然后，根據尾渦的下降速度計算尾渦不同探測時刻t的渦核高度yt：

yt=y+ωt

(18)

最后，根據垂直距離重新計算每一個探測時刻t的尾渦距離Rlt和探測角度φlt：

(19)

3.2 尾渦環量的計算

根據式(19)計算出的渦核位置更接近大氣中運動尾渦的實際位置，利用該結果反演出的尾渦環量也會更加精確。將式(19)替換式(16)中的Rl和φl(l=1,2)，得到只考慮尾渦下沉而不考慮消散的動態尾渦環量計算公式為：

(20)

然而，尾渦在下沉運動的同時，其環量也要逐漸衰減。公式(15)中，用于計算速度的環量應為不同時刻的環量值，即

(21)

式中：Γ10和Γ20為兩個尾渦的初始環量;p(t)為尾渦消散函數;ti和tk分別為探測到i點和k點的時刻。

由式(21)推導出同時考慮尾渦下沉和尾渦消散的動態尾渦環量計算公式：

(22)

3.3 動態尾渦特征參數計算仿真

下面利用仿真的動態消散、下沉的尾渦流場作為參照組，與通過激光雷達回波計算得出的特征參數進行對比，分析改進算法的準確性。

首先，只考慮尾渦下沉的影響，利用式(20)估算出尾渦的環量大小，結果如圖8中菱形和圓形線所示，圖中實線為根據文獻[15]得出的結果，虛線為根據尾渦消散模型(式(4)～式(5))得出的結果。由圖8可見，估算出的尾渦環量偏小，這是因為公式(20)沒有考慮尾渦的消散情況。但是與圖7對比可以看出，即使只改進尾渦的渦核位置確定方法，最后的尾渦環量估算精度也有明顯提高。

圖8 只考慮尾渦下沉影響的優化算法計算結果

Fig.8Calculationresultsofoptimizationalgorithmonlyconsideringtheeffectofvortexdescent

然后，同時考慮尾渦消散和下沉的影響，通過式(22)估算出尾渦的初始環量。利用不同時刻的估算值，結合式(4)～式(5)的尾渦消散模型計算出不同時刻的尾渦環量值，其結果如圖9所示。從圖9來看，利用式(22)估算出的尾渦初始環量更接近尾渦消散模型和文獻[15]的研究結果。

將上述實驗結果中不同時刻估算的初始環量值求和取平均得到最終反演出的尾渦初始環量值。只考慮尾渦下沉影響因素時計算出的左、右尾渦初始環量分別為267.64、260.21 m2/s，對比通過機型參數計算出的尾渦初始環量290 m2/s，其誤差分別約為7%與10%；而綜合考慮尾渦消散和下沉因素時計算出的左、右尾渦初始環量分別為281.25、281.40 m2/s，其誤差約為3%。前面2.5.2節的研究結果表明，不考慮尾渦下沉和消散影響的環量計算誤差在25%左右。由此可以看出，與尾渦消散相比，尾渦下沉對尾渦環量的估算精度影響更大。這是因為，渦核位置影響到計算使用的數據點與渦核的距離r的大小，從而涉及每一個數據的計算，累計誤差較大，導致環量計算偏差較大，而所選取用于計算的數據點處于渦核附近，這時尾渦還未進入或者剛剛進入快速消散階段，環量衰減較小，故對最終的環量計算精度的影響程度相對較低。

圖9同時考慮尾渦下沉和消散影響的優化算法計算結果

Fig.9Calculationresultsofoptimizationalgorithmconsideringtheeffectsofvortexdescentanddecay

4 結語

在實際情況中，飛機尾渦產生后是逐漸消散并且不斷下沉的，而靜態尾渦特征參數計算方法對尾渦環量的估算偏差較大。本文綜合考慮了尾渦的消散和下沉因素，提出了利用激光雷達回波計算動態變化的飛機尾渦特征參數的方法。首先使用特定門限值提取尾渦流場的速度包絡，解算出尾渦徑向速度分布，然后利用該徑向速度計算出渦核位置與尾渦環量。該方法能夠更加準確地反演出動態尾渦環量。利用仿真的動態尾渦流場激光雷達探測回波進行計算分析顯示，相較于尾流消散，渦核下沉運動對尾渦環量的計算精度影響更大，所以對渦核位置的精確計算能夠有效提高尾渦整體特征參數的估算精度。