0 引 言
在航空气象领域，通常将飞机高度在600 m 以下，风向风速急剧变化的这一大气现象称为低空风切变，其具有变化快、强度大、危害性高等典型特征［1］。如果在起飞或降落阶段飞机遭遇风切变，而此时飞行员缺乏充足的时间来调整飞机状态，可能会导致严重飞行事故发生［2］。因此，有关风切变检测的问题一直以来都是民航领域关注热点。作为对于低空风切变检测流程的关键，风场速度估计的结果将直接影响整个风切变检测结果的准确程度，因此，准确估计风速至关重要［3-5］。
机载气象雷达是飞行中不可或缺的设备，可以对飞机航路上气象状况进行实时监测，及时获取气象数据发现危险天气，以帮助飞行员提前规避风切变等危险大气现象，从而保障飞行过程的安全［6］。随着航空领域事业的发展和进步，航空安全的实际使用需求对机载气象雷达的探测范围、扫描精度等性能提出了更高的要求。目前，与相对成熟的平面相控阵雷达相比，安装在飞机表面与飞行器贴合的共形阵机载气象雷达具有许多优点，如共形阵天线孔径的增大能提高目标探测时方位向的分辨率，天线扫描范围的扩大有助于实现全方位扫描等［7］。共形阵机载雷达所具备的优点可以更好地适应航空气象领域发展的需求。因此，对共形阵机载气象雷达体制下目标检测的课题意义重大。
机载气象雷达在检测风切变过程中风速估计是关键步骤，但由于地杂波分布范围较广，强度大，导致雷达接收到的风切变信号完全被地杂波信号所覆盖，进而导致风速估计结果的准确性受到影响。目前，针对低空风切变风速估计方法主要基于色加载知识辅助的空时自适应处理（Space-Time Adaptive Processing, STAP）风速估计方法［8］、基于直接数据域-扩展因子化法（Direct Data Domain-Extended Factored Approach, DDD-EFA）的风速估计方法［9］、基于同伦稀疏STAP 的风速估计方法［10］和基于改进辅助通道的风速估计方法［11］等。以上方法在传统线阵或面阵体制下都能准确估计风速，但对于阵元在空间上呈现两维或三维分布的半球共形阵，还未有文献对其风切变检测技术和风速估计方法进行研究报道。
基于此，本文提出了一种基于D3AR 风速估计方法。该方法将空时自回归算法（Space-Time Autoregression，STAR）思想引入到直接数据域（Direct Data Domain, DDD）算法中，首先将待检测单元的数据分别从空域、时域以及空时域进行信号对消处理，然后将处理后的数据矩阵描述为空时自回归模型并估计模型参数，再通过构造与杂波子空间正交的空间来实现对杂波的抑制，最后通过提取待检测单元的最大多普勒频率来估计风场速度。根据仿真结果显示，该方法有效地实现了地杂波抑制，并且能够精确估计风速。
1 半球共形阵雷达信号模型
1.1 地杂波信号
基于Ward 模型［12］仿真半球共形阵体制下的杂波回波信号。首先依据载机雷达参数确认雷达所照射区域内单个距离单元的宽度，然后将单个距离单元按照方位角划分为多个杂波块并计算每个杂波块的回波信号，通过将单个距离单元内的所有杂波块相加，得到杂波回波信号。
图1为半球共形阵机载气象雷达几何模型，假设半球共形阵共有M个半径分别为Rm的圆弧（顶点不放阵元），每一层圆弧上均匀分布着N个阵元，且相邻两个圆弧之间的间距为d，θ 为方位角，φ 为俯仰角，ψ为空间锥角，V为载机速度，由几何关系可知波束矢量K( θ,φ )=[cosθcos φ sin θ cos φ sin φ]T。
图1 半球共形阵机载气象雷达几何模型
如图1（b）所示，令右侧第一段圆弧为半球形阵列的第一层，每一层的第一个阵元位于X轴正方向，则由半球共形阵机载雷达阵列几何模型可以推导出第m(m = 1,2,…,M )行第n(n = 1,2,…,N)列的阵元位置坐标Pm,n为
式中，Rm 为第m 层圆弧的半径，
，β 为共形阵机载雷达圆弧所对应的圆心角，α为图1（c）YOZ平面所示的夹角。
根据阵元坐标以及波束矢量推导得到半球共形阵雷达地杂波回波信号的空间角频率wsc：
式中，φl 为第l(l = 1,2,…,L )个距离单元杂波回波的俯仰角，L为回波总距离门数。
半球共形阵雷达地杂波回波信号的时间角频率wtc为
式中，fr为脉冲重复频率，λ为波长。
假设半球形阵列天线的第m 行第n 列个阵元在第k个脉冲下接收的第l个距离单元的杂波回波数据为Cl(m,n,k)，则
式中，Rl 为第l 个距离单元的斜距，θ 的取值范围为[ 0,π ]。
1.2 低空风切变信号
由于低空风切变是体分布式目标，所包含的气象粒子具有微粒性、叠加性和随机性等特点，因此采用撒点法［13］对风切变场回波数据进行仿真，通过将单个距离单元内的所有散射点相加，得到该距离单元的风场回波信号。则半球形阵列天线的第m 行第n 列阵元在第k 个脉冲下接收的第l 个距离单元的风切变信号回波数据为Sl(m,n,k)，且有
式中，Q表示一个距离单元内散射点总数，Rq为第q个散射点的斜距，Aq 为由雷达方程推导得到的第q个散射点回波的幅度［14］：
式中，Z 为风场反射率因子，Pt 为机载雷达的发射功率，G为天线的增益。
式（5）中，wss(θq,φq )表示第q 个散射点的空间角频率，wts(θq,φq )为时间角频率，且有
式中，θq 为第q 个散射点的方位角，φq 为俯仰角，vq为径向速度。
1.3 回波信号
当半球共形阵机载气象雷达用于探测风切变时，其回波信号由地杂波C、风切变S 和噪声组成，则半球形阵列天线接收的第l 个距离单元的总回波数据Xl为
假设在雷达工作时，在一个CPI 内发射K 个脉冲，用Xl表示第l个距离单元的半球共形阵机载气象雷达所接收的回波数据矩阵：
2 基于D3AR 低空风切变风速估计方法
基于D3AR 低空风切变风速估计方法首先对待检测单元的回波数据分别从空域、时域、空时域进行信号对消处理，然后将处理后的样本数据通过AR 模型构造与杂波子空间正交的空间来实现对杂波的抑制，最后通过提取待检测单元的最大多普勒频率来估计风场速度。下面详细讨论以下3个部分：信号对消处理、估计STAR 模型系数矩阵以及低空风切变风速估计。
2.1 信号对消处理
DDD 方法仅利用待检测距离单元的数据信息进行样本数据的获取，当待检测距离单元含有低空风切变信号目标时，此时通过AR 模型来构造的空时滤波器在抑制杂波的同时会造成风切变信号损失。因此，为了防止风切变信号被抑制导致无法准确估计风速，首先需要对待检测距离单元数据分别从时域、空域以及空时域进行信号对消处理，具体过程如下：
定义风切变信号Sl(m,n,k) 为
根据载机速度、脉冲重复频率等先验信息，假设风切变信号方位已知，fl 为归一化多普勒频率（fl的范围为[]-1,1）。此时，对于低空风切变信号，相邻两脉冲间的相位差为
相邻两阵元之间的相位差为
式中，ΔPm,n 为相邻两阵元位置坐标差。因此待检测单元可通过利用两阵元（脉冲）之间的相位差分别从空域、时域以及空时域进行信号对消处理。
对第l 个待检测单元的和回波数据从空域利用两阵元相消滤除风切变信号时，可以得到如下数据矩阵：
式中，Yls为第l个待检测距离单元从空域利用两阵元相消滤除风切变信号后的数据矩阵，其维度为(MN - 1) × K，“*”表示共轭运算。
对第l 个待检测距离单元的回波数据从时域进行对消处理时，当待更新的fl 在[]-1,1 中取某一值时，利用两脉冲对消后得到如下数据矩阵：
式中，Ylt为第l个待检测距离单元从时域利用两脉冲相消进行信号滤除处理后的数据矩阵，其维度为MN ×(K - 1) 。
同理，利用两脉冲、两阵元相消对第l 个待检测距离单元回波数据从空时域进行信号对消处理后的数据矩阵为
式中，Ylst为第l个待检测距离单元从空时域利用两阵元两脉冲相消进行信号滤除处理后的数据矩阵，其维度为(MN - 1) ×(K - 1) 。
更新fl 的取值并执行式（13）～（15），当待检测距离单元利用信号对消后的样本数据矩阵经过AR 模型所构造的空时自适应处理器的输出功率最大时，此时Yls、Ylt和Ylst3 个数据矩阵中的风切变信号被有效消除。
2.2 估计STAR模型系数矩阵
假设当待更新的fl 在[-1,1 ]范围中所取某值是该待检测距离单元风场目标的归一化多普勒频率，此时由该单元通过信号对消处理后得到的3个样本数据矩阵Yls、Ylt以及Ylst中就只包含杂波和噪声，为了方便计算，选取Yls 中的(K - 1) 列数据、以及Ylt 中的(MN - 1) 行数据，再将3 个样本数据矩阵排列成(MN - 1 )(K - 1) × 1 维矢量并分别记为Y1、Y2 和Y3。由于任意随机过程都可描述为AR 模型，则Y1、Y2和Y3满足如下自回归模型［15］：
式中，t = 1,2,3，k = 1,2,…,K - P，Ai为D ×(MN - 1)维的空时自回归系数矩阵，D 代表AR 模型的空域阶数，P 为时域阶数。AR 模型的时域阶数和空域阶数未知，需要通过样本数据对其进行估计，由于已有许多文献介绍了AR 模型阶数估计方法［15-16］，本文不再叙述。模型阶数D、P 确定后，再通过训练样本估计自回归系数矩阵。
记A=[A0 A1…AD-1]，A的维数为D×(MN-1) P，将处理后的训练样本单元的数据重新改写为
式中，et 为(MN - 1) P ×(K - P )维矩阵，则处理后的所有样本数据矩阵为
式中，E为(MN - 1) P × 3(K - P )维矩阵，因此，根据式（17），训练样本数据矩阵E满足如下方程：
可以看出，式（20）中A的解是由矩阵E的零空间确定，但由于E 满秩不存在零空间，所以根据最小均方误差准则，对于自回归系数矩阵A 的求解可以转换为最小二乘准则来估计［16］：
为了确保不产生平凡解，A 必须满足AAH = I,且样本矩阵E 的D 个小奇异值对应的左奇异矢量构成矩阵A的列。
2.3 低空风切变风速估计
完成自回归系数矩阵A 的估计后，利用A0,A1,…,AD - 1构造如下滤波器系数矩阵：
式中，B 为D(K - P )×(MN - 1 )(K - 1) 维矩阵,且满足如下方程：
矩阵B 与杂波和噪声所张成的子空间正交，根据匹配滤波理论，得到D3AR 算法的空时自适应权矢量为［17］
式中，
为低空风切变风场信号空间导向矢量Ss中的前MN-1 个分量，
为时间导向矢量St中的前K-1个分量，且有
式中，ws 为归一化空间角频率，fl 上文中假设的归一化多普勒频率。由此可以得到信号匹配结果为
通过更新fl的取值，可以得到对应的空时自适应处理器权矢量，进行信号匹配后求解输出功率。当其最大时，此时fl的取值即为对应距离单元的多普勒频率估计值。这是因为当处理器的输出功率最大化时，该待检测距离单元经过信号对消处理后的样本数据矩阵中风切变信号被有效滤除，再通过AR 模型所构造的空时滤波器就可以使杂波信号得到有效抑制，同时风切变信号得到积累；如果处理器的输出功率不是最大，说明该待检测距离单元在信号对消处理过程中风切变信号没有被有效滤除，导致得到的样本数据经过AR 模型所构造的空时滤波器在抑制杂波的同时会造成风切变信号损失。即可以通过式（28）来估计得到该待检测距离单元
：
进而通过式（29）得到风速估计值
̂为
最后，依照上述方法分别估计每一个距离单元的风速值，完成整个风场的风速估计。
2.4 方法流程及步骤
图2展示了基于D3AR的半球共形阵风切变风速估计方法的整体流程。
图2 基于D3AR风速估计方法基本流程
步骤1：分别从空域、时域以及空时域利用两脉冲（阵元）进行信号对消处理；
步骤2：将处理后的训练样本单元描述为空时自回归模型；
步骤3：根据样本单元估计空时自回归模型阶数，并构造滤波器模型系数矩阵；
步骤4：根据匹配滤波理论，构造与杂波子空间正交的空间来实现对杂波的抑制，最后通过搜索该待检测距离单元的输出功率最大时所对应的多普勒频率，完成风速估计；
步骤5：更新待检测单元，重复上述步骤估计所有待检测距离单元的风速值，完成整个风场速度估计。
3 仿真结果及分析
3.1 仿真参数设计
假设风场目标处于8.5～16.5 km 处，表1 所示为仿真参数值。
表1 载机及雷达仿真参数
参数载机高度∕m信噪比∕dB主瓣方向∕（°）阵元行数距离分辨率∕m载机速度∕（m·s-1）参数值600 5 90 8 150 87.5参数圆心角∕（°）杂噪比∕dB波长∕m阵元列数采样脉冲数脉冲重复频率∕Hz参数值（90,0）40 0.05 24 64 7 000
3.2 仿真结果分析
图3 为半球共形阵机载气象雷达回波信号第87 号距离门空时谱仿真结果，由图可以看出，地杂波信号的空时二维谱在空间的分布表现为椭圆形状，风切变信号在空间表现为一条窄带，且地杂波信号的功率明显大于风切变信号，使得风切变信号完全被淹没，从而导致难以被检测。
图3 第87号距离单元空时谱
图4 为D3AR 方法、最优STAP 方法、DW-STAP方法以及STAR 方法的第87 号距离单元的改善因子对比图。从图可以看出，上述方法均在主杂波区域即零频区形成凹口，造成性能损失，与其余方法相比，本文方法在主杂波区域形成的凹口更窄、更浅，杂波抑制性能明显提高。
图4 改善因子对比图
图5 为D3AR 方法与最优STAP 方法、DWSTAP 方法、以及STAR 方法风速估计结果对比图。从图中可以看出，本文所提方法在不利用参考单元样本数据的情况下仍然能较准确地估计风速。对比图显示，最优STAP 法的风速估计结果不理想，主要是由于半球共形阵特殊阵列流型导致杂波非平稳性更强，使得所需的独立同分布样本数不满足RMB 准则，最优STAP 方法性能受损［5］；经DW补偿后，仍然不能得到理想效果，这是由于DW方法仅对时域进行了多普勒频率补偿，没有对空间频率进行补偿，所以经过DW 补偿后，仍不能准确估计风速［18］；STAR 方法需要用到与待检测单元相邻单元的数据来构造AR 模型建立空时二维滤波器进行杂波抑制以及估计风速［17］，其结果的准确性依赖于IID 样本数量。而D3AR 方法仅利用待检测单元的数据且不需要进行距离依赖性矫正就能准确估计风速。
图5 风速估计结果
本文所提方法主要包括求滤波器权矢量和空时滤波两部分的运算量，表2 为D3AR 方法、最优STAP 方法、DW-STAP 方法以及STAR 方法运算量对比，其中L1 为最优STAP 方法和DW-STAP 方法中所需的训练样本数。从表中可以看出本文所提方法的运算量比最优STAP方法、DW-STAP法以及STAR法都要小。
表2 不同风速估计方法运算量对比
方法最优STAP DW-STAP STAR D3AR滤波器权矢量()MNK3 +()MNK2L1()MNK3 +()MNK2L1()MNP3 +()MNK3D()K - P - 1()P()MN-13+()MN-12( )K-12D( )K-P空时滤波MNK MNK MNK()MN-1 ()K-1
表3 给出了D3AR 方法、最优STAP 方法、DWSTAP 方法以及STAR 方法的均方根误差对比，可以看出，本文所提D3AR方法的均方根误差最小。
表3 不同风速估计方法均方根误差
风速估计方法最优STAP DW-STAP STAR D3AR均方根误差∕（m·s-1）18.180 9 14.908 1 7.704 7 1.143 0
4 结束语
针对半球共形阵体制下进行低空风切变检测时会受到强低杂波信号的干扰，导致风切变信号难以检测的问题，提出了一种基于D3AR 的半球共形阵风切变风速估计方法。该方法将空时自回归算法思想引入到直接数据算法思想中，仅用待检测单元滤除风切变信号后的样本数据来构建空时自回归模型，再通过构造与杂波子空间正交的空间来实现对杂波的抑制，最后通过提取待检测单元的最大多普勒频率来估计风场速度。仿真结果表明，该方法可以在半球共形阵体制下仅用待检测距离单元的数据准确估计风速，且该方法仅在时域进行滑窗处理，避免了空域滑窗对阵列结构的严格要求，因此适用于任意阵列，但由于该方法中AR 模型参数的选择会影响风速估计结果的准确性，因此需进一步研究更为准确的AR 模型参数估计算法。
［1］DESHPANDE M D, STATON L. Determination of Windspeed Within a Weather Storm Using Airborne Doppler Radar［C］∕∕IEEE Proceedings of the Southeastcon,Williamsburg,VA,USA:IEEE,1991:508-519.
［2］WILSON J W, WAKIMOTO R M. The Discovery of the Downburst: T. T. Fujita's Contribution［J］. Bulletin of the American Meteorological Society,2001,82（1）:49-62.
［3］李海，谢瑞杰，谢伶莉，等.复杂地形环境下基于GAMPSTAP 的低空风切变风速估计方法［J］.电子与信息学报,2023,45（2）:576-584.
［4］李海,毕金枝,孟凡旺,等.机载柱形共形阵低空风切变风速估计方法［J］.雷达科学与技术,2022,20（6）:651-657.
［5］呼延泽.基于毫米波LFMCW 雷达的低空风切变解模糊方法研究［D］.天津：中国民航大学，2021.
［6］RTCA∕DO-220. Minimum Operational Performance Standards for Airborne Weather Radar with Forward-Looking Windshear Capability［S］.Washington:RTCA,2016.
［7］张泽鹏.共形阵列稀布与方向图综合算法研究［D］.南京:南京理工大学,2021.
［8］李海，谢伶莉，张志宁.复杂地形环境下基于CL-KASTAP 的低空风切变风速估计方法［J］.电子与信息学报，2021，43（12）:3647-3655.
［9］LI Hai, FENG Kaihong, CHENG Weijie ,et al. Estimation of Low-Altitude Wind-Shear Wind Speed Based on DDDEFA Under Aircraft Yawing［C］∕∕2021 CIE International Conference on Radar, Haikou, China:IEEE, 2021:3106-3111.
［10］李海，程伟杰，谢瑞杰.基于同伦稀疏STAP 的低空风切变风速估计［J］.系统工程与电子技术，2022，44（4）:1174-1181.
［11］LI Hai, YANG Wenheng, HU Yanze, et al. Low-Altitude Wind-Shear Wind Speed Estimation Based on Improved Auxiliary Channel［C］∕∕2021 CIE International Conference on Radar,Haikou,China:IEEE,2021:3112-3116.
［12］WARD J. Space-Time Adaptive Processing for Airborne Radar Data Systems［R］. Lexington, Massachusetts: Lincoln Laboratory of MIT,1994:25-45.
［13］焦中生，沈超铃，张云.气象雷达原理［M］.北京：气象出版社，2005.
［14］李静，何姣阳，李田家，等.基于FLUENT 的下击暴流三维风场建模［J］.成都信息工程大学学报，2021，36（5）:508-511.
［15］DUAN Keqing, XIE Wenchong, WANG Yongliang. A New STAP Method for Nonhomogeneous Clutter Environment［C］∕∕2010 the 2nd International Conference on Industrial Mechatronics and Automation, Wuhan, China:IEEE,2010:66-70.
［16］DUAN Keqing,FEI Gao, XIE Wenchong ,et al. Parameters Estimation Method for STAR Based on Clutter Degree of Freedom［C］∕∕2009 WRI Global Congress on Intelligent Systems,Xiamen,China:IEEE,2009:211-215.
［17］PARKER P,SWINDLEHURST A. Space-Time Autoregressive Filteringfor Matched Subspace STAP［J］.IEEE Trans on Aerospace and Electronic Systems,2003,39（2）:510-520.
［18］魏小丹.非均匀环境下机载雷达空时自适应处理方法研究［D］.西安:西安电子科技大学，2022.