在气象环境预报保障和大气科学研究中,准确描述和理解复杂大气流体运动结构变化和不稳定机理至关重要(穆穆等,2003)。然而,近几十年来,即使随着计算机技术的高速发展和气象观探测仪器性能不断提升,要准确量化地提取复杂大气流体运动结构依然非常困难。另一方面,在对大气流动进行稳定性分析时,需要对控制大气运动的偏微分方程组进行大幅度简化或者计算大型矩阵的逆和特征值,众所周知大型矩阵计算代价巨大,且在大气运动稳定性分析中不存在普适性方法。而DMD(Dynamic Mode Decomposition)方法作为从整体稳定性分析(Mezi$\acute{c}$,2005;Rowley et al.,2009)的基础上发展起来的一种低维系统分解技术,不需要依赖任何具体的基函数,只需要一组时空数据就能分解得到流场多尺度运动模态的空间结构及对应频率(Schmid and Sesterhenn,2008);相比于现有的大数据处理算法,DMD方法具有模型简单、数值稳定性良好、计算效率高等特点,能够给出准确的分析结果并为非线性动力系统提供简化的线性模型(Tirunagari et al.,2017;Mohan et al.,2018)。

DMD方法由Schmid and Sesterhenn(2008)提出,基本思想是直接从实验或数值模拟获得的流场数据中提取出流场的动力学信息。Schmid(2010)对DMD方法在数学上进行了详实的理论推导,并将其分别对在数值方法和实验得到的流场数据进行分解,从中提取动力学模态信息,验证了DMD方法的可靠性。Muld et al.(2012)应用DMD方法对高速列车周围流场的数值模拟结果进行处理,提取了关键性流场涡团结构。Tu et al.(2014)等介绍了DMD与其他线性系统辨识方法的关联和等价性。Rowley and Dawson(2017)综述了DMD方法在模型降阶、流场分析和控制等问题上的应用。潘翀等(2010)率先开展了对DMD方法的应用研究,将DMD方法应用在NACA0015翼型加装Gurney襟翼后的尾流流场分析中。研究表明,DMD方法能够捕获流场中的各阶高频谐波,为研究复杂流场中不同尺度的相干结构提供了重要信息。张宾等(2015)对旋转圆柱绕流的实验数据进行了DMD处理,发现了旋转圆柱与尾流结构共振表现出低频涡结构特征,对相关工程领域避免事故的发生具有重要的指导意义。叶坤等(2016)采用DMD方法对超燃冲压发动机凹腔流动的稳定性进行分析,进一步理解了凹腔自激震荡的机理。综上所述,作为一种在不同领域广泛应用的模态分解方法,DMD方法的优势在于不仅能够获得流场的动力学信息,还可以给出流场特征结构的频率和增长率等信息(秦文瑾等,2016;邵昱昌和王彤,2017)。DMD方法适用于分析复杂流体运动(洪树立和黄国平,2017),同时在不同领域动力学系统建模方面具有广阔应用前景(寇家庆和张伟伟,2018;Tirunagari et al.,2017)。但是,将动力学模态分解方法应用于大气运动流场分析,尤其是平流层急流类型的天气系统分析还未见相关的公开报道(潘留杰等,2018;许琪等,2018)。

本研究首次将动力学模态分解方法应用于大气运动流场数据分析中。基于ERA5再分析数据集,首先采用DMD方法对一段时间内200 hPa急流运动速度场进行分解,提取了流场的主要模态以及频率,分析了急流运动在时间和空间上的主要特征,且根据特征值获得的模态增长/衰减率判断出相应模态的稳定性。其次建立了200 hPa大气急流运动演化的动力学降阶模型,以重构和预测流场动力学过程,评估了DMD方法在分析大气运动流场特征的能力。

1 资料和方法

1.1 DMD方法基本原理

DMD方法是基于Koopman理论(Koopman,1931)发展起来的一种数据驱动的矩阵分解方法,结合了主成分分析(Shlens,2005)和傅立叶变换的优点(Tirunagari et al.,2017),能够同时提取出数据序列的空间模态及其对应的频率。

对于任意给定的动力系统以恒定采样时间间隔Δ

其中 $\mathrm{ }{\mathit{v}}_j$代表第j个时刻的流场状态值,例如某个时刻的大气运动矢量场;$\mathrm{ }{\mathit{V}}_1^N$表示由待分析的所有时刻流场快照构成的矩阵,下标1表示序列中的第一个元素,下标N表示进入序列中的最后一个元素,即矩阵${\mathit{V}}_1^N$的第一列和最后一列。通常对向量${\mathit{v}}_j$在M个空间网格点进行模拟计算,即${\mathit{v}}_j\in {\mathbf{R}}^M$。

根据式(1)中的快照序列数据集可以构造出另外两个不同的快照矩阵${\mathit{V}}_1^{N-1}$,$\mathrm{ }{\mathit{V}}_2^N$

其中 ${\mathit{V}}_1^{N-1}$为包含从1到第N-1个流场的序列集合;矩阵${\mathit{V}}_2^N$为包含从2到第N个流场的序列集合。对于获取的大气运动流场数据,通常情况下空间网格点数目要远大于时间快照数目,即$M\gg N$。

考虑到流场变化的连续性,假定流场后一时刻的速度场能用前一时刻的速度场进行近似线性(由于大气运动是缓变的,因此在短时间内满足线性变化假设)表示,即:

由此得到快照矩阵

当快照数目增加到足够多,可以假定流场快照之间是线性独立的。由此快照矩阵的最后一个向量${\mathit{v}}_N$可以表示为前N-1个向量的线性组合,即:

表示为矩阵形式

其中 $\mathit{ }\mathit{c}=\left(c_1{c}_2\cdots c_{N-1}\right)$表示线性相关系数;R为残差向量;$\mathrm{ }{\mathit{e}}_{N-1}^{}$为单位向量。

由以上可得

上式中矩阵S为酉矩阵,且仅有c矩阵未知。

当残差最小时,酉矩阵S的特征值近似表达了矩阵A的特征值,从而可通过分析低维矩阵S的特征值和特征向量来研究原高维度系统的时空特征。矩阵S可通过对流场快照矩阵${\mathit{V}}_1^{N-1}$进行QR分解获得,但该方法容易产生病态矩阵,无法准确地提取流场的模态。因此Schmid and Sesterhenn(2008)提出通过奇异值分解对高阶算子进行相似变换,得到系统的低阶表达。该方法具有更好的数值稳定性,现已成为动力学模态分解的标准算法。

1.2 大气运动流场的DMD分析算法

DMD算法数学表达式简单,计算易于实现,为许多流体力学问题的机理分析提供了新的有力工具。本文将DMD方法应用于大气运动流场分析,无论是高分辨率数值预报模式产生的预报场数据,还是全球变分资料同化系统产生的分析场都只是描述了某一个时刻大气运动的状态,即严格来说是瞬时状态场(实际该时刻前后状态变化不大),相当于对大气运动状态进行相机拍照后存储下来,因此可以借用流体力学中的流动快照概念,命名为大气运动快照。本文实际上是利用DMD方法对一段时间200 hPa大气急流运动快照进行分析,采用标准DMD算法实现大气运动流场的动力学模态分解,详细说明如下:

1)首先从大气二维运动矢量场V(u,v)中获取在时间上连续均匀分布的N个时刻的数据序列,构建快照矩阵X。其中u、v为水平速度分量:

其中 $\mathrm{ }{\mathit{x}}_j$为第j个时刻的大气运动矢量场;$\mathrm{ }{\mathit{x}}_j$和$x_{j+1}$之间的时间步长为Δt;向量${\mathit{u}}_M^N$中上标代表第N个时刻,下标代表第M个空间网格点。在标准DMD算法中要求快照间时间间隔必须相等。对于三维速度矢量场,可以将垂直速度分量按相同的规则添加到上述公式中。

2)快照矩阵${\mathit{X}}_1^{N-1}$包含矩阵X的前N个时间序列,矩阵${\mathit{X}}_2^N$则由矩阵X的后N个时间序列构成:

根据DMD方法的原理,通过系统矩阵A,快照矩阵${\mathit{X}}_2^N$可以表示为:

3)矩阵A是一个含有数据量极大的矩阵,难以直接计算。因此可以给出矩阵A的相似矩阵$\widehat{A}$来替代高维矩阵A。为寻求相似变换的正交子空间,可通过对矩阵${\mathit{X}}_1^{N-1}$进行奇异值分解:

奇异值分解得到的矩阵U和W为酉矩阵,$\mathrm{ }{\mathit{W}}^{\mathrm{T}}$是矩阵W的转置矩阵,$\mathrm{ }{U}^{\mathrm{T}}$是矩阵U的转置矩阵。其中$\mathit{U}\in {\mathit{R}}^{M\times r}$,包含数据序列的空间结构;$\mathbf{ }\mathbf{\Sigma }\in {\mathbf{R}}^{r\times r}$,为奇异值的对角矩阵;$\mathit{ }\mathit{W}\in {\mathbf{R}}^{N\times r}$,描述了数据的时间结构;$\mathrm{ }\widehat{A}\in {\mathbf{R}}^{r\times r}$,r为矩阵${\mathit{V}}_1^{N-1}$的秩。

4)矩阵$\widehat{A}$可通过最小化${\mathit{X}}_2^N$和$A{X}_1^{N-1}$差值的Frobenius范数来确定:

结合式(14)和(15),式(16)可表示为:

因此,矩阵A的最优估计$\widehat{\mathit{A}}$可以表示为:

由于矩阵$\widehat{\mathit{A}}$是A的相似变化,因此矩阵$\widehat{\mathit{A}}$包含矩阵A的主要特征值。记矩阵$\widehat{\mathit{A}}$的第j个特征值为${\lambda }_j$,对应的特征向量为$w_j$,则第j个DMD模态为${\mathrm{\Phi }}_j=U{w}_j$。

5)根据提取出的大气运动流场对应的DMD模态,可以进一步估计大气运动流场的演化过程。将DMD模态作为基,原始快照矩阵可展开为:

其中 $\mathrm{ }{\omega }_j$对应矩阵$\widehat{\mathit{A}}$的第j个特征值${\lambda }_j$,$\mathrm{ }{\omega }_j$=$\mathrm{l}\mathrm{o}\mathrm{g}\left({\lambda }_j\right)/\mathrm{\Delta }t;\mathit{\alpha }\left\{{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_r\right\}$为DMD模态对应的幅值,可通过最小二乘法求解方程$x_1=\mathrm{\Phi }\alpha$求得:

其中 $\mathbf{ }\mathbf{\Phi }\left\{{\mathrm{\Phi }}_1{\mathrm{\Phi }}_1\cdots {\mathrm{\Phi }}_r\right\}$为快照矩阵${\mathit{X}}_1^{N-1}$对应的DMD模态矩阵;向量${\mathit{x}}_1$为快照矩阵${\mathit{X}}_1^{N-1}$的第一列。由于式(19)在时间和模态上都是离散的,可通过外推对大气运动矢量场采样时间外的流场动态发展过程进行近似预测。实际应用中,选择前几阶模态即可实现对大气运动流场的准确重构和预测。

以上,为对实现大气运动流场动力学模态分解算法的简单描述。获得的最重要的结果是模态${\mathrm{\Phi }}_j$和其对应的特征值${\lambda }_j$。$\mathrm{l}\mathrm{o}\mathrm{g}\left({\lambda }_j\right)/\mathrm{\Delta }t$的实部对应模态的增长率,表示模态随时间的演化和发展情况,即是增长还是衰减;虚部对应模态的频率。在进行稳定性判断时,如果增长率为正,则对应的模态不稳定;如果增长率为负,则对应的模态为稳定模态;若增长率为零则对应的模态为周期性模态。同样也可以将特征值置于复平面上进行判断,对于中性稳定的流动,大部分特征值都落在复平面上距离原点幅值为1的单位圆附近,对应的模态为稳定模态;位于单位外的特征值对应的模态是不稳定模态;位于单位圆上的特征值对应的模态是周期性模态。

2 数值实验结果及分析

使用的资料为欧洲中期天气预报中心(European Center for Medium-Range Weather Forecasts,ECMWF)ERA5水平风场高分辨率再分析资料,资料的水平分辨率为0.25°×0.25°,时间分辨率为1 h(赵天保等,2010;高庆九等,2018),只分析垂直方向高度为200 hPa的平流层急流流场,提取流场的时间段为:2013年1月7日00时—12日23时,瞬态流场数量共计为144,分析区域为泛亚洲区(30°~150°E,25°~50°N)。

2.1 结果与分析

图1为200 hPa急流运动流场DMD模态的特征值的分布情况。几乎所有的特征值都落在单位圆附近,这表明DMD分解的200 hPa急流速度场数据是位于吸引集周围的,也说明DMD分解得到的模态是中性稳定的(Hemati et al.,2014)。由于处理的数据序列为实数,图1中特征值呈共轭对称分布,对称中心对应的频率为零。图1a中,横轴为实部,纵轴为虚部。空心点表示在单位圆内,对应的模态为稳定模态;实心点表示在单位圆外,对应不稳定模态。可以看出绝大部分模态为稳定模态,且其对应的特征值靠近单位圆,不稳定模态极少。特征值取对数后的分布(图1b)。其中,横轴为虚部,纵轴为实部。空心点对应的模态为稳定模态,实心点对应的模态为不稳定模态。稳定模态可以理解为对应该时间段内变化缓慢、逐渐消退的天气过程;不稳定模态则对应该时间段内逐渐发展增强的大气流场结构,对天气系统的变化发展起主导作用(苗春生等,2018;张柳等,2020)。

图2为200 hPa急流运动流场DMD模态幅值与对应频率和增长率之间的关系。其中频率和增长率均经过归一化处理。根据式(20),模态幅值α代表了模态对初始时刻流场$x_1$的贡献。在物理意义

上,幅值α的绝对值大小代表了各个模态携带着的能量/强度。图2a中频率对应的模态幅值存在一个明显的峰值,随着模态频率的增大对应的幅值逐渐减小,即高频模态对初始时刻流场的贡献逐渐减小。图2b中,空心点对应的增长率为负,对应的模态衰减;实心点对应的增长率为正,对应模态增长。可以看出绝大部分模态的增长率为负,并集中在零附近,表明该部分模态随时间缓慢衰减;极少数模态增长率为正。

表1为200 hPa急流运动流场前6阶模态的增长率以及频率。其中第1阶模态增长率较小,频率为零,其模态代表了流场中与时间无关的信息,近似于平均流场。第2~6阶模态为共轭模态,其中第3、5阶模态对应的增长率为正,为不稳定模态,对流场的变化发展起主导作用。其余模态为稳定模态,在增长率上略有差异,频率较小,呈缓慢衰减趋势。

DMD方法能够分解出流场不同z物理频率代表的运动模态。在空间上低频模态表示了较大尺度的运动特征、高频模态表示的是小尺度特征。不同的模态代表了流场不同空间尺度的运动规律,这些流动模态叠加就形成了流场所体现出的主要特征。不同尺度的DMD模态与大气科学中的大、中、小尺度没有严格的对应关系,后者对应不同空间尺度的具体天气系统。本文初次尝试应用DMD方法对大气运动流场进行分解,仅对一个区域200 hPa的急流运动流场数据集进行分析,当利用DMD方法对全球大气运动场数据进行分析时,得到的不同动力模态特征可能反映了全球大、中、小尺度大气运动的特征,但考虑到大气多尺度运动的复杂性,分解得到的不同大气动力模态与全球大、中、小尺度大气运动不可能具有完全的一一对应关系。

200 hPa急流运动流场前六阶DMD模态实部的空间形态分布如图3所示,模态虚部并未展示,虚部与实部存在90°的位相差,模态均经过归一化处理。图中伪彩色图反映流场垂直涡度的分布,箭头代表模态流场的速度矢量,特征值的实部为增长率,虚部为频率。总的来说,获得的六个模态分别展示了200 hPa急流的多尺度运动特征,具体分析如下。对于DMD模态而言,第一阶模态对应的增长率为负,频率为零,故其对应的DMD模态只有实部。该模态流场结构包含时间平均的流场信息,展现了急流流场中占主导地位的流动结构。如图3a所示,第一阶模态流场的空间分布呈现出波动状流动结构,以纬向西风为主,这种特征在东西方向上呈一致的

延伸,在80°E、25°~40°N以东区域存在一支较强的西风急流,急流轴位于25°~35°N附近,第一阶模态的空间分布特征基本符合我们对西风急流空间结构的认识。第二阶模态实部也为负,为衰减模态,相比第一阶模态,第二阶模态在流场大尺度运动的背景上出现较大尺度的涡旋结构,并呈现出气旋、反气旋交替分布的特点。需要注意的是模态中涡的旋转方向并不代表流场中真实涡的旋转方向,模态的空间分布结构只是体现动力学模态分解对急流运动流场运动规律的提取,这些结构无法从大气运动瞬态流场中直接体现出来,因而DMD方法获得的流场结构具有非直观性。第三阶模态的空间分布与第二阶模态类似,涡旋结构明显,但其特征值实部为正,为不稳定发展模态。不稳定模态对流场从稳定平衡状态到不稳定发展状态起主导作用,具有重要意义。第四阶模态中,随着频率的增加涡旋结构的尺度明显减小、数量增多,空间结构更为复杂,其对应的增长率为负,为衰减模态。第五阶模态空间分布同第四阶模态类似,但其增长率为正,为不稳定增长模态。第六阶模态的频率更高,流场的结构尺度进一步减小,空间分布也更为杂乱。在物理上这些高频模态代表了嵌入大尺度流场涡旋结构这一动力过程中的典型各阶小尺度扰动,这些模态可能源于大尺度涡旋流动对环境流体的卷携作用,它们虽然在相干性上不及大尺度模态,但刻画了流场的扩张。

由上述分析可得,DMD方法可以通过对一段时间大气流场数据序列的分析,提取出该时间段内流场不同频率、不同尺度的相干结构,这些结构无法从大气运动瞬态流场中直接体现出来,但是其对大气流场内部能量的转移、大气系统的不稳定性都有很大的影响。因此,将DMD方法应用在大气运动分析上,有助于对大气运动特征有更进一步的理解。

2.2200 hPa急流流场的重构和预测

为了进一步观察DMD方法对大气200 hPa急流运动流场特征的提取效果,基于动力模态分解获得的200 hPa急流运动流场的前六阶模态,根据式(19)建立大气急流运动流场的降阶模型,对流场进行重构。图中伪彩色图反映流场垂直涡度的分布,箭头代表急流运动流场的速度矢量。如图4所示,左侧一列分别T=12 h、60 h、120 h的真实流场,右侧为相应时刻的重构流场。可以看出三个时刻重构流场与真实流场的空间结构基本吻合,充分说明DMD方法能够准确地提取出大气流场的主要结构,并且分解得到的前几阶模态就包含了原始流场的大部分信息,从而在给定的时间段内能够给出原始流场的准确重构。

根据式(19),利用DMD模态对200 hPa急流运动流场重构后的低维表示可以对大气动力学流场的动态发展过程进行近似预测,预测流场的时间范围为2013年1月13日—16日,共96 h。图5给出了采样时间段和预测时间段模态系数(模态振幅/能量)随时间的变化,由于第一阶模态增长率和频率都近似为零,因此没有给出。其中图5a为前6阶模态中稳定模态的模态系数随时间的变化情况,可以看出稳定模态的模态系数随时间均表现为趋于收敛,这与表1中对应模态增长率为负相对应;由表1可得稳定模态增长率的绝对值为DM₂>DM₄>DM₆,对应图中可以看出第二阶模态系数的衰减速度最快,其次为第四阶、第六阶模态。图5b为不稳定模态的模态系数随时间的变化情况,可以看出不稳定

模态的模态系数缓慢地趋于发散,呈现出增长趋势;由表1可得不稳定模态的增长率绝对值为DM₃>DM₅,对应图中可以看出第三阶模态随时间的增长速率明显大于第五阶模态。

3 讨论和结论

介绍了动力学模态分解方法(DMD)的基本理论及算法,提出在大气科学数据中引入DMD方法进行大气运动流场的多尺度分析,以200 hPa急流运动流场为研究对象,引入DMD方法进行分析,得出如下结论:

1)DMD方法成功捕捉到了200 hPa急流天气系统运动变化过程中的主要模态和对应频率以及流场在不同频率下的流动特点,包括不同频率下模态的空间分布、能量/强度的分布及其稳定性;并基于获得的前6阶模态建立了急流流场的降阶模型,实现了流场的准确重构。

2)相比传统的模态分解方法,DMD方法可以有效获得多尺度大气运动中每个频率运动模态的结构信息,实现了从复杂大气运动流场中提取单一频率拟序结构的能力。

3)通过DMD分解方法可以获得大气运动流场不同模态对应的特征值,模态特征值的实部反映流场能量/强度在时间层面的增长和衰减程度,特征值的虚部可以反映流场结构在时间层面的形态变化,通过特征值分析即可判断大气运动模态的稳定性。

4)DMD方法分解得到的不同的模态代表了流场不同空间尺度的运动规律,零频模态代表了大气运动流场的基本流动结构,大气运动流场可以看作在这个基本流动结构上叠加不同频率的动力学模态,由此,将复杂的流动问题转化为若干主要流动模态随时间的演化。

参考文献