基于移動擴展模型的無人機三維 空地非平穩信道建模

摘 要： "無人機空地信道幾何統計模型通常采用固定的速度和移動方向，無法描述UAV空地傳播環境的非平穩特性。面向基于雙圓柱散射體的UAV三維空地信道模型，提出在空地建模中引入具有動態速度和移動方向的高斯馬爾可夫模型，模擬無人機在現實環境中的變化；為了使該模型模擬不同的運動，引入擴展因子切換模型的運動軌跡；此外，運用幾何分析法，根據時變的速度和運動方向推導了空時相關函數和多普勒功率譜密度，并研究了無人機水平和垂直運動對信道統計特性的影響。仿真結果表明，擴展的模型僅會導致時域的非平穩，對空域沒有影響，并驗證了擴展模型的通用性。

關鍵詞： "無人機； 空對地； 高斯馬爾可夫模型； 非平穩

中圖分類號： "TN929.5 """文獻標志碼： A

文章編號： "1001-3695（2022）02-036-0537-06

doi：10.19734/j.issn.1001-3695.2021.07.0319

3D air-to-ground non-stationary channel modeling of UAV "based on extended mobility model

Liu Min， Zhang Zhizhong， Deng Bingguang

（School of Communications amp; Information Engineering， Chongqing University of Posts amp; Telecommunications， Chongqing 400065， China）

Abstract： "Geometry-based stochastic modeling（GBSM） for UAV air-to-ground（A2G） channel usually assumes fixed velocity and moving direction，which can’t describe the non-stationary characteristics of the UAV air-ground propagation environment.Facing the UAV 3D air-to-ground channel model based on double cylindrical scatterer，this paper proposed to introduce a Gauss Markov model with dynamic speed and moving direction into the A2G to simulate the changes of the UAV in the real environment.In order to make the model simulate different motions，it introduced an expansion factor to switch the motion trajectory of the model.In addition，this paper used the geometric analysis method to derive the space-time correlation function and Doppler power spectral density based on the time-varying speed and direction of motion，and studied the influence of the UAV horizontal and vertical movement on the channel statistical characteristics.The simulation results show that the extended model only cause non-stationarity in the time domain without affecting the air space，and verifies the generality of the extended model.

Key words： "unmanned aerial vehicle（UAV）； A2G； Gauss Markov model； non-stationary

0 引言

近年來，無人機已經在軍事和民用領域發揮著重要的作用［1］，如天氣監測、交通控制、貨物運輸、通信中繼等。UAV通道具有一些獨特的通道特性，如3D部署、高機動性、時空非平穩性以及機身陰影［2］。與傳統的“車—車”“移動—移動”模型不同，UAV通常在三維平面內移動，并且以相對較低的高度飛行，這時應考慮到建筑物和路邊障礙等散射體條件。因此，傳統的信道模型無法直接應用于UAV通信場景［3，4］。

目前，UAV的信道建模和測量工作已經大量展開，從研究結果來看， UAV信道模型可以分為確定性模型［5］、非幾何統計性模型（non-geometric statistical model，NGSM）［6］和基于幾何的統計性模型（geometry-based stochastic model，GBSM）。確定性建模通常具有較高的準確性，但需要大量數據來表征特定的傳播環境；而非幾何統計性建模是基于測量的，具有一定的復雜度。與上述建模方法不同，基于幾何的統計性建模具有較高的準確性和較低的復雜度，因此該方法能廣泛用于模擬各種無線信道。

空對地信道的接收端和發送端運動引起的非平穩性在文獻［7］的測量活動中可以觀察到。然而，在過去有關無人機信道建模的工作中，非平穩性往往被忽略。文獻［8］提出了三維窄帶單球模型，并研究了相關的統計特性；文獻［9］提出了適用于UAV多輸入多輸出（multiple-input and multiple-output，MIMO）瑞利信道的三維單圓柱模型，并研究了UAV的運動方向和位置對空時相關函數的影響；文獻［10］提出一種用于無人機通信的3D橢圓柱MIMO信道模型，并且空時相關函數被無人機的移動性所影響。然而上述研究均假設信道是廣義平穩的，即信道統計特性不隨時間變化。根據UAV信道測量［11］表明，廣義平穩假設只在非常短的時間間隔內有效。在現實的無人機通信場景中，UAV和地面端（ground side，GS）的運動導致參數隨時間變化，這將會導致UAV傳播時的非平穩性。因此，研究UAV的非平穩性至關重要。

當前UAV通信的局限性是未考慮收發端的移動模型。在無人機網絡的仿真中，移動模型能準確描述移動節點在真實環境中的運動軌跡［12，13］。文獻［14］提出了一系列的移動模型，然而在大多數UAV信道模型中，收發端的速度和移動方向被簡單設為固定值，這與實際的UAV場景不符，并會導致不準確的信道統計特性。文獻［15］利用時變的布爾變量描述視距路徑、地面反射和散射路徑的動態生滅過程來反映信道的非平穩特性；文獻［16］用無人機旋轉引起的時變傳播距離來反映空對地信道的非平穩性，忽略了速度和運動方向對其的影響；文獻［17］對于寬帶截斷橢圓信道，地面端僅考慮了固定不變的速度和移動方向； 文獻［18］提出一種適用于UAV-MIMO信道的3D半球模型，其考慮了由于UAV和散射體的快速移動而導致的非平穩傳播。然而，這些基于GBSM的研究都沒有考慮UAV的移動模型對信道統計特性的影響。在UAV信道仿真中，需要精確的移動模型來描述UAV終端的移動行為。先前的工作缺少用速度以及運動方向的變化來表征信道的非平穩性，過度簡化了移動軌跡，導致無法準確描述無人機的移動軌跡。

基于此，本文參考文獻［19］中的雙圓柱模型，考慮了視距、單反射、雙反射和地面反射四種不同的多徑分量，并研究了無人機水平運動和垂直運動對時域和空域的影響，得出了新的結論。由于上述文獻假設信道是廣義平穩的，沒有考慮信道的非平穩性。本文首先引入高斯馬爾可夫移動模型來體現信道的非平穩性，并根據動態的速度和運動方向推導了空時相關函數和多普勒功率譜密度；并且在城市無人機空地場景下對高斯馬爾可夫模型進行擴展，擴展的模型具有廣泛性，能描述不同的運動軌跡；根據擴展的模型推導了無人機四種不同的運動。最后通過蒙特卡羅［20］仿真方法比較了不同的運動軌跡對無人機時空相關函數和多普勒功率譜密度的影響，該結論為非平穩信道的研究提供了有效的幫助，并驗證了擴展移動模型的通用性。

1 信道模型

1.1 系統模型

本文采用雙圓柱UAV-MIMO通信系統。該模型假設大量散射體分別位于發送端（ T X ）和接收端（ R X ）的兩個半徑為 R T 和 R R 的圓柱表面上。本文沒有像文獻［19］一樣假設固定的速度和運動方向，而是考慮了具有動態速度和運動方向的高斯—馬爾可夫移動模型。雙圓柱UAV-MIMO模型如圖1所示，該模型的發送端配備有 M T根天線，接收端配備有M R根天線，其中M T=M R=2，兩個圓柱體分別位于無人機和地面端。假設N 1個散射體分布在半徑為R T的圓柱體表面上，用S（n 1）表示第n 1個（n 1=1，2，…，N 1）有效散射體，類似地，假設N 2個散射體分布在半徑為R R的圓柱體表面上，并且用S（n 2）表示第n 2個（n 2=1，2，…，N 2）有效散射體。同時考慮了由地面引起的地面反射，地面反射發生在地面端半徑為R的圓柱底部，并且0≤R≤R R。地面端的散射體S（n 3）表示第n 3個（n 3=1，2，…，N 3）有效散射，半徑R的分布用概率密度函數f R=2R/R2 R 表示。表1定義了圖中所涉及的參數。

根據圖1中基于幾何的3D模型，針對空對地的通信環境考慮了四個非視距分量：a）單反射射線SBU，經 T X 端的圓柱表面的散射體 S（n "1） 反射到達 R X ；b）單反射射線SBG，波在到達 R X 之前，經 R X 周圍的圓柱表面散射體 S（n "2） 反射；c）地面反射射線SBGR，經 T X 端的圓柱底面散射體 S（n "3） 反射再到達 R X ；d）雙反射射線DB，波從 T X 的圓柱表面的散射體 S（n "1） 散射到 R X 周圍的圓柱表面，再經過 R X 的圓柱表面 S（n "2） 散射到達 R X 。第 p 根UAV天線和第 q 根地面端天線之間的信道沖激響應是視距、SBU、SBG、SBGR、DB分量之和，可以寫為

h pq（t）=hLOS pq（t）+hSBU pq（t）+hSBG pq（t）+hSBGR pq（t）+hDB pq（t） ""（1）

hLOS pq（t）= "KΩ pq K+1 ""e -j ""2π "λ ε pq "e j2 π t v T λ （ sin "βLOS T "sin "ξ T） ""e j2 π t v T λ [（ cos （αLOS T-γ "T） "cos "βLOS T）× cos "ξ T] e j2 π t v R λ ［（ cos （αLOS R-γ ""R） "cos "βLOS T） ""（2）

hSBU pq（t）= "η SBUΩ pq K+1 ""lim "N "1→+∞ "1 "N 1 "∑ N 1 n "1=1 "e jφ（n "1） e -j 2 π "λ （ε pn "1+ε n "1q）

e j2 π t v T λ （ sin "β（n "1） T "sin "ξ T） e j2 π t v T λ [（ cos （α（n "1） T-γ ""T） "cos "β（n "1） T）× cos ξ T]

e j2 π t v R λ ［（ cos （α（n "1） R-γ ""R） cos "β（n "1） R） ""（3）

hSBG pq（t）= "η SBGΩ pq K+1 ""lim "N 2→+∞ "1 "N 2 "∑ N 2 n 2=1 "e jφ（n 2） e -j 2 π "λ （ε pn 2+ε n 2q）

e j2 π t v T λ （ sin "β（n "2） T "sin "ξ T） e j2 π t v T λ [（ cos （α（n "2） T-γ ""T） "cos "β（n "2） T）× cos "ξ T]

e j2 π t v R λ ［（ cos （α（n "2） R-γ ""R） cos "β（n "2） R） """（4）

hSBGR pq（t）= "η SBGRΩ pq K+1 ""lim "N 3→+∞ "1 "N 3 "∑ N "3 n 3=1 "e jφ（n "3） e -j 2 π "λ （ε pn "3+ε n "3q）

e j2 π t v T λ （ sin "β（n "3） T "sin "ξ T）

e j2 π t v T λ [（ cos （α（n "3） T-γ ""T） "cos "β（n 3） T）× cos "ξ T]

e j2 π t v R λ ［（ cos （α（n "3） R-γ ""R） "cos "β（n "3） R） ""（5）

hDB pq（t）= "η DBΩ pq K+1 ""lim "N 1，N 2→+∞ "1 "N 1N 2 nbsp;∑ N 1，N 2 n 1，n 2=1 "e jφ（n "1，n "2）

e -j 2 π "λ （ε pn "1+ε n "1n "2+ε n "2q）

e j2 π t v T λ （ sin "β（n "1） T "sin "ξ T） e j2 π t v T λ [（ cos （α（n "1） T-γ ""T） "cos "β（n 1） T）× cos "ξ T]

e j2 π t v R λ ［（ cos （α（n "2） R-γ ""R） "cos "β（n "2） R） ""（6）

式（2）～（6）中的 ε pq、ε pn "1、ε pn "2、ε pn "3、ε n "1q、ε n "2q、ε n "3q、ε n "1n "2 分別表示鏈路 p-q、p-s（n "1）、p-s（n "2）、p-s（n "3）、s（n "1）-q、s（n "2）-q 、 s（n "3）-q 和 s（n "1）-s（n "2） 的距離； λ 是載波波長； Ω pq 是 p-q 鏈路的總能量； K 是萊斯因子。

此外， η SBU、η SBG、η SBGR、η DB 分別表示各個能量分量占總能量 Ω pq/K+1 的比例，并且滿足 η SBU+η SBG+η SBGR+η DB=1 。相位 φ（n "1）、φ（n "2）、φ（n "3）、φ（n "1，n "2） 是均勻分布在[-π，π）上的獨立隨機變量。根據空間距離公式、散射體和天線之間的距離， ε pq、ε pn "1、ε pn "2、ε pn "3、ε n "1q、ε n "2q和ε n "3q 表示成

ε pq= （DTR x）2+（DTR y）2+（DTR z）2 """（7）

ε pn "1= （DTn "1 x）2+（DTn "1 y）2+（DTn "1 z）2 """（8）

ε n "1q= （Dn "1R x）2+（Dn "1R y）2+（Dn "1R z）2 """（9）

ε pn "2= （DTn "2 x）2+（DTn "2 y）2+（DTn "2 z）2 """（10）

ε n "2q= （Dn "2R x）2+（Dn "2R y）2+（Dn "2R z）2 """（11）

ε pn "3= （DTn "3 x）2+（DTn "3 y）2+（DTn "3 z）2 """"（12）

ε n "3q= （Dn "3R x）2+（Dn "3R y）2+（Dn "3R z）2 """"（13）

式（7）中 DTR x=Δ T "cos "θ T-D-Δ R "cos "θ R ， DTR y=Δ T "sin "θ T- sin "θ R，DTR z=D "tan "β 0 ；式（8）中 DTn "1 x=Δ T "cos "θ T-R T "cos "α（n "1） T，DTn "1 y=Δ T "sin "θ T-R T "sin "α（n "1） T，DTn "1 z=D "tan "β 0+H 0-R T "tan "β（n "1） T；式（9）中Dn "1R x=R T "cos α（n "1） T-Δ R "cos "θ R-D，Dn "1R y=R T "sin "α（n "1） T-Δ R "sin "θ R，Dn "1R z=R T "tan "β（n "1） T-H 0；式（10）中DTn "2 x=-D+Δ T "cos "θ T-R R "cos "α（n "2） R，DTn "2 y=Δ T "sin "θ T-R R "sin "α（n "2） R，DTn "2 z=D "tan "β 0-R R "tan "β（n "2） R+H 0 ；式（11）中 Dn "2R x=R R "cos "α（n "2） R-Δ R "cos "θ R ， Dn "2R y=R R "sin "α（n "2） R-Δ R "sin "θ R，Dn "2R z=R R "tan "β（n "2） R-H 0；式（12）中DTn "3 x=-D+Δ T "cos "θ T- cos "α（n "3） RR，DTn "3 y=Δ T· "sin "θ T

，

DTn "3 z=D· "tan "β 0+H 0；式（13）中Dn "3R x=R "cos "α（n "3） R-Δ R· "cos "θ R，Dn "3R y=R "sin "α（n "3） R-Δ R "sin "θ R，Dn "3R z=H 0 。

參數 Δ T 是第 p 根天線元件和UAV天線陣列的中心之間的距離，參數 Δ R 是第 q 根GS端的天線元件和GS端天線陣列的中心之間的距離，對于均勻線性陣列，它們定義為

Δ T= 1 2 （M T+1-2p）δ T ""（14）

Δ R= 1 2 （M R+1-2q）δ R ""（15）

對于單反射射線，離開角與到達角具有一定的幾何關系，并且可以相互轉換。對于SBU射線，離開角與到達角有如下關系：

cos "α（n "1） R≈-1 ""（16）

sin "α（n "1） R≈ "R T D "sin α（n "1） T 1- R T D "sin α（n "1） T """（17）

cos "β（n "1） R≈ cos "β 0+ R T D "sin "β 0 "cos "β 0A（n "1） T ""（18）

sin "β（n "1） T≈ sin "β 0- R T D ""cos 2 β 0A（n "1） T ""（19）

其中： A（n "1） T≈ tan "β（n "1） T "cos "β 0+ cos "α（n "1） T "sin "β 0 。

對于SBG射線，離開角與到達角有如下關系：

cos "α（n "2） T≈1 ""（20）

cos "α（n "2） T≈ "R R D ""sin "α（n "2） R 1+ R R D ""cos "α（n "2） R """（21）

cos "β（n "2） T≈ cos "β 0+ R R D "sin "β 0 "cos "β 0A（n "2） R ""（22）

sin "β（n "2） T≈ sin "β 0- R R D ""cos 2 β 0A（n "2） R ""（23）

其中： A（n "2） R≈ tan "β（n "2） R "cos "β 0+ cos "α（n "2） R "sin "β 0 。

對于SBU射線，離開角與到達角有如下關系：

cos "α（n "3） T≈-1 ""（24）

cos "α（n "3） T≈ "R D "sin α（n "3） R 1+ R D "cos α（n "3） R """（25）

cos "β（n "3） T≈ cos "β 0- R D "sin "β 0 "cos "β 0A（n "3） R ""（26）

sin "β（n "3） T≈ sin "β 0+ R D "cos 2 β 0A（n "3） R ""（27）

其中： A（n "3） R≈ tan "β（n "3） R "cos "β 0- cos "α（n "3） R "sin "β 0 。

實際上，提出模型是理論的參考模型。假設散射體的數量是無限的，即 N 1，N 2，N 3→∞ 。在這種情況下，連續的表達式 α R 、 β R 、 α T 、 β T 能用離散的表達式 α（n i） R 、 β（n i） R 、 α（n i） T 、 β（n i） T 代替。假設方位角和俯仰角是獨立的，聯合的概率密度函數為 f（α，β）≈f（α）f（β） 。

本文用von Mises PDF分布來描述方位角 α T 和 α R 。

f（α）= "e k "cos "（α-α μ） 2 π I 0（k） "- π≤ α ≤π "（28）

其中： I 0（·） 是第一類零階修正Bessel函數； α Tμ 、 α Rμ 表示散射體在水平面上分布的中心角，并且 α μ ∈［-π，π）。參數 k 控制中心角 β μ 的信號傳播，其隨著環境的非各向同性散射而增加。俯仰角 β T 、 β R 用余弦分布表示為

f（β）= "π "4β m ""cos（ π "2 ""β-β μ β m ） |β-β μ|≤β m ≤ π 2 """（29）

其中： β∈［β μ-β m，β μ+β m］；β m 是中心仰角的最大偏移量； β μ 是中心角。

最后，假設地面端的散射體在 O′ 和 R R 地面上均勻分布，用 f（R）=2R/R2 R 來描述半徑為 R 的分布。

1.2 移動模型

移動模型的隨機性能反映復雜多變的無人機通信，并且模擬無人機不規則的運動，從而體現信道的非平穩性特性。本文將文獻［12］中的二維高斯馬爾可夫移動性模型擴展到三維平面，模擬無人機的運動軌跡，并用二維模型模擬地面端的運動。為了提高模型的通用性，本文擴展了高斯—馬爾可夫移動模型，并用擴展的模型描述UAV和GS的運動軌跡。通過引入擴展因子 χ s 來切換模型，進一步擴展模型的移動軌跡。該模型能表示四種不同的運動，當 χ s=1 時，模型能模擬布朗運動和不規則運動；當 χ s=0 時，模型能描述勻速運動和勻變速運動。在該模型中，假設發送端的速度和方向是時間上相關的高斯馬爾可夫過程。應用式（30），根據 t i-1 時刻速度和方向的值及隨機變量，可以計算出 t i 時刻的運動軌跡。

v T（t i）=χ v Tv T（t i-1）+（1-χ v T）μ v T+χ s（X i-1 1-（χ v T）2 ） ξ T（t i）=χ ξ Tξ T（t i-1）+（1-χ ξ T）μ ξ T+χ s（Y i-1 1-（χ ξ T）2 ） γ T（t i）=χ γ "Tγ T（t i-1）+（1-χ γ "T）μ γ "T+χ s（Z i-1 1-（χ γ "T）2 ） ""（30）

其中： i=1，2，3，… ；時間間隔 t=t i-t i-1；χ v T/ξ T/γ "T 是調整隨機性的參數，并且 0≤χ v T/ξ T/γ "T≤1，當χ v T/ξ T/γ "T =0時，該模型類似于布朗運動，當 χ v T/ξ T/γ "T =1時，該模型變成線性模型。因此可以通過改變 χ v T/ξ T/γ "T 的值改變模型的隨機性。 v T（t i）、ξ T（t i）和γ T（t i） 分別為時刻 t i 的速度和運動方向；類似地， v T（t i-1）、ξ T（t i-1）和γ T（t i-1） 對應于時刻 t i-1 。當 i→∞ 時， μ v "T、μ ξ "T和μ γ "T 分別為速度和方向的平均值； X、Y、Z 是獨立不相關的高斯隨機變量，其均值為0，方差為1； χ s 是切換模型的參數，當 χ s=1 ，模型為高斯馬爾可夫移動模型，當 χ s=0 ，模型沒有了隨機性。

可以進一步用初始速度和方向表示式（30）為

v T（t i）=χi v "Tv T（t 0）+χ s（ 1-（χ v T）2 ""∑ i-1 j=0 ""χi-j-1 v "TX j）+（1-χi v T）μ v T ξ T（t i）=χi ξ "Tξ T（t 0）+χ s（ 1-（χ ξ T）2 ""∑ i-1 j=0 ""χi-j-1 ξ TX j）+（1-χi ξ "T）μ ξ T

γ T（t i）=χi γ "Tγ T（t 0）+χ s（ 1-（χ γ "T）2 ""∑ i-1 j=0 ""χi-j-1 γ "TX j）+（1-χi γ "T）μ γ "T """（31）

無人機的運動用該模型表示， v T（t 0） 是無人機 t 0 時刻的初始速度， ξ T（t 0） 是無人機 t 0 時刻運動方向的俯仰角， γ T（t 0） 是無人機 t 0 時刻移動方向的水平角，以上參數參見圖1。上述的移動模型適用范圍廣泛，可以表示四種不同的移動軌跡。

軌跡1 當 χ v T/ξ T/γ "T=0，χ s=1 時， v T（t i）=μ v T+X i-1，μ v T 為定值， X i-1 是獨立不相關均值為0、方差為1的高斯隨機過程。可以發現， v T（t i） 隨著 X i-1 的變化而變化， X i-1 的變化具有不確定性，因此 v T（t i） 的變化也具有不確定性。這種具有完全不確定性質的隨機運動叫做布朗運動。

軌跡2 當 χ v T/ξ T/γ "T=1 ， v T（t i）=v T（t i-1） ，移動模型的速度和運動方向不隨時間變化，這種運動稱為勻速運動。

軌跡3 當 0≤χ v T/ξ T/γ "T≤1 ， χ s=1 時， v T（t i）=χ v T·v T（t i-1）+（1-χ v T）μ v T+X i-1 1-（χ v T）2 "。該運動軌跡受 χ v T/ξ T/γ "T 和高斯隨機變量 X 的影響，并且 χ v T/ξ T/γ "T 能控制隨機性的程度。

軌跡4 當 χ s=0 時，模型喪失了隨機性，并且生成了更簡單的移動性軌跡，例如速度可以改為與時間有關的線性函數。用泰勒級數展開 χ s=0 時的速度表達式，具體如下：

v T（t i）=χi v T·v T（t 0）+（1-χi v T）μ v "T=

［v T（t 0）-μ v "T］·χi v T+μ vT=［v T（t 0）-μ v "T］· exp （i log e χ v T）+μ vT=

［v T（t 0）-μ vT］×［1+ i log e χ v T 1！ + （i log e χ v T）2 2！ +…］+μ v "T

當 χ v→1 ，只需要考慮泰勒級數的前兩項而舍棄其他項，仍可以確保其精度，從而得出

v T（t i）=［v T（t 0）-μ v "T］×［1+ i log e χ v T 1！ ］+μ v "T= i·（［v T（t 0）-μ v "T］ log e χ v T）+v T（t 0） ""（32）

因此得到一個勻變速模型， ［v T（t 0）-μ v "T］ log e χ v T 相當于加速度。可以發現當 χ v T=0.999 9 時，對于式（32）具有足夠的精度。同理可得到地面端 v R（t） 和 γ R（t） 的表達式。

1.3 空時相關特性

對于兩個任意的復數衰落包絡 h pq（t）和h p′q′（t） ，空—時相關函數定義為

R pq，p′q′（δ T，δ R，t ，Δ t）= E[h* pq（t）h p′q′（t +Δ t）] "Ω pqΩ p′q′ """"（33）

其中：（·）表示復共軛運算； E[·] 表示統計期望運算符。

由式（1）可以得到四個分量的疊加：

R pq，p′q′（δ T，δ R，t ，Δ t）=RLOS pq，p′q′（δ T，δ R，t， Δ t）+RSBU pq，p′q′（δ T，δ R，t ，Δ t）+

RSBG pq，p′q′（δ T，δ R，t ，Δ t）+RSBGR pq，p′q′（δ T，δ R，t ，Δ t）+RDB pq，p′q′（δ T，δ R，t ，Δ t） ""（34）

各部分表達式為

RLOS pq，p′q′（δ T，δ R，t ，Δ t）= K K+1 "e -j 2 π "λ （ε pq-ε p′q′）

e j2 π t（f LOS（t +Δ t）-f LOS（t））· e j2 πΔ tf LOS（t +Δ t） d α T d β T ""（35）

RSBU pq，p′q′（δ T，δ R，t ，Δ t）" = η SBU K+1 ∫β T "2 β T "1∫ π "-π f（α T）f（β T） e -j 2 π "λ [（ε pn "1+ε n "1q）-（ε p′n "1+ε n "1q′）]

e j2 π t（f SBU（t +Δ t）-f SBU（t））· e j2 πΔ tf SBU（t +Δ t） d α T d β T" " （36）

RSBG pq，p′q′（δ T，δ R，t ，Δ t） = η SBG K+1 ∫β R "2 β R1∫ π "-π f（α R）f（β R）

e -j 2 π "λ [（ε pn "2+ε n "2q）-（ε p′n "2+ε n "2q′）]

e j2 π t（f SBG（t +Δ t）-f SBG（t））·

e j2 πΔ tf SBG（t +Δ t） d α R d β R" " （37）

RSBGR pq，p′q′（δ T，δ R，t ，Δ t）=

η SBGR K+1 ∫R R 0∫ π "-π f（α R）f（R）e-j 2 π "λ [（ε pn "3+ε n "3q）-（ε p′n "3+ε n "3q′）] "e j2 π t（f SBGR（t +Δ t）-f SBGR（t））· e j2 πΔ tf SBGR（t +Δ t ）d α R d β R" " （38）

RDB pq，p′q′（δ T，δ R，t ，Δ t） =

η DB K+1 ∫β R "2 β R "1∫ π "-π ∫β T "2 β T "1∫ π "-πe j2 πΔ tf DB（t +Δ t）

e -j 2 π "λ [（ε pn2+ε n2q）-（ε p′n2+ε n2q′）] e j2 π t（f DB（t +Δ t）-f DB（t））

f（α R）f（β R） f（α T）f（β T） d α R d β R d α T d β R" " （39）

其中： f（t）= v T（t） λ ［（ cos（ α T-γ T（t）） "cos "β T） ×cos "ξ T（t）+ "sin "β T "sin "ξ T（t）］+ v R（t） λ ""cos "（α R-γ R（t）） "cos "β R， f LOS（t）、 f SBU（t）、 f SBG（t）、 f SBGR（t）、 f DB（t） 的表達式依此類推。

多普勒功率譜密度 S（f，t） 由時間相關函數的傅里葉變化得來，可以表示為

S（f，t） =∫∞ "-∞R pq，p′q′（0，0，t ，Δ t）· e- j 2π f Δ t dΔ t" " （40）

由于式（36）～（40）中的積分具有復雜的整體，很難用符號積分提供閉合的解。為了解決這個問題，本文使用蒙特卡羅進行數值積分，該方法在積分誤差和操作時間方面具有一定的優點，并且能較好地確定結果。因此利用此方法仿真了時間相關函數和多普勒功率譜密度。

2 數值結果與分析

本章將研究無人機的水平運動和垂直運動對空時相關函數的影響，所需設置的參數如下： D=100，H 0=5，λ=0.1，R T=5，R R=3，β 0=30°，θ T=90°，θ R=90°，k T=10，k R=3，α Tμ=45°，α Rμ=180°，β Tμ=45°，β Rμ=60°，β Tm=30°，β Rm=30°。無人機的初始速度設置為v T（t 0） =10 m/s，地面端的速度為 v R（t 0） =5 m/s，本文分別考慮無人機水平和垂直運動的情況。參數設置為 γ T（t 0）=γ R（t 0）=0°，μ v "T=μ v "R =10 m/s， μ γ "T=μ γ "R=0°。無人機水平運動方向的參數ξ T（t 0）=0°，μ ξ "T=0°。垂直運動的參數ξ T（t 0）=90°，μ ξ "T=90° 。所有軌跡的其他參數設置：高斯隨機變量 X、Y、Z 的標準方差分別設置為0.3、0.1、0.1， t step =0.05 s，觀察周期為0～2 s。

圖2為速度和運動方向與時間的相關圖。勻速運動的參數設置為 χ v T/ξ T/γ "T=1 。圖2（a）～（c）為勻速運動時速度和運動方向隨時間變化的曲線圖，可以發現速度和運動方向都不隨時間變化。勻變速運動的參數設置為 χ s=0 ， χ v T/ξ T/γ "T=0.999 9，v T（t 0） =0 m/s。由圖2（d）～（f）可以看出 v T 隨時間 增加，且非常接近線性函數，因此該模型可以近似看做加速度為 ［v T（t 0）-μ vT］ log e χ v T 的勻變速運動。隨機運動的參數設置為 χ s=1 ， χ v T/ξ T/γ "T=0.9 ；布朗運動的參數設置為 χ s=1 ， "χ v T/ξ T/γ "T=0 。對比圖2（g）～（i）和（j）～（l）可以觀察到，隨著 χ v T/ξ T/γ "T 減小，隨機度增加，速度和運動方向變化更加劇烈，這表明運動軌跡變得更加隨機且難以預測。

通過設置 δ T=δ R=0 消除空域分量來獲得時間相關函數。圖3展現了文獻［19］中廣義平穩信道不同的運動方向對時間相關性的影響。可以看出UAV水平運動時，時間相關函數隨時延的增加下降的速度比垂直運動快，這意味著UAV水平移動對時間相關函數的影響比垂直運動大。并且，DB分量的時間相關函數隨時延增大而減小的速度比其他分量更快，這說明DB分量對時延的變化更敏感。產生這個現象是因為DB分量比其他分量受更多散射體的影響。這些結果表明，UAV的運動方向對時間相關函數有著影響，UAV的水平運動可能導致時間相關函數減小。

圖4為Δ t =0時各個分量的水平和垂直運動的空間相關函數。可以清楚地觀察到，勻速運動時的運動方向對空間相關函數沒有影響，產生這個現象的原因是通過設置Δ t=0 ，消除了包含Δ t 的項。

勻速運動時的不同運動方向的多普勒功率譜密度如圖5所示，可以看出，垂直運動相對于水平運動的多普勒有所加寬，這意味著運動方向對多普勒密度有所影響。

由于無人機在現實環境中受風、氣壓等許多因素的影響，導致速度和運動方向呈現不規則的變化，所以隨機移動性模型更接近無人機的真實運動軌跡。圖6是在文獻［19］中引入動態的速度和運動方向的時間相關函數。為了清楚起見，設置 χ v T/v R=0.995，χ ξ "T/γ "T/γ "R=1 。從圖中可以看出 t=0.5 時，水平運動比垂直運動下降得略快。但是與垂直方向上UAV在 t=1.5 處的相關函數相比，垂直方向上UAV的時間相關函數隨時延的增加而減小得更快；垂直運動在 t=2 時的時間相關函數隨時延增加而下降的速度比垂直運動要快得多。可以看出，在不同時刻相關函數的形狀通常是不同的，這是由移動模型的隨機性引起的。對比圖3～6的SBU分量的時間相關函數可以發現，移動模型的隨機性對信道非平穩性有著顯著影響。在非平穩的信道中，由于隨機性的影響，UAV在水平和垂直運動中相關函數的下降速度之間的關系變得不一致。

圖7為無人機勻加速運動對時間相關函數的影響。UAV和地面端的初始速度設為0 m/s， χ v T/v R=0.999 9 。由圖7可以看出時間相關性隨時間的變化而變化，說明該運動能模擬信道在時域的非平穩性并且不會改變垂直和水平運動相對下降的速度。隨著時間的增加，也就是速度增加，時間相關性下降得越快。因此在空對地通信中，必須重視UAV快速的移動所引起的信道響應的變化。

圖8展示了UAV隨機運動對多普勒功率譜的影響，是在 t =0.5 s時SBU分量 χ v T 不同時的多普勒功率譜密度。可以觀察到，多普勒功率譜密度的形狀隨著 χ v T 的減小而發生變化的程度越明顯，證明移動模型的較強隨機性導致非平穩性增加。

通過設置Δ t=0 消除包含項，研究隨機運動的空間相關函數。結果發現，無論速度是否隨時間變化，隨機運動中的空間相關函數與勻速運動中的相同，這表明UAV的速度變化僅將非平穩性引入時域。

圖9將文獻［19］的理論值、文獻［21］的測量值與本文引入高斯—馬爾可夫模型后的時間相關函數進行對比。時間相關函數的測量數據來自文獻［21］中低海拔UAV通信，此處使用飛行距離代替相關時間，以符合文獻［21］中的描述。可以看出在雙圓柱模型中引入高斯馬爾可夫模型后，無人機的時間相關函數與測量值的擬合度比文獻［19］的擬合度更高，證明引入了移動模型的信道更符合實際的通信場景。

3 結束語

本文提出將高斯—馬爾可夫移動模型用于雙圓柱UAV-MIMO系統，來表示非平穩的空對地信道，并擴展該模型，讓其具有廣泛性。擴展后的模型不僅能調整模型的隨機性，還能描述隨機運動、線性運動等不同情況下的移動軌跡。基于雙圓柱模型，本文推導了時變的空時相關函數和多普勒功率譜密度，研究了UAV水平運動和垂直運動對統計特性的影響。仿真結果表明，動態運動僅會導致時域的非平穩，對空域沒有影響，并可能導致對時間相關函數的影響。通過與測量值的對比發現，引入移動模型更符合無人機實際的通信場景。這些研究可以作為具有逼真的移動軌跡的UAV-MIMO空對地非平穩信道建模的參考。在未來空天地一體化通信網絡中，無人機可以作為衛星對地通信的中繼站，在未來的研究中，還需要將信道模型推廣到中繼站對地的通信場景中。

參考文獻：

［1］ "Shakhatreh H，Sawalmeh A H，Al-Fuqaha A， et al. Unmanned aerial vehicles（UAVs）：a survey on civil applications and key research challenges［J］. IEEE Access ，2019， 7 ：48572-48634.

［2］ Wang Chengxiang，Huang Jie，Wang Haiming， et al .6G wireless channel measurements and models：trends and challenges［J］. IEEE Vehicular Technology Magazine ，2020， 15 （4）：22-32.

［3］ 賈汝冰.無人機場景下的寬帶信道建模與仿真［D］.北京：北京交通大學，2019. （Jia Rubing.Modeling and simulation of broadband channel in UAV scenario［D］.Beijing：Beijing Jiaotong University，2019.）

［4］ Jiang Hao，Zhang Zaichen，Wu Liang， et al .A non-stationary geometry-based scattering vehicle-to-vehicle MIMO channel model［J］. IEEE Communications Letters ，2018， 22 （7）：1510-1513.

［5］ Greenberg E，Levp P.Channel characteristics of UAV to ground links over multipath urban environments［C］//Proc of IEEE International Conference on Microwaves，Antennas，Communications and Electronic Systems.Piscataway，NJ：IEEE Press，2017.

［6］ Matolak D W，Sun Ruoyu.Air-ground channel characterization for unmanned aircraft systems，part I：methods，measurements，and models for over-water settings［J］. IEEE Trans on Vehicular Technology ，2017， 66 （1）：26-44.

［7］ Matolak D W，Sun Ruoyu.Unmanned aircraft systems：air-ground channel characterization for future applications［J］. IEEE Vehicular Technology Magazine ，2015， 10 （2）：79-85.

［8］ Jin Kun，Cheng Xiang，Xiao Huge， et al .Three dimensional modeling and space-time correlation for UAV channels［C］//Proc of the 85th IEEE Vehicular Technology Conference.Piscataway，NJ：IEEE Press，2017.

［9］ Zeng Linzhou，Cheng Xiang，Wang Chengxiang， et al. A 3D geometry-based stochastic channel model for UAV-MIMO channels［C］//Proc of IEEE Wireless Communications and Networking Conference.Pisca-taway，NJ：IEEE Press，2017.

［10］ Jiang Hao，Zhang Zaichen，Wu Jiang， et al .Three-dimensional geometry-based UAV-MIMO channel modeling for A2G communication environments［J］. IEEE Communications Letters ，2018， 22 （7）：1438-1441.

［11］ Cai Xuesong，Gonzalez-Plaza A，Alonso D， et al .Low altitude UAV propagation channel modelling［C］//Proc of the 11th European Conference on Antennas and Propagation.Piscataway，NJ：IEEE Press，2017：1443-1447.

［12］ Biomo J D M M，Kunz T，St-Hilaire M.An enhanced Gauss-Markov mobility model for simulations of unmanned aerial Ad hoc networks［C］//Proc of the 7th IFIP Wireless and Mobile Networking Confe-rence.Piscataway，NJ：IEEE Press，2014.

［13］ "He Ruisi，Ai Bo，Stüber G L， et al. Mobility model-based non-stationary mobile-to-mobile channel modeling［J］. IEEE Trans on Wireless Communications ，2018， 17 （7）：4388-4400.

［14］ Batabyal S，Bhaumik P.Mobility models，traces and impact of mobility on opportunistic routing algorithms：a survey［J］. IEEE Communications Surveys amp; Tutorials ，2015， 17 （3）：1679-1707.

［15］ 王亞文，朱秋明，陳小敏，等.基于幾何隨機的無人機空地信道模型［J］.數據采集與處理，2020， 35 （1）：118-127. （Wang Yawen，Zhu Qiuming，Chen Xiaomin， et al .Novel geometry-based stochastic model for UAV-based air-to-ground channels［J］. Journal of Data Acquisition amp; Processing ，2020， 35 （1）：118-127.）

［16］ Ma Zhangfeng，Ai Bo，He Ruisi， et al .Impact of UAV rotation on MIMO channel characterization for air-to-ground communication systems［J］. IEEE Trans on Vehicular Technology ，2020， 69 （11）：12418-12431.

［17］ Chang Hengtai，Bian Ji，Wang Chengxiang， et al .A 3D wideband geometry-based stochastic model for UAV air-to-ground channels［C］//Proc of IEEE Global Communications Conference.Piscataway，NJ：IEEE Press，2018：206-212.

［18］ Zhu Qiuming，Jiang Kaili，Chen Xiaomin， "et al .A novel 3D non-stationary UAV-MIMO channel model and its statistical properties［J］. China Communications ，2018， 15 （12）：147-158.

［19］ Cheng Xiang，Li Yiran，Wang Chengxiang， et al .A 3D geometry-based stochastic model for unmanned aerial vehicle MIMO Ricean fading channels［J］. IEEE Internet of Things Journal ，2020， 7（ 9）：8674-8687.

［20］ Mathews J，Fink D K.Numerical methods using MATLAB［M］.4th ed.New Jersey：Prentice-Hall Inc.，2005.

［21］ Simunek M，Fontán F P，Pechac P.The UAV low elevation propagation channel in urban areas：statistical analysis and time-series generator［J］. IEEE Trans on Antennas and Propagation ，2013， 61 （7）：3850-3858.