1 引言

自Rossby(1940)和Ertel(1942)提出位涡概念以来，位涡一直受到中外气象学家的高度重视，并在大气科学众多领域得到广泛应用，涉及从边界层过程(Zhang, et al, 2009)、小尺度地形环流(Rotunno, et al, 1999)及雷暴结构(Gray, 2001; Gold et al, 2008a, 2008b, 2008c; Nielsen-Gammon, et al, 2008b; Russell et al, 2008, 2012; Chagnon, et al, 2009; Weijenborg, et al, 2017)、中尺度锋生动力学(Thorpe, 1990; Davies, et al, 1998; Chen, et al, 2003)到天气尺度扰动的形成与发展(Wu, et al, 2000; Romero, 2001; Lackmann, 2002; van Delden, 2003; Ahmadi-Givi, 2004)、大尺度的极涡变化(易明建等, 2009)，乃至统计气候(赵亮等, 2009; Čampa, et al, 2012)等各个方面。从位涡的守恒性和可反演性出发，到等熵位涡思想(寿绍文, 2010; 陶祖钰等, 2012; 周小刚等, 2014)、分段位涡反演等一系列位涡反演方法(赵玉春等, 2007; Nielsen-Gammon, et al, 2008a)，再到位涡方程及位涡收支分析(高守亭等, 2002; 赵兵科等, 2008; Tory, et al, 2012; 吴迪等, 2015)、湿位涡(吴国雄等, 1995)及其在强降水天气过程中的应用(王建中等, 1996; 任余龙等, 2007; 葛晶晶等, 2012; 张研等, 2015)等，与位涡有关的研究成果越来越多并逐步开始向指导业务方面转化(Morgan, et al, 1998; Brennan, et al, 2008)。特别是随着数值天气预报模式的逐渐普及，高时、空分辨率的模式输出资料为进一步开展与位涡有关的理论研究与业务应用提供了便利条件(黄亿等, 2009; 赵小平等, 2014)。而位涡的多尺度研究与应用又对其计算精度提出了更高的要求。

很多数值天气预报模式(如WRF、ARPS)的垂直结构采用地形追随坐标系，其并非传统定义位涡的z坐标系或p坐标系，因此，在利用模式输出资料计算位涡时，需要考虑坐标系转换问题及其带来的计算误差。比较普遍的计算方式是将模式坐标系格点上的有关变量插值到z坐标或p坐标之中再做计算(下简称为插值法)，但是由于风场和位温场的量级及其空间分布不同，特别是在台风、温带气旋、锋面、中尺度对流系统等重点关注的天气系统附近，气象要素水平梯度往往比较大，并经常伴随着较强的垂直切变，且由于中国的三级阶梯地势和复杂的地形环境，在地表海拔高度起伏变化大的地区，地形追随坐标系与z(p)坐标系的贴合程度就会变得较差，受这两点因素的共同影响，插值过程会给计算结果引入更多的误差。位涡的计算误差会直接对其在中小尺度对流性天气精细化结构分析中的应用产生影响，此外，误差随着时空不断累积，将会增强有关统计分析的不确定性，对结果产生不良影响。Cao等(2011)利用坐标变换方程从p坐标系静力位涡即静力平衡条件下埃特尔(Ertel)位涡的简化形式出发，推导了其在地形追随坐标系下的表达形式。而后直接通过绝热条件下的f平面地形追随坐标系原始方程组，导出了相同的计算形式，并证明了这一位涡形式亦满足地形追随坐标系中的守恒性质，为改善位涡计算误差问题提供了一种新的思路，即通过动力约束的方法，可以直接利用模式的地形追随坐标系输出资料计算p坐标系静力位涡，从而减少插值误差，但其推导的位涡形式并未全面考虑到垂直运动的影响。而随着天气研究向更小尺度发展，同时考虑三维风场和位温分布的全型埃特尔位涡能够更加完整反映中小尺度天气系统的结构特征，在今后的业务及研究中存在着更为广泛的应用前景。

2 计算方法推导 2.1 地形追随坐标系与z(p)坐标系的变换

任一气象要素可表示为三维空间变量和时间的函数

(1)

式中，垂直坐标可选为气象中常用的z坐标或p坐标。而z(p)坐标又可在地形追随坐标系中表示为

(2)

式中，η为地形追随坐标系垂直坐标。例如η在中尺度WRF模式中的定义为η=(p_h-p_t)/(p_s-p_t)，其中，p_h为静力气压，p_s和p_t分别为地面气压和顶层气压，容易看出地面η值为1，并随高度递减，顶层η值为0。

利用复合函数求导法则，由式(1)、(2)分别求取气象要素A对地形追随坐标系变量x、y、η、t的偏导，并整理后可得

(3)

(4)

(5)

(6)

以上4式即为地形追随坐标系与z(p)坐标系的变换关系，可以将式(3)(对z坐标)的变换关系用图 1形象地表示。图中δA_z=A(G)-A(O)，δA_η=A(M)-A(O)，两边同时除以δx可得，取极限δx→0即可得式(3)，式(4)—(6)均可采用类似方式进行形象理解，不再赘述。

通过以上推导得出的地形追随坐标系与z(p)坐标系的变换关系，可以很方便地将z(p)坐标系中定义的物理量改写为地形追随坐标系中的表达形式，从而直接计算其在模式网格点上的准确值。下面即利用这一变换关系写出地形追随坐标系下直接计算z坐标系埃特尔位涡的表达式。

2.2 埃特尔位涡在地形追随坐标系的表达形式

z坐标系中埃特尔位涡可表示为

(7)

式中，ζ_αz为z坐标系三维绝对涡度矢量，▽_zθ为z坐标系三维位温梯度矢量，其展开式为

将展开式中风场与位温x、y方向上的偏导数用式(3)、(4)替换，可得

式中，。从模式中容易得到地形追随坐标系格点的高度信息，故为计算方便考虑，不必对z方向各项进行变换。

因此，z坐标系埃特尔位涡在地形追随坐标系中可简单表示为

(8)

式中，ζ_αη^*由于在z方向上并未做变换，可称之为伪地形追随坐标系三维绝对涡度矢量；▽_η^*θ可称之为伪地形追随坐标系三维位温梯度矢量(z方向未做变换)；C为地形修正项。

综上容易看出，在利用垂直方向为地形追随坐标的模式输出资料对埃特尔位涡进行计算时，若欲采用式(7)(即插值法)，则必须将模式输出资料预先插值到z坐标系上，无法避免地带来了前文所述的误差问题；而式(8)(即直接法)则不必考虑这一预插值过程，及其带来的误差，直接得到模式格点上的准确值，并且式(8)与式(7)在物理意义上是等价的，便于利用位涡的一系列性质作进一步研究与应用。下面将通过地形追随坐标系数值模式模拟的实际个例采用以上两种方法分别计算埃特尔位涡，检验直接法的优越性，并对插值法所带来误差的分布进行分析。

3 个例数值模拟与计算效果检验 3.1 资料与研究方法

分别采用插值法(即式(7))和直接法(即式(8))利用模式输出的风场和位温资料对2016年7月18—21日发生在华北地区的一次黄淮气旋爆发过程中强降水时段20日01时(北京时，下同)的埃特尔位涡进行了计算。地形追随坐标系格点的风场和位温资料来自中尺度WRF模式。模拟采用3重嵌套，水平网格距分别为27、9、3 km，垂直方向分为51个η层，自下而上分层由密至疏，顶层气压为50 hPa。初始场和侧边界条件使用NCEP的0.5°×0.5° GFS资料，采用MYJ边界层方案、Noah陆面过程方案，微物理过程采用Lin方案，在27和9 km网格采用Kain-Fritsch积云对流方案。模式积分时间为2016年7月19日02时—20日20时，共42 h，3个网格结果输出间隔分别为60、20、5 min。

研究所关注的区域，即模式3 km网格区域，其东西向格点数为211，南北向格点数为223，格点上的地形数据由WRF中提供的最高分辨率30″(约0.9 km)地形资料插值而来。具体范围如图 2所示，包括京津冀地区、山西东部、河南北部和山东西北部地区，从西向东分布着黄土高原、太行山脉、华北平原及山东丘陵，地形背景复杂多变，地形的起伏特征及其梯度变化在3 km网格中得到了较好的反映。20日01时，黄淮气旋中心位于冀鲁豫三省交界处附近，中心最低气压990.6 hPa，气旋中心附近气象要素空间变化较为剧烈，与起伏的地形影响叠加在一起，为插值法计算的埃特尔位涡引入了更多误差，可满足研究所需。

首先利用式(7)(插值法)和式(8)(直接法)通过中央差分分别计算z坐标系格点和WRF-η坐标系格点上的埃特尔位涡，其中采用插值法计算时，需预先将WRF模式输出的η坐标系上的风场和位温场资料线性插值(视其为计算误差的主要来源)到z坐标系(这里采用的z坐标系垂直共分98层，两层间隔为恒定200 m)。之后为了方便进行误差分析，分别将插值法算得的z坐标格点上的埃特尔位涡(以下简写为PV_z)插值到η坐标系，并将直接法算得的WRF-η坐标格点上的埃特尔位涡(以下简写为PV_η)插值到z坐标系(由于位涡单位的量级远远小于风场和位温，因此，认为与插值法计算位涡时所引入误差相比，这一过程中引入的误差是次要的)，再以直接法结果(PV_η)为参考标准，分别在η坐标系和z坐标中比较插值法与直接法结果的差异(PV_η-PV_z)，从而说明直接法在减少计算误差方面的优越性。下面分别以0.89665、0.5256η面，及2600、4800、7000、9000 m高度上的结果为例，对插值法误差的水平分布进行分析。

3.2 插值法误差在η坐标系视角中的分析

图 3为0.89665和0.5256η面上的PV_η和PV_η-PV_z的分布；图 4所示为0.89665和0.5256η面上的PV_z和地形修正项C的分布。比较图 3a、b和图 4a、b发现，与直接法相比采用插值法并未对埃特尔位涡的水平分布形势产生明显的改变。在低层(0.89665η面上)，仅在(40°N，118°E)和(37°N，115°E)附近存在小块区域的负值位涡分布差别，而大部分地区的位涡正、负间隔形态，大致的梯度走向与大值区位置均无明显差别。在中层0.5256η面上，两张图更是高度一致，并未发现肉眼可见的区别。

然而通过PV_η-PV_z图发现，这两种方法所得位涡分布形势的高度一致性似乎掩盖了其在数值上的较大差异。图 3c、d中均存在成片的大于±0.5 PVU的差异区，甚至可达±2 PVU，与±20 PVU的位涡值域相比，这一量级的偏差已经无法被忽略。在低层(0.89665η面上)，有大片沿太行山脉呈南北走向分布的量级大于±1 PVU的偏差区域存在，显示出了在大气低层, 插值法误差对地形起伏具有较强的敏感性；而在中层(0.5256η面上)，在气旋中心附近高位涡梯度区分布着大于±0.5 PVU的偏差区域，偏差量级较0.89665η面要小，分布形态与位涡的分布形态类似，说明由低层到中层，插值法误差对地形的敏感性逐渐减弱，而对气象要素空间分布的敏感性增强，总体上表现出误差减少的趋势。

地形修正项C与地形的分布较为一致，在地表起伏较大的太行山脉和山东丘陵附近，从低层到中层均表现出不大于±0.1 PVU的修正值。值得注意的是，在忽略地形修正项C后PV_η和PV_z的表达式形式一致，但从结果来看PV_η-PV_z的量级远大于地形修正项C的量级，即若不考虑实际物理意义，在η坐标系中采用z坐标系同一公式计算位涡的结果，两者之间的偏差仍可达到(PV_η-PV_z-C)~(PV_η-PV_z)。从而说明了插值法在误差引入方面具有的可观效果，和直接法在精确计算方面的优越性。

在η坐标系中分析插值法误差，需要将插值法算得的PV_z再次插值到η坐标系，可能会进一步增加PV_z的误差，为了公平起见，下面将直接法算得的WRF-η坐标格点上的PV_η插值到z坐标系，并对其结果进行分析。

3.3 插值法误差在z坐标系视角中的分析

图 5、6分别为2600、4800、7000、9000 m高度上的PV_η、PV_η-PV_z的分布，代表了低层、中层、中高层和高层的情况。在低层2600 m高度上，偏差的分布与0.89665η面类似，主要偏差沿太行山脉呈南北走向分布，但范围和强度较0.89665η面有所减弱，可能说明将PV_η插值到z坐标系过程中引入的误差在一定程度上抵消了插值法计算PV_z的误差。中层4800 m偏差的分布与0.5256η面类似，均表现出受地形的影响减弱，主要偏差位于气旋中心附近高位涡梯度区，但偏差的量级高于0.5256η面达到了±1 PVU。中高层7000 m高度，PV_η-PV_z的量级较4800 m减弱，为±0.5 PVU，到达最小水平，且偏差的分布表现出随机性和平均化。至高层9000 m高度，偏差分布更加平均化且量级增加到±1 PVU水平。

通过以上分析发现，在z坐标系，PV_η-PV_z表现出与η坐标系中情况类似的量级分布，以及对于地形和高位涡梯度的敏感性，共同说明插值法在误差引入方面具有可观效果。下面利用本节所得数据，通过气旋中心附近偏差的剖面，对插值法引入误差的垂直结构作进一步分析。

3.4 插值法误差垂直结构

图 7、8分别为过(36.5°N，114°E)的PV_η、PV_z、PV_η-PV_z和(PV_η-PV_z)/PV_η百分数的纬向、经向垂直剖面。与3.3节相同，采用z坐标系视角，其中PV_η已插值到z坐标系。剖面选取如图 2白线所示，经纬向交叉点位于太行山脉以东、气旋中心以西地形坡度较大地区，便于分别考察地形和位涡梯度对误差分布的影响。

比较图 7a、b和图 8a、b发现，与η面图类似，采用插值法并未对埃特尔位涡的垂直分布产生明显的改变。纬向剖面中仅在(115.2°E，7 km)处存在PV_z大于PV_η的强正值中心，(114.0°E，10 km)附近存在小块区域的负值位涡分布差异；经向剖面由于取自大地形坡度地区，可以比较明显地看到8 km高度以上位涡的细微结构出现差别，但仍不足以对位涡分布的分析构成挑战。

进一步分析PV_η-PV_z的垂直分布。图 7c中的纬向剖面自西向东分别穿越了太行山脉、华北平原和山东丘陵地区。可以看出在大气低层，太行山脉和山东丘陵附近分布着量级大于±0.5 PVU的偏差而华北平原却鲜有分布。其中在整个太行山脉上空，广泛分布的低层偏差可从2 km高度延伸到4 km左右，表明直接法误差对大片分布的起伏地形更具有敏感性。到了4 km高度以上的大气中层，偏差对地形效应的敏感性减弱，大于±1 PVU量级的偏差开始集中于高位涡梯度区附近，表现出对气象要素空间分布较强的敏感性。从中高层8 km开始，偏差分布的随机性和平均化开始体现，10 km以上偏差开始趋于整层分布。图 7d中的经向剖面位于太行山东面的陡坡一侧，以受地形的影响为主，可以看出低层偏差广泛分布于4 km高度以下，4—7 km存在一个相对弱偏差层，8 km高度以上，大偏差开始趋于整层分布。对PV_η-PV_z的垂直分布分析的结果与前两节一致。

在对误差进行考察时，相对误差往往比绝对误差更能反映数据的可信度。从图 8c、d中的相对偏差百分数可以看出，中层附近对气象要素空间分布敏感存在较大绝对偏差的区域，其相对偏差百分数基本小于±10%；而在受地表起伏影响的中低层，表现为绝对偏差和相对偏差均很大，相对偏差百分数可以达到±50%。这说明和与气象要素空间分布有关的误差相比，由地形起伏而引入的误差，会对位涡计算的可信度产生更大的影响。特别是与统计气候相关的应用中，在地形起伏较大区域的边界层中，采用插值法计算的位涡数据，其较大的绝对、相对误差随时空累积后，会在很大程度上降低数据的可信度，为后续研究的开展制造麻烦，需尽可能避免。于是采用直接法在减少位涡计算误差中的优越性得到了进一步凸显。

3.5 插值法误差的统计分析

本部分通过对插值法均方根误差(相对于直接法结果)的分布进行分析，从统计学角度印证以上各节所得结果，对直接法在改善位涡计算误差方面的效果进行量化，并进一步讨论误差的时间稳定性。图 9(图 10)，为不同层次经向(纬向)平均的插值法均方根误差的纬向(经向)分布，图 11为不同海拔地区水平平均的插值法均方根误差的垂直分布，图 12为不同层次、海拔地区平均的均方根误差的时间变化。

从图 9可以看出，在400—3000 m的大气低层，自西向东随着地形平均高度的降低，插值法均方根误差逐渐减小，在高原处普遍大于0.7 PVU，而在平原处通常小于0.3 PVU。-90、0、90 km三处位置的均方根值类似，而在地形快速变化的-90 km处，插值法均方根误差远大于地形变化和缓的0 km和地形平坦的90 km。说明在低层，起伏的地形在插值法误差的产生中起了决定性作用。在-270、-90 km附近分布着两处位涡均方根峰值，在这些地方受到地形因素和气象要素分布的共同影响，插值法均方根误差也达到了1.5 PVU以上的峰值。而在5000—6000 m(中层)、8000—9000 m(高层)，插值法均方根误差受地形影响减弱，总体保持在0.4 PVU水平。

在图 10中，地表平均海拔高度的南北变化较为平缓，而在气旋中心附近气象要素变化剧烈，位涡均方根值快速增加。可以看出在400—3000 m高度，插值法均方根误差普遍大于0.4 PVU，但变化趋势与位涡均方根值不能很好对应；5000—6000 m处，均方根误差有所减小，总体保持大于0.3 PVU水平，但在气旋中心附近变化趋势与位涡均方根值保持一致，表现出了对气象要素分布的敏感性；至8000—9000 m(高层)，均方根误差分布更加平均化，且量级增大至0.4 PVU水平。插值法均方根误差的水平分布与前面各小结结果取得了一致。

太行山脉东侧坡度大，地表海拔由黄土高原的1000 m迅速降低至200 m左右，而河北中南部地区的地表海拔通常小于200 m(图 2)，因此以200 m作为分界，能够更好地表现出受地形影响地区和平原地区，插值法误差垂直分布的不同特点。图 11中平原地区，2000 m以下插值法均方根误差随高度升高迅速衰减，在2000—4000 m层次误差达到最小，4000—8000 m误差呈波浪形变化(在平原地区z坐标面与η坐标面贴合程度较高，当插值后z坐标面正好取在η面附近时，计算误差较小)，8000 m以上误差开始显著增长。而受地形影响较大的地区，8000 m以下均方根误差明显大于平原地区，且在4000 m以下差异较大，4000—10000 m误差略有增大，10000 m以上误差开始快速增大，但小于平原地区(平原地区高层网格垂直方向更加稀疏，误差更大)。

从均方根误差的时间变化(图 12)可以看出，模式启动初期经历了起转阶段(大约为6 h)，位涡均方根值和插值法均方根误差均呈现快速增长，19日08时之后，进入稳定变化阶段。插值法均方根误差与位涡均方根值的时间变化趋势能够较好地对应，插值法均方根误差随位涡均方根值变化的敏感性表现为：受地形影响大的地区(黑线)>平原地区(红线)，低层(实线)>高层(点划线)>中层(虚线)。黑色实线(地面海拔高于200 m地区的400—3000 m情况)对应的均方根值变化较为平缓，而均方根误差却随之产生了较大振幅变化；红色虚线(地面海拔低于200 m地区的5000—6000 m情况)对应的均方根值变化较为剧烈，而均方根误差的变化反而较为平缓，形成了鲜明对比。这可能是由于当低层的地形起伏处同时存在较强的位涡分布时误差增长较大(图 9)，而中层存在较强的位涡分布时误差增长较小(图 10)所致。同样可以看出，不同位置插值法均方根误差的大小差异表现为：受地形影响大的地区(黑线)>平原地区(红线), 低层显著而中高层不显著。垂直方向上均方根误差的分布整体表现出低层(蓝色实线)>高层(蓝色点划线)>中层(蓝色虚线)的形势。均与前文分析一致，并随时间表现出了一定的稳定性。

因此总的来看(图 9、10、11、12)，模式进入稳定运行阶段后，在4000 m以下的中低层，利用插值法计算位涡将引入0.5 PVU水平的误差，4000—8000 m高度该误差约为0.3 PVU，8000 m以上误差增大，至对流层顶达到1.2 PVU。

4 结论与讨论

首先推导了地形追随坐标系与z(p)坐标系的变换关系，并利用这一关系，得到与z坐标系形式(插值法)相差了地形修正项C的埃特尔位涡在地形追随坐标系中的表达形式(直接法)。同时，使用WRF模式对发生在华北地区复杂地形背景下的一次黄淮气旋爆发过程进行了模拟，利用模式η坐标系输出资料，采用以上两种方法分别计算了强降水时次埃特尔位涡，分别在η坐标系视角下、z坐标系视角下和垂直方向对插值法引入误差的分布进行了分析，并利用均方根误差对结果进行了统计，主要结论如下：

(1) 在大气低层，插值法误差对地形起伏具有较强的敏感性，并随高度逐渐减弱。而对气象要素空间分布敏感性较强的误差主要分布在大气中层。中高层以上，误差增大且分布趋于平均化。

(2) 对气象要素空间分布敏感的区域，存在较大的绝对误差；而对地形起伏敏感的区域，绝对误差和相对误差都比较大。

(3) 利用直接法计算位涡可有效减少插值法引入的计算误差(均方根)，中低层可达0.5 PVU水平，中高层可达0.3 PVU水平。

值得注意的是，本研究模拟结果采用的水平分辨率较高、垂直分层较细，而在更粗的网格空间中，例如在文中所选网格的高层区域，插值法所引入的误差将进一步增大，会对某些网格较粗的业务数值预报系统中的应用产生更大影响。此外，为了便于对地形敏感型误差和气象要素分布敏感型误差进行综合考察，本研究选取的研究范围为7经距×7纬距、水平500 km尺度的中尺度区域，在此尺度下插值法引入误差尚不足以影响到对位涡分布的分析。而对于几十千米或更小尺度的对流系统，选取更小的研究范围后误差将更容易在图上体现，此时插值法误差对位涡分布形势的影响将变得不容忽略。

在插值法中，z坐标系网格的选取也会对计算结果的误差产生影响，文中所采用的为等距200 m垂直间隔z坐标系网格。对不同垂直间隔的不等距z坐标系网格条件下，插值法误差的分布情况有待进一步研究。同时利用2.1节中推导的地形追随坐标系与z(p)坐标系之间的变换关系，可以进一步分析其他不同量级的z(p)坐标系定义物理量，在两种计算方法下的误差情况。

通过文中的分析与讨论，直接利用地形追随坐标系资料计算z坐标系下定义的埃特尔位涡方法的优越性得到了证明。应用在某些中小尺度对流性天气精细化结构的分析及气候统计研究中，可以有效减少插值法误差对结果产生的不良影响，提高计算和分析准确度。