引言

Zou等^[1]提出的BDA(Bogus Data Assimilation)方法使台风路径预报得到了改善，之后不少学者对BDA方案作了进一步研究并取得了较好结果^[2-5]。这种传统BDA方案构造的是一个轴对称台风，不能反映个别台风具体特征，有必要改用更为精细的非对称台风Bogus模型来做出改进。BDA方法一般只考虑体现动力性质的风场和气压场约束，但台风的温湿场非对称结构同等重要的影响着台风的演变趋势^[6]。BDA方法也不能调整台风周围大环境场，但台风与周围的环境场相互作用对台风的演变也有明显影响^[7]。而ATOVS卫星辐射亮温资料与温湿廓线有关，同时具有空间覆盖广的特点。国外一些研究及应用表明^[8-10]，利用三维或四维变分方法对NOAA极轨卫星TOVS/ATOVS特别是AMSU微波辐射率资料进行直接同化，对数值预报水平的提高有较大贡献, 对模式预报水平的影响已经达到或超过传统的探空观测资料。毕宝贵等^[11]对TOVS和云导风反演资料进行变分同化得到合成分析场也表明比仅用常规探空资料更接近实况，弥补了实况资料的不足。万齐林^[12]用Grapes三维变分同化系统对校准后的ATOVS反演产品进行同化分析也表明了ATOVS反演资料的价值。潘宁等^[13-14]对ATOVS资料直接变分同化技术做了详细介绍，并用MM5-3DVAR系统进行同化实验得到一些极具参考意义的结论。张华等^[15-16]对西北太平洋台风威马逊作AMSU微波资料3DVAR直接同化试验，对台风环流及湿热场产生了正面影响。任强等^[17]则更为精细地考虑了台风数值预报中AMSU微波资料同化时云雨区的资料处理问题。从这些工作所启示，对ATOVS资料进行同化正好能弥补BDA方法中温湿场及台风周围大环境场调整上的不足。另一方面，台风降水的预报也是一大难点，一些工作表明了台风降水的复杂性^{[18-20, 33-34]}, 而对与大气湿度廓线密切相关的ATOVS资料进行同化正好可以为台风降水预报提供必要参考。

文中针对登陆台风韦帕，在BDA方案中同化一种非对称海面风压场Bogus资料的同时，引进快速辐射传输模式RTTOV8，利用四维变分同化技术，把多星多轨道ATOVS卫星辐射亮温资料同化到模式中去，优化初始场台风湿热结构及台风外围环境场，以获取更精细的初始台风环流和温压湿场结构，提高台风预报效果。考虑到韦帕台风陆面移行特性，试验中引入了Noah陆面过程方案。

1 使用资料介绍 1.1 非对称Bogus资料

在以台风中心为原点的极坐标中，考虑海面摩擦作用，假定成熟台风呈稳定状态，并在惯性项中考虑台风移动对曲率半径的影响，从台风域内海面任意空气质点的水平运动方程可导出切向风和径向风的表达式：

$ v_\theta ^3 + {a_1}v_\theta ^2 + {a_2}{v_\theta } + {a_3} = 0 $

(1)

$ {v_r} = (-d\partial p/\partial \theta + rk{v_\theta }\cos \varphi )/({v_\theta } + {v_s}\sin \alpha + fr + rk\sin \varphi ) $

(2)

式中k为摩擦系数，${a_1} = 2({v_s}\sin \alpha + fr + rk\sin \varphi ), {a_2} = {({a_1}/2)^2} + {r^2}{k^2}{\cos ^2}\varphi-rd\partial p/\partial r $, ${a_3} =-f{r^2}d\partial p/\partial r-({v_s}rd\partial p/\partial r)\sin \alpha-(k{r^2}d\partial p/\partial r)\sin \varphi - (rkd\partial p/\partial \theta )co{\mathop{\rm s}\nolimits} \varphi $。φ为摩擦阻力偏离实际风矢量反方向的夹角，d为比容。式(1) 求解可采用邢富冲^[21]的卡当方程求解方法。v_θ, v_r的求算关键是∂p/∂r和∂p/∂θ的求取。引入藤田公式^[22]，考虑移动台风的非对称性质，并进行风廓线约束推广，得出如下形式：

$ p(r, \theta ) = {p_m}(\theta )- \Delta p(\theta )/{(1 + a{[r/{r_0}(\theta )]^b})^{1/2}} $

(3)

其中p_m为台风外围环境气压，p_c为观测的台风中心海平面气压, Δp=p_m-p_c，a为表征最大风速半径位置的参数，其与方向有关。b为最大梯度风控制参数。p_m为环境气压，其受到赤道高压或副热带高压等大尺度系统的影响时具有明显的方向性，若台风只受副高引导沿高压边沿移动，此时靠近副高一侧p_m比远离副高一侧大，因此对p_m估算如下：

$ {p_m}(\theta ) = {p_\infty }-{p_a}\sin (\theta-{\alpha _0}-{\beta _1}) $

(4)

p_∞为平均环境气压，取其经典值1010 hPa^[23]。α₀为台风移向，β₁=arcsin(δsinβ₀)，β₀为实际移向与内力方向的夹角，逆时针方向为正。内力方向取为西北偏北向^[24]，取δ=1/8^[25]。而后采用Kuihara等^[26]的方法确定出分析台风中心位置及半径，此半径上的最大气压和最小气压之差即为p_a，再利用改进的MM5滤除方案^[27]，得出消去分析台风的基本风场。

为了得到台风常数r₀，选取原分析场台风外围的一条既能反映台风系统又能反映周围天气系统影响的闭合等压线L(θ)，进行傅里叶级数插值拟合后并代入式(3) 有

$ {r_0}(\theta ) = L(\theta ){\{ [\left( {\frac{{\Delta p(\theta )}}{{{p_m}(\theta ) - {p_L}}}} \right)^2 - 1]/a\} ^{ - 1/b}} $

(5)

p_L为L(θ)对应的气压。结合式(3) 即可得出∂p/∂r和∂p/∂θ，代入式(1) 和(2) 即可求出台风域内风场，然后将此风场与扣除移速后的基本场叠加得到最终风场，这样，基本场风速中的非对称成分得以部分地保留下来。气压场则由式(3) 提供。文中不考虑海面非线性耗散作用，取a=1.2，k=1.7×10^-4[23]，φ=38°^[28], p_c=948 hPa, p_a=3.5 hPa，调整参数使台风域内最大风速等于观测得到的近地面最大风速即可。

1.2 ATOVS卫星辐射亮温资料处理

卫星观测资料为ATOVS-L1C资料。同化时前需做一些预处理，包括AMSU-B观测资料稀疏、临边检查和海岸检查、极值检查等^[20]。同化通道HIRS为5~7、11~16通道，AMSUA为4~8通道、AMSUB为3~4通道。HIRS资料剔除云量大于7成的廓线^[29]，AMSUA资料的云检测使用降水概率法^[30]，AMSUB资料的云检测采用Bennartz散射因子法^[31]。

2 数值试验的设计 2.1 同化系统的设计

在MM5 4D-Var系统的基础上，引入了快速辐射传输模式RTTOV8，使之能够直接同化卫星辐射亮温探测资料。目标泛函可以定义为：

$ J = {J_B} + {J_{{\rm{Bogus}}}} + {J_{BT}} $

(6)

其中J_B表示模式控制变量X与背景场X_B的偏差，J_Bogus表示模式模拟量与Bogus资料的偏差，J_BT表示模式模拟的卫星辐射亮温值与实际观测的卫星辐射值之间的偏差。J_BT、J_Bogus具体表达式可参见文献[1]。J_B具体表达式可参见文献[13-14], 其中观测算子即应用RTTOV-8正演计算亮温资料的过程，其主体是RTTOV-8模式，但还包括了一些水平插值算子、垂直插值算子、通道选择算子和质量控制算子等。文中ATOVS资料的各通道权重系数给定的原则是使各个通道对代价函数的贡献大致均等。

同化卫星资料时，首先给定初始迭代值y₀⁰，在预先设定的同化时间窗[t₀，t_N]内正向积分MM5预报模式，每6个时步存储一次预报模式的状态矢量。同时在有卫星观测的时刻，在卫星观测点正演模拟ATOVS参与变分同化通道的辐射值，计算目标函数中背景场项与卫星资料强迫项的值；当预报结束，从同化窗的截止时刻开始，在同化时间窗[t₀，t_N]内反向积分MM5伴随模式，这里伴随模式的初值取为：${\delta ^*}y({t_N}) = {W_N}(y_N^0-{\hat y_N}), {\hat y_N}$为观测量。在此积分过程中，伴随模式的基本状态量由MM5预报模式正向积分时存储的状态矢量y_N⁰提供。而卫星强迫项的梯度则在具有卫星观测的时刻，带入相应的RTTOV-8正演模拟值与NOAA卫星实际观测的偏差$ {W_N}(y_N^0-{\hat y_N})$反向积分RTTOV-8的伴随模式得到。当反向积分到t₀时刻得到的值δ^*y(t₀)即为目标函数相对于初始控制变量的梯度∇J(y₀)。随后，按照共轭梯度法确定的目标函数的下降方向，用单向搜索法确定目标函数沿共轭方向充分下降的最优步长，同时利用有限记忆拟Newton方法确定新的初始值y₀¹。重复以上步骤，就可以得到一系列的y₀^k，直到按照预先给定的精度逼近最优的(y₀)_p。最后以(y₀)_p作为MM5预报模式的初始场进行预报，检验卫星资料对模式预报精度的改进。

RTTOV8的43个等压面的温度、水汽以及地表温度、地面附近风速和气压、云顶气压和云量等要素，主要由MM5模式提供。对于温湿廓线中MM5模式层顶以上部分，利用NESDIS标准气候廓线库合成不同纬度带、不同季节的气候参考廓线来填补^[29]。臭氧取自模式本身提供的参考廓线。

2.2 试验方案的设计

选取0713号台风韦帕个例，以2007年9月17日18时(UTC)作为初始时刻进行四维变分同化试验并作48小时模拟预报。模拟采用的数值模式为非静力版本的MM5模式及其伴随模式和RTTOV8模式及其伴随模式。分析资料为6小时间隔的1°×1°的NCEP分析资料、台风实测资料来自中国国家气象中心的分析实况。试验区域以25.5°N, 124.5°E为中心，水平网格为75 × 91，网格距为54 km，垂直方向为分布不均匀的23层。同化窗口为172分钟，同化迭代次数为25，优化算法采用有限内存准牛顿算法^[32]，积云参数化为GRELL方案，边界层参数化为MRF方案。4种数值试验方案如表 1所示。

ATOVS资料轨道只覆盖约十几个纬度的区域，而同一卫星经过同一地点的时间间隔较长不宜同时同化。为了扩大可用资料量，采用多星多轨道的同化方案，在0时刻加入17日17时44分的NOAA-18卫星ATOVS资料，73分钟和171分钟分别加入19时13分和20时51分的NOAA-16卫星的ATOVS资料，ATOVS卫星资料近乎覆盖了整个模式计算区域。另外，重点考察同化ATOVS资料的效果，下文只涉及方案1，3，4的初始场和预报场试验结果分析。

3 数值试验结果分析

图 1给出了初始时刻沿台风中心的温度场和经向水平速度场纬向剖面图。从图中看到，方案1的温度线几乎东西向平直，而方案3和方案4皆出现了暖心结构，其中方案4因为加入卫星资料，等温线出现了东西非对称性及一些中尺度扰动，由于台风中心附近的ATOVS资料因云检测原因被剔除掉了，故卫星资料所携带的信息是从外围向中心传播的。从水平风场分布来看，后两种方案最大风速比方案1高，东西两侧数值差异增大，出现更明显的非对称性。而方案4的风场与方案3相比在高低层都变弱了，原因为ATOVS辐射亮温资料与模式变量的高度非线性制约了目标泛函的总体下降程度使风场的调整达不到单独同化Bogus资料时的程度, 而非ATOVS资料本身的影响作用。同样可以看出卫星资料是从外围影响台风环流结构的。

图 2给出了方案3、4初始时刻600 hPa等压面上的云雨水含量。加入了大量ATOVS卫星资料的方案4云雨水含量远大于方案3。从图中也看出在卫星资料覆盖了整个模式区域情况下，优化得出的云雨水也只出现在台风周围的局部区域，在南边尤其台风西南侧出现了大片的水汽输送带。这说明多星多轨道同化方式能较好地优化出云雨水含量的合理分布。再与2007年9月18日02时(北京时)的风云二号气象卫星云图(图略)相比较来看，云雨水分布与云图分布具有一定的对应性。从云雨水分布图中还看到，选取多星多轨道的卫星资料进行同化，对整个大环境场结构都有所调整，也反映出了大尺度云系结构及中小尺度云雨水分布及水汽输送状况。因Bogus资料的抑制和云检测等原因，方案4中台风中心附近区域上空的云雨水含量没有得到调整。

图 3给出了试验模拟12小时得到的经过台风中心的位涡纬向垂直剖面。3个方案模拟12小时后在台风中心附近都维持一个涡柱结构。而方案1中位涡结构散乱，不能较好体现台风结构；方案3则是分别在中层和低层存在两个高中心，低层中心数值达到9，且高于上层；方案4由于大量卫星资料的引入，其结构与方案3有所不同，其结构细腻，体现出更多的中尺度信息，且在高层、中层和低层都分别出现了高值中心，而高层中心数值大于低层。位涡作为一个综合反映大气动力学和热力学性质的物理量，可定性地反映出台风的发展趋势^[35-37]。当对流层中高层存在一个位涡高值中心时将会有利于对流层低层气旋性涡度的发展即有利于台风的发展。因此方案3不利于台风的发展。实际情况中台风仍处加强时期，因而方案4的位涡分布相对合理。然而方案4中的位涡峰值小于方案3，其原因与图 1风场的差异相对应。另外，不同于初始时刻ATOVS资料只影响台风外围的情况，模拟12小时后卫星资料信息已经从台风外围传播到了台风内部。

图 4是4种方案间隔6小时的路径预报误差及强度演变误差。几种同化方案的路径预报效果都优于控制试验。前12小时内方案2、3、4的路径基本一致，之后经台湾岛时产生分离，24小时前后登陆，之后3个同化方案预报路径都继续偏离实况，没有深入内陆，而是在海岸附近擦边而过。然而方案3因为考虑了Noah陆面过程方案，其误差小于方案2。方案4加入ATOVS资料后台风的路径预报改善不明显，说明表征动力学变量的Bogus资料一定程度上控制着路径预报模拟效果。强度误差图中，由NCEP分析资料得到的强度预报效果较差；方案2的结果有了相当大的改善，整个过程预报误差较为平稳；考虑Noah陆面过程之后，方案3的预报强度稍比方案2弱，这恰好与实况更为接近。方案4加入大量ATOVS卫星资料后，基本能保持BDA方法的强度控制优势。同化方案试验中，初始时刻台风处于成熟期，于是模拟试验都存在预报中后期台风过度发展的现象，这是台风后期路径预报不准确的原因之一。

图 5是方案1、3、4模拟得到的24小时降水量(大于3 cm)。方案1降水分布零散且雨量偏小；方案3的降水强度比方案1增大，降水落区与Bogus资料相对应；方案4的降水范围比方案3扩大，降水落区发生改变，降水强度进一步加大。与2007年9月19日02时(北京时)的风云二号气象卫星云图(图 5d)比较得知，方案4的雨带分布较好地符合了卫星云图的云带分布，且降水落区与方案1较为相似，原因为NCEP分析资料事实上已经过众多观测资料的大尺度同化、分析及订正，Bogus资料的加入使这种订正失效，而在中尺度区域模式中同化ATOVS又使得大尺度背景场的订正得以部分恢复，同时降水强度较方案3稍有增加，降水落区发生改变。由于海上降水观测资料的缺乏及篇幅有限，此处不再探讨ATOVS资料同化对降水状况的影响。

5 结论

(1) 提出的非对称Bogus方法包括：从考虑摩擦的水平运动方程导出卡当方程形式的风速公式并求解；对台风外围闭合等压线进行傅里叶级数插值拟合；计算副高影响下各个方向的环境气压；对藤田公式作非对称推广和风廓线约束；Bogus风场吸收滤除弱涡旋并扣除台风移速后的基本场成分。该方法全程以非对称形式计算，避免了假定台风轴对称时引入非对称量的难题，也避免了假定台风呈椭圆时所带来的难以确定风场长短轴方向的难题，还避免了给定最大风速半径时的不确定性。该方法考虑非对称性质并吸收基本场信息及估算副高影响，只针对台风初值化设计，其他场合未必适用。

(2) BDA方法对台风韦帕临近登陆及登陆后的路径预报误差较大，而考虑Noah陆面过程后路径预报效果有很大提高。

(3) 同化ATOVS资料对台风韦帕动力场影响较小，能保持BDA的强度控制特性和路径预报效果；同化ATOVS资料能有效调整台风韦帕初始场湿热场结构，并加入中尺度结构信息，对台风韦帕的模拟降水强度及降水落区产生较大影响。

(4) 成熟台风在模拟后期往往过度发展，这是BDA方法延长预报时效需要解决的问题之一。

(5) 怎样有效地把ATOVS资料尤其HIRS红外资料应用到台风数值预报中的云雨区中，这仍是一个需要解决的难题。