航空自组网是指利用空中节点(如民航客机)构建的移动自组织网络，为节点间或节点与地基站间提供航空通信/导航/监管/空管(CNS/ATM)等服务。与卫星通信相比，航空自组网具有时延低、部署容易、数据传输率高等特点^[1]，成为学术界近年的研究热点^{[2, 3, 4, 5]}。现有研究成果主要集中于参照航空通信环境设计更优的空空路由或接入协议，较少关注空空链路的物理层传输体制与信道特性^{[2, 3]}。空空链路作为航空自组网的唯一数据传输通道，在FAA和EUROCONTROL联合提出的L波段数字航空通信系统(L-DACS)中具有非常重要的研究价值^[6]。为确定空空链路物理层的各关键参数，如符号周期、子载波间隔、分集距离等，研究其散射环境的二阶统计特性，如时域自相关、多普勒扩展等，以期为下一代航空移动通信的高效空空组网提供必要的物理层信息。

信道二阶统计特性与诸多通信指标有关，尽管研究较为复杂, 却一直受到业界的广泛关注^{[7, 8]}。文献[9]首次推导了二维场景下基站固定、终端移动、散射体静止时接收信号包络的时域自相关函数和功率谱, 为后续通信系统的设计与性能分析提供了理论和技术支撑。文献[10]将该研究结果推广到三维场景。然而航空自组网中作为收发端的空中节点均在三维立体空域内作高速运动(如民航客机速度可达1 000 km/h)^[11]，这类无线通信环境为典型的移动-移动(mobile-to-mobile, M2M)散射环境。除地表上静止的各类散射体以外(如高山)，大气层中还广泛分布着相对于空中节点作高速的随机运的各类散射体(如云、雨、雾、尘埃、电离层反射等)，由此造成的多径效应将表现出明显的非稳态特性。为分析航空自组网空空链路的二阶统计特性，必须将为地面蜂窝网建立的Clarke模型推广到三维，并融入运动节点与非稳态散射体的影响。以车联网为应用背景，文献[12, 13]推导了二维非稳态M2M散射场景的时域自相关函数，但未将其结果推广到三维。文献[14, 15, 16, 17, 18]研究了收发节点相对静止或匀速直线运动、散射体分布于给定三维区域内(如圆柱体、球体、椭球体等)的M2M散射场景的时域自相关函数，但没有考虑散射体相对于收发节点作随机运动的情况。上述研究为航空自组网非稳态三维M2M散射场景的时域自相关特性的分析提供了可参考的解决思路。针对航空移动通信环境，文献[19, 20]通过实测分析了空空链路几个典型场景的接收信号强度的统计特性(即信道的一阶统计特性)。文献[21]采用有限马尔科夫模型研究了三维M2M散射场景下接收信号幅度的统计特性。遗憾的是上述研究均没有研究空空链路的二阶统计特性，也没有考虑散射体相对于空中节点作随机运动对接收信号统计特性的影响，不能全面评估航空自组网非稳态三维M2M散射信道的统计特性。

结合航空移动通信环境的三维非稳态特征，本文推导了立体空域内散射体相对于收发节点随机运动的M2M散射场景的时域自相关函数的数学表达式。该表达式为高维积分，本文采用蒙特卡洛方法求得其数值解。用射线跟踪法建立仿真验证模型进一步证明该表达式的正确性，并用该模型对不同节点运动速度、不同散射体速度和不同散射体数量情况下的时域自相关函数进行了分析。

1 理论分析 1.1 非稳态三维散射模型

图 1为非稳态三维散射场景示意图。设收发端周围有${N_S}$个运动散射体${S^n}$($n = 1,2, \cdots ,{N_S}$)，图 1详细展示了第n条散射路径的主要参数：

1) 任意时刻t，坐标原点为收端A_R所在空间位置，三维坐标的_x轴为收发节点的直视路径，x轴的正方向由发端A_T指向收端A_R。

2) 发端相对于收端的速度大小为${v_{\rm{T}}}$，其运动方向的水平角和仰角分别为${\alpha _{\rm{T}}}$和${\beta _{\rm{T}}}$；由于本文模型考虑的是相对速度，故可以模拟收发节点运动的M2M场景。

3) 发端A_T到第n个散射体Sⁿ的发射角的水平角和仰角分别为$\alpha _{\rm{T}}^n$、$\beta _{\rm{T}}^n$；第n个散射体Sⁿ到收端A_R的到达角的水平角和仰角分别为$\alpha _{\rm{R}}^n$、$\beta _{\rm{R}}^n$。

4) 第n个散射体Sⁿ以固定速度、沿随机方向相对于收发节点运动，其运动速度、运动方向的水平角和仰角分别为${v_{{S^n}}}$、$\alpha _v^{{S^n}}$、$\beta _v^{{S^n}}$，以此模拟散射体随机运动造成的非稳态散射场景。

5) $\alpha _{\rm{R}}^{{\rm{LOS}}} = {\rm{\pi }}$表示直视路径水平角。

1.2 非稳态三维散射场景的时域自相关函数

第1.1节所描述的三维非稳态散射模型将文献[12]提出的二维模型扩展到三维。研究表明由多径造成的复信道增益可表示为^[22]：

$ \mu (t) = \sum\limits_{n = 1}^{{N_S}} {{c_n}{{\rm{e}}^{{\rm{j}}(2{\rm{\pi }}{f_n}t + {\theta _n})}}} $

(1)

式中，c_n表示衰减因子；f_n表示由收发节点、散射体移动产生的多普勒频移；θ_n表示路径的相移。将二维场景推广到三维场景的难点在于把仰角β引入到f_n的表达式中。沿用文献[12]的思路，把f_n分成$f_n^{\rm{T}}$、$f_n^{\rm{TS}}$、$f_n^{\rm{SR}}$共3个部分，即：

$ {f_n} = f_n^{\rm{T}} + f_n^{{\rm{TS}}} + f_n^{{\rm{SR}}} $

(2)

图 2为按照文献[10]的三维模型思路建立的散射体到收端的散射路径示意图。设散射体Sⁿ在三维坐标下的初始位置坐标为(0,0,0)，设电磁波传播方向为$(\alpha _{\rm{R}}^n,\beta _{\rm{R}}^n)$，散射体运动方向为$(\alpha _v^{{S^n}},\beta _v^{{S^n}})$，经过单位时间后散射体移动至$({x_0},{y_0},{z_0})$。多普勒频移是由收发端相对运动给电磁波传播带来光程差造成的，则第n条路径的多普勒频移可表示为^[10]：

$ \begin{array}{*{20}{c}} {f_n^{{\rm{SR}}} = }\\ {\frac{1}{\lambda }({x_0}\cos \alpha _{\rm{R}}^n\cos \beta _{\rm{R}}^n + {y_0}\sin \alpha _{\rm{R}}^n\cos \beta _{\rm{R}}^n + {z_0}\sin \beta _{\rm{R}}^n)} \end{array} $

(3)

又因为单位时间内接收端的运动距离可表示为$\left( {{v_{{S^n}}}\cos \alpha _v^{{S^n}}\cos \beta _v^{{S^n}},{v_{{S^n}}}\cos \beta _v^{{S^n}}\sin \alpha _v^{{S^n}},{v_{{S^n}}}\sin \beta _v^{{S^n}}} \right)$，将其代入式(3)并依次代替$({x_0},{y_0},{z_0})$，三角函数化简后，可得：

$ \begin{array}{*{20}{c}} {f_n^{{\rm{SR}}} = \frac{{{v_{{S^n}}}}}{\lambda }[\cos (\alpha _v^{{S^n}} - \alpha _{\rm{R}}^n)\cos \beta _{\rm{R}}^n\cos \beta _v^{{S^n}} + }\\ {\sin \beta _v^{{S^n}}\sin \beta _{\rm{R}}^{{S^n}}]} \end{array} $

(4)

式中，λ为载波波长。令$\beta _{\rm{R}}^n{\rm{ = 0}}$，β_v^Sn，则$f_n^{{\rm{SR}}} = \frac{{{v_{{S^n}}}}}{\lambda }[\cos (\alpha _v^{{S^n}} - \alpha _{\rm{R}}^n)]$，退化为二维情况。类比式(4)可分别得到$f_n^{\rm{T}}$和$f_n^{\rm{TS}}$表达式，故f_n表示为：

$ \begin{array}{l} {f_n} = \frac{{{v_{\rm{T}}}}}{\lambda }[\cos ({\alpha _{\rm{T}}} - \alpha _{\rm{T}}^n)\cos \beta _{\rm{T}}^n\cos {\beta _{\rm{T}}} + \sin \beta _{\rm{T}}^n\sin {\beta _{\rm{T}}}]{\rm{ + }}\\ \frac{{{v_{{S^n}}}}}{\lambda }[\cos (\alpha _v^{{S^n}} - \alpha _{\rm{T}}^n)\cos \beta _{\rm{T}}^n\cos \beta _v^{{S^n}} + \sin \beta _{\rm{T}}^n\sin \beta _v^{{S^n}}] + \end{array} $

(5)

文献[12]提出的自相关函数求解办法可得到非稳态三维散射场景下的时域自相关函数表达式为：

$ a(\tau ) = 2\sigma _0^2\int {\int_{S \times S} {\int {\int_{S \times S} {\int {\int_{S \times S} {{{\rm{e}}^{{\rm{j}}2{\rm{\pi }}{f_n}\tau }}p(\alpha _v^{{S^n}},\beta _v^{{S^n}})p(\alpha _{\rm{T}}^n,\beta _{\rm{T}}^n,\alpha _{\rm{R}}^n,\beta _{\rm{R}}^n)p({v_s})} } } } } } {{\rm{d}}_{\alpha _v^{{S_n}}}}{{\rm{d}}_{\beta _v^{{S_n}}}}{{\rm{d}}_{\alpha _{\rm{T}}^n}}{{\rm{d}}_{\beta _{\rm{T}}^n}}{{\rm{d}}_{\alpha _{\rm{R}}^n}}{{\rm{d}}_{\beta _{\rm{R}}^n}} $

(6)

把式(5)代入式(6)，可得到非稳态三维信道场景下的时域自相关函数为：

$ \begin{array}{c} a(\tau ) = 2\sigma _0^2\int_{{\rm{ 0}}}^{{\rm{ }}2{\rm{\pi }}} {\int_{{\rm{ 0}}}^{{\rm{ \pi }}} {\int_{{\rm{ 0}}}^{{\rm{ }}2{\rm{\pi }}} {\int_{{\rm{ 0}}}^{{\rm{ \pi }}} {\int_{{\rm{ 0}}}^{{\rm{ }}2{\rm{\pi }}} {\int_{{\rm{ 0}}}^{{\rm{ \pi }}} {{{\rm{e}}^{{\rm{j}}\left\{ {2{\rm{\pi }} \cdot \frac{{{v_{{S^n}}}}}{\lambda }\left[ {\cos (\alpha _v^{{S^n}} - \alpha _{\rm{T}}^n)\cos \beta _{\rm{T}}^n\cos \beta _v^{{S^n}} + \sin \beta _{\rm{T}}^n\sin \beta _v^{{S^n}} + \cos (\alpha _v^{{S^n}} - \alpha _{\rm{R}}^n)\cos \beta _{\rm{R}}^n\cos \beta _v^{{S^n}} + \sin \beta _{\rm{R}}^n\sin \beta _v^{{S^n}}} \right]\tau } \right\}}}} } } } } } \times \\ {{\rm{e}}^{{\rm{j}}\left\{ {2{\rm{\pi }}\frac{{{v_{\rm{T}}}}}{\lambda }\left[ {\cos ({\alpha _{\rm{T}}} - \alpha _{\rm{T}}^n)\cos \beta _{\rm{T}}^n\cos {\beta _{\rm{T}}} + \sin \beta _{\rm{T}}^n\sin {\beta _{\rm{T}}}} \right]\tau } \right\}}}{\kern 1pt} \frac{{\sin \beta _{\rm{T}}^n}}{{4{\rm{\pi }}}}\frac{{\sin \beta _{\rm{R}}^n}}{{4{\rm{\pi }}}}\frac{{\sin \beta _v^{{S^n}}}}{{4{\rm{\pi }}}}{{\rm{d}}_{\alpha _v^{{S^n}}}}{{\rm{d}}_{\beta _v^{{S^n}}}}{{\rm{d}}_{{\alpha _{\rm{T}}}^n}}{{\rm{d}}_{\beta _{\rm{T}}^n}}{{\rm{d}}_{\alpha _{\rm{R}}^n}}{{\rm{d}}_{\beta _{\rm{R}}^n}} \end{array} $

(7)

式(6)到式(7)的推导过程为：

1) 设发射角和到达角相互独立，故联合概率密度函数$p(\alpha _n^{\rm{T}},\beta _n^{\rm{T}},\alpha _n^{\rm{R}},\beta _n^{\rm{R}})$可表示为$p(\alpha _n^{\rm{T}},\beta _n^{\rm{T}})$与$p(\alpha _n^{\rm{R}},\beta _n^{\rm{R}})$的乘积。

2) 设发端在三维球面内等概率发出电磁波信号，则球面上单位面积能接收到电磁波辐射的概率为${\rm{d}}p = {\rm{d}}\Omega /4{\rm{\pi }}$，其中${\rm{d}}\Omega = \sin \beta _{\rm{T}}^n{\rm{d}}\beta _{\rm{T}}^n{\rm{d}}\alpha _{\rm{T}}^n$，故$\alpha _{\rm{T}}^n$和$\beta _{\rm{T}}^n$的联合概率密度函数$p(\alpha _{\rm{T}}^n,\beta _{\rm{T}}^n) = \sin \beta _{\rm{T}}^n/4{\rm{\pi }}$^[23]。

3) 设电磁波信号在三维球面内等概率到达收端，故到达角的水平角$\alpha _{\rm{R}}^n$和仰角$\beta _{\rm{R}}^n$的联合概率密度函数$p(\alpha _{\rm{R}}^n,\beta _{\rm{R}}^n) = \sin \beta _{\rm{R}}^n/4{\rm{\pi }}$。

4) 设散射体以恒定速度在三维空间内沿任意方向作随机运动，参考步骤2)，散射体速度的水平角$\alpha _v^{{S^n}}$和仰角$\beta _v^{{S^n}}$的联合概率密度函数满足$p(\alpha _v^{{S_n}},\beta _v^{{S_n}}) = \sin \beta _v^{{S_n}}/4{\rm{\pi }}$；散射体运动速度${v_s}$的概率密度函数$p({v_s}) = 1$。

如果多径数目趋于无穷，故复信道增益等价为一个高斯过程，根据中心极限定理得到该过程的方差$2\sigma _0^2 = 1$，但实际情况中多径数目不可能趋于无穷，因此式(7)中保留了该常数项。式(7)为一个相当复杂的六重积分，难以用符号积分给出其精确值，但可采用数值积分求解。传统的辛普森或梯形法在解决高维积分时，因求解时间会随积分重数增加而迅速增加，严重影响计算效率。蒙特卡洛数值积分方法是一种以概率统计理论为指导的一类数值计算方法，具有较高的运算精度且计算效率不受积分重数影响，因此本文采用蒙特卡洛方法计算该六重积分的数值解。

2 基于射线跟踪法的仿真验证模型

为进一步验证式(7)的正确性，本文建立图 1所示散射方式的射线跟踪仿真验证模型。由文献[24]可得：多普勒频移${f_n} = 1/(2{\rm{\pi }}\Delta t)\Delta \varphi $，其中$\Delta \varphi = 2{\rm{\pi }}\Delta l/\lambda $，Δl表示Δt时间内的路程差，将其代入式(1)，可简化为：

$ \mu (t) = \sum\limits_{n = 1}^{{N_S}} {{c_n}{{\rm{e}}^{{\rm{j}}\left( {\frac{{2{\rm{\pi }}}}{\lambda }{\rm{\Delta }}l} \right)}}} $

(8)

结合式(8)与射线跟踪法，本文可给出非稳态三维散射场景的仿真验证办法，具体实现流程如下：

1) 设置初始参数：包括收端坐标，发端坐标，载波频率f_c，散射体个数N_s，散射体运动速度v_Sⁿ，发端相对于收端的运动速度大小v_T、相对速度水平角${\alpha _{\rm{T}}}$和仰角${\beta _{\rm{T}}}$。

2) 计算散射体Δt时刻的三维坐标：设散射体最初随机分布在收端(原点)周围，不局限于任何球形或柱形空间，随机产生第n个散射体运动方向的水平角α_v^Sn和仰角β_v^Sn，最小时间间隔t_d内散射体的运动路程${\rm{\Delta }}S = {V_{{S^n}}}{t_d}$，把该路程分解到(x,y,z)3个方向上，其中，${\rm{\Delta }}{S_x} = {\rm{\Delta }}S\cos \beta _v^{{S^n}}\cos \alpha _v^{{S^n}}$，${\rm{\Delta }}{S_y} = {\rm{\Delta }}S\cos \beta _v^{{S^n}} \cdot $ $\sin \alpha _v^{{S^n}}$，${\rm{\Delta }}{S_z} = {\rm{\Delta }}S\sin \beta _v^{{S^n}}$，再把Δt(Δt=Nt_d，即N个时间间隔)时间内(x,y,z)方向上的路程分别进行累加，加上散射体初始位置坐标，可得每个散射体Δt时刻后的三维坐标。

3) 计算Δt时刻的路程差Δl：计算第n个散射体Sⁿ与发端A_T、收端A_R之间的欧式距离$\left| {{S^n}{A_{\rm{T}}}} \right|$和$\left| {{S^n}{A_{\rm{R}}}} \right|$，两段距离相加得到Δt时刻A_T与A_R的散射路程长度$\left| {{A_{\rm{T}}}{A_{\rm{R}}}} \right|$；若Δt=0时，A_T与A_R的散射路程长度为$\left| {{A_{\rm{T0}}}{A_{\rm{R0}}}} \right|$，那么$\Delta l = \left| {\left| {{A_{\rm{T}}}{A_{\rm{R}}}} \right| - \left| {{A_{{\rm{T0}}}}{A_{{\rm{R0}}}}} \right|} \right|$，具体可参见图 1。

4) 计算时域自相关函数：把Δl带入式(8)得到多径信号Δt时刻的合成信号，结合自相关函数的定义$a(\tau ) = E\left\{ {\mu *(t)\mu (t + \tau )} \right\}$即可求得非稳态三维散射场景的时域自相关函数。

采用射线跟踪模型可避开高微积分求解难的问题，从几何学角度直观地模拟非稳态三维散射场景。该模型可模拟三维空域内空中节点的相对运动和散射体以不同密度随机分布的情况；令$\beta _{\rm{T}}^n = 0$，$\beta _{\rm{R}}^n = 0$，$\beta _v^{{S^n}} = 0$，该模型还可模拟非稳态散射二维场景。

3 数值仿真结果

数值仿真验证部分的初始参数为：收端坐标(0,0,0)，发端坐标(−2 000,0,0)(单位：m)，载波频率f_c=900 MHz，散射体个数${N_{\rm{S}}} = 1{\rm{ }}000$，发端相对运动方向的水平角${\alpha _{\rm{T}}} = {\rm{\pi }}/4$，仰角${\beta _{\rm{T}}} = {\rm{\pi }}/3$。除特别说明，均按以上参数设置。

令$\beta _{\rm{R}}^n = 0$、$\beta _v^{{S^n}} = 0$、$${v_{\rm{T}}} = 0{\rm{ }}{\kern 1pt} {\kern 1pt} {\rm{m/s}}$$，式(7)退化为文献[12]中的二维非稳态散射场景下的时域自相关函数。图 3是文献[12]和本文式(7)的数值结果，分别采用了蒙特卡洛高维积分法和第3节的射线跟踪法。3条曲线具有较好的重合度，尤其是当$\tau $较大时，验证了式(7)的正确性，也说明本文所建三维模型可将文献[12]的二维场景作为一个特例。

令${v_{{S^n}}} = 100{\rm{ m/s}}{\rm{ }}{v_{\rm{T}}} = 100{\rm{ m/s}}$，图 4为三维场景下射线跟踪法对式(7)的验证结果。两者具有较强的一致性，进一步验证了三维场景下式(7)的正确性。

令${v_{\rm{T}}} = 10{\rm{ m/s}}$，图 5为散射体速度从${v_{{S^n}}} = 0{\rm{ m/s}}$到${v_{{S^n}}} = 20{\rm{ m/s}}$的时域自相关函数。令${v_{\rm{T}}} = 500{\rm{ m/s}}$，图 6为散射体速度从${v_{{S^n}}} = 0{\rm{ m/s}}$到${v_{{S^n}}} = 400{\rm{ m/s}}$变化时的时域自相关函数。前者模拟空中节点相向飞行(相对速度较小)的情况，后者则对应空中节点相背运动(相对速度较大)的场景。当收发节点之间的相对速度较小时，尽管散射体速度比较小，但仍可影响接收端信号的时域自相关函数；当收发节点之间的相对速度较大时，时域自相关函数变化明显加快，散射体速度大到一定程度才能对接收端时域自相关函数产生影响。

航空移动通信中散射体数目会随着天气和飞行高度变换而变化。大气层中云层稀薄或节点高空飞行时，散射体数目较少；而当大气层中云层比较厚或节点近地飞行时，散射体数目将明显增加。令${v_{\rm{T}}} = 500$m/s，${v_{{S^n}}} = 300$m/s，图 7散射体个数N_S从10～到10 000时的时域自相关函数。散射体密度越大，时域自相关函数衰减得越快，表明运动散射体确实是产生多普勒频移的重要源头。

上述曲线给出了不同节点运动速度、不同散射体运动速度和不同散射体数量情况下，非稳态三维散射场景时间自相关函数的变化趋势。该信息可估计信道建模的重要参数。对时间自相关函数$\left| {a(\tau )} \right|$进行傅里叶变换可得到信道的多普勒功率谱$S({f_d})$，多普勒功率谱衰减到一定程度的频带宽度为信道的多普勒扩展${B_d}$，由此可得到信道的相干时间${T_c} = 1/{B_d}$^[25]。研究成果可为下一代航空移动通信的高效空空组网提供必要的物理层信息。如下一代航空移动通信的OFDM、MIMO系统的符号周期应小于相干时间以减小时间选择性衰落，时域导频间隔应小于相干时间以保证信道估计质量，天线之间的相干距离$d = v{T_c}$(v为电磁波传播速度)，发射分集时间应大于相干时间以达到多径分集抗衰落的作用。

4 总 结

本文将二维非稳态散射场景的自相关函数求解办法推广到了三维场景，给出该数学表达式，建立射线跟踪法的仿真模型，通过与蒙特卡洛高微积分数值解进行比较验证了该表达式的正确性。数值仿真表明三维非稳态散射场景的时域自相关函数的衰减速度随节点运动速度、散射体运动速度和分布密度的增加而增加。该模型与求解办法具有较强的推广性：1) 通过降维可等效为二维散射场景；2) 可在三维空域内模拟空中节点相对运动；3) 可模拟空中节点周围散射体不同密度随机分布的情况。研究成果可为下一代航空移动通信的高效空空组网提供必要的物理层信息。若融入地面通信环境的基本特征，还可用于车联网等典型地面移动自组网的通信链路的统计特性分析。本文后期将进一步研究融入MIMO体制的非稳态三维散射场景下的信道二阶统计特性。