确定大气边界层顶高度的新方法及数值实验*

周建印 项杰 黄思训2)†

1) (国防科技大学气象海洋学院, 南京 211101)

2) (国家海洋局第二海洋研究所, 卫星海洋环境动力学国家重点实验室, 杭州 310012)

提出了一种确定大气边界层顶高度的数值微分新方法, 该方法使用了正则化技术, 把对弯角廓线求导数的数值微分问题转化为求目标泛函极小值的问题, 采用双参数模型函数方法来选择正则化参数, 最后利用最大梯度法确定边界层顶高度.首先通过两个数值实验验证了新方法的有效性, 实验结果显示, 随着掩星资料噪音的增多, 由差分法和结合L曲线方案的数值微分方法得到的边界层顶高度波动增大, 而通过双参数模型函数方法得到的高度很稳定, 这说明新方法能够很好地过滤噪音, 从而保留廓线的主要信息.随后基于2007—2011年1, 4, 7, 10月的COSMIC弯角数据, 利用新方法分析了全球海洋大气边界层顶高度的季节特征, 并与用掩星资料自带的大气边界层顶高度数据zbalmax得到的季节分布进行对比.结果表明, 两者的季节分布特征十分一致: 海温相对周围海域高的区域, 边界层顶高度较高, 反之, 边界层顶高度较低; 在暖流经过的海域, 边界层顶高度较高, 在寒流经过的海域, 边界层顶的高度相对较低.

1 引 言

行星边界层是对流层底层的大气, 下垫面主要通过热量辐射传输、摩擦阻力、水汽通量、污染物排放以及地形对边界层发生作用[1].边界层的结构复杂多变, 边界层顶高度随时空变化而变化, 通常定义为混合层(常值通量层)与自由大气过渡的高度[2], 高度变化范围在几百米到几千米之内.边界层顶高度在数值预报模式中的边界层参数化方案以及气溶胶的反演中有着重要的作用, 具体表现为: 如果边界层顶高度太低, 则边界层与云层解耦,从而抑制了热量、水分和湍流动能从海洋表面到云层的垂直传输, 并可能加速云的消散; 如果边界层顶高度太高, 则会形成积云而不是层状云[3], 所以确定边界层顶高度是一项十分重要的工作.

边界层顶往往存在湿度、温度、折射率、气溶胶粒子以及云滴的急剧变化[4], 无线电掩星探测获得的弯角和折射率廓线在边界层顶处往往具有大的垂直梯度, 因此, 最大梯度所在的高度被认为是边界层顶高度[5].在计算边界层顶高度的观测数据方面, 由于边界层顶高度的变化在几百米左右, 而被动红外和微波探测的分辨率在1—2 km, 故使用分辨率大约在20 m的掩星资料计算边界层顶高度更加精确, 而且全球定位系统掩星探测(global positioning system radio occultation, GPS-RO)技术具有全天候、全球覆盖、长期稳定等优点, 使得在分析全球边界层顶高度时更具有优势.

前人在计算边界层顶高度上, 提出了许多方法.Holzworth[6]通过从最高表面温度沿干绝热线随高度上升到与最近观测到的温度廓线的交点, 来得到最大混合层高度; Coulter[7]比较了由激光雷达、声雷达以及温度廓线得到的混合层顶高度, 发现激光雷达与声雷达探测的高度基本一致.Van Pul等[8]由无线电探空资料通过气块法(parcel method)得到了边界层顶高度.在 20 世纪, Bianco和Wilczak[9]基于模糊逻辑方法的混合深度识别算法来确定边界层深度; Lokoshchenko[10]提出了一种方法, 即把参数积分(如不稳定能量)最大值所在的高度作为混合层高度; Balsley等[11]使用最小风切变和最大Brunt-Väisälä频率对夜晚稳定边界层顶高度进行计算; Sokolovskiy等[12]把COSMIC掩星探测折射率梯度的最大下降点的高度确定为边界层顶高度, 2007年Sokolovskiy等[13]又利用COSMIC掩星弯角资料对边界层顶高度进行了计算; Baars等[14]修正了小波协方差方法, 利用雷达资料得出边界层顶高度; Seidel等[15]用无线电掩星资料比较了计算边界层顶高度的6个气象要素(如温度、湿度、折射率等), 结果显示折射率和位温的梯度最小值得出的边界层顶高度较高;Dai等[16]使用多个外场观测项目中飞机测量的廓线资料, 重新计算了边界层顶高度; Yan等[17]把带有正则化技术的数值微分方法应用于COSMIC掩星弯角资料来确定大气边界层顶高度.

正则化技术是苏联科学院院士吉洪诺夫(Tikhonov A N)[18]为了获得线性不适定问题稳定的近似解而于1963年提出的.随着研究的深入,正则化的方法得到了蓬勃的发展, 目前常见的正则化方法主要有Tikhonov正则化方法[19], Lavrentiev正则化方法[20], 迭代正则化方法(如Tikhonov迭代方法[21]、Landweber迭代方法[22]), 以及基于谱理论的正则化方法(如TSVD (谱截断)正则化方法[23])等.前人在大气科学领域中应用正则化方法做了很多有意义的工作[24−28], 而其中正则化方法的一个关键环节是正则化参数的选择, 常用的正则化参数选择方法有L曲线准则[29−31]、广义交叉验证准则[32,33]、Morozov 偏差准则[34]等.其中,Morozov偏差准则是一种应用非常广泛的后验选取策略, 但是此原则仍然不够理想, 于是1984年Morozov[35]提出了吸收的Morozov准则(damped Morozov principle), 它是通过求解吸收的 Morozov偏差方程来得到正则化参数, 但该方法在进行数值求解时, 收敛速度慢且计算量巨大.于是在1998年, Kunisch和Zou[36]提出单参数模型函数法来求正则化参数, 从此模型函数法为减少计算量指出了一条明确的道路.

在Yan 等[17]提出的带有正则化技术的数值微分中, 正则化参数的选择采用了L曲线方法(因此, 该方法在下文中简称为L曲线方法), 但是此方法计算量大, 且L曲线方法的收敛性问题在理论上并没有得到解决.因此, 本文提出一种计算大气边界层顶高度的新的数值微分方法, 其中仍然应用了正则化技术, 但是正则化参数的选择采用了双参数模型函数方法.本文第2节首先介绍了研究所用到的数据, 然后介绍双参数模型函数方法的原理;第3节进行数值实验, 并与另外两种方法的结果作比较, 以验证方法的合理性; 第4节应用双参数模型函数方法得到全球海洋边界层顶高度的季节分布特征, 并与利用掩星资料自带的边界层顶高度得到的结果作对比; 第5节是总结.

2 数据与方法

2.1 数 据

FORMOSAT-3/COSMIC是美国与中国台湾合作的一个GPS掩星探测项目, 本文使用的弯角廓线数据来自COSMIC项目, 包括2007—2011年的 1, 4, 7, 10月份的全球的弯角数据 (https://cdaac-www.cosmic.ucar.edu/cdaac/products.html),总样本为1074349条廓线, 其中平均每个月的数量在50000—70000条, 4和7月份比1和10月份数量多 30000—40000 条 (图1); 同时, atmPrf文件本身还自带了利用弯角数据得到的大气边界层顶高度数据zbalmax.考虑到数据不可避免会带有缺省缺测等不足, 所以从以下两方面进行质量控制.1)筛选数据.考虑到边界层顶基本都在500 m以上, 5 km 以下, 所以本文采用这样的弯角数据,即廓线的最低高度在0.5 km以上, 最高高度在5 km 以下.2)控制数据长度.选择有效数据的长度至少在5个以上, 以减小由于数据缺少而带来的误差.

图1 2007−2011年1, 4, 7, 10月份掩星廓线数量的时间分布Fig.1.Temporal distribution of the number of occultation profiles in January, April, July and October, 2007−2011.

2.2 双参数模型函数方法

大气边界层顶常常会有温度或水汽压的急剧变化(温度急剧上升或水汽含量急剧下降), 因此,可以采用最大梯度法来确定边界层顶高度.近些年来, 使用“开环”跟踪技术之后, GPS 掩星资料在低层大气的准确度有了显著提高, 使得GPS掩星资料可以提供更多有价值的低对流层(尤其是离地面2 km以下)大气的信息, 但是由于水汽的影响,低对流层(尤其是离地面2 km以下)大气廓线数据仍然含有负偏差.Yan 等[17]提出的数值微分方法中, 正则化参数的选择采用了L曲线法, 由于L曲线法具有前面提到的两个缺陷, 所以本文采用双参数模型函数方法来确定弯角梯度.

2.2.1 模型的构建

模型构建的基本思路与Yan 等[17]相同, 为了文中叙述的方便, 下面对模型的构建作简要介绍.设弯角为连续可微的函数 α (z) , 弯角廓线的梯度为 φ (z)= α′(z) , z 为高度, 范围为 z0⩽z⩽zn.

在高度范围 I =[z0,zn]中取i=1,2,···,n, 根据牛顿莱布尼兹公式可以得到

对于(1)式右侧, 可以用辛普森法则近似为

其中, 2h=zi+1−zi−1=2(zn−z0)/n, 本文将廓线等分为 n 份, 每一份高度 h 为 10 m, 误差为error= −h5/90φ(4)(ξ),ξ∈ [zi−1,zi+1], 如 果 弯 角梯度变化剧烈, 比如存在锋面过程等剧烈变化的天气时, 则误差不容忽视.

将(2)式代入(1)式可以得到

经整理后可得到弯角梯度应满足的线性方程:

其中

此时, 问题归结为求线性方程(4)的解.

2.2.2 双参数模型函数方法的应用

考虑模型(4), 如果观测资料b带有高频成分,则在反演大气弯角垂直梯度时, 结果会产生较大的误差.所以, 为了克服 (4) 式的不适定性, 本文把(4)式转化为求解如下目标泛函极小值的问题:

其中X是满足泛函极小值的解, α ,β 是待求的双正则化参数, 而矩阵L为一阶导数算子, 表示对x进行约束, 使得x的梯度变化不至于过于剧烈,起到平滑的效果.

本文与Yan 等[17]文章中的不同之处是采用双参数的模型函数方法[36−39]来确定目标泛函中两个参数的最优解, 进而求解出弯角梯度, 其基本技术路线如下:

令目标函数 F (α,β) 为

此时吸收的Morozov偏差方程为

式中 γ >1,µ>1 为吸收系数, δ 为误差水平.因为(8)式是关于 α ,β 的高度非线性方程, 用通常的迭代方法求解结果并不理想(比如Newton法和拟Newton法), 原因是他们只具有局部收敛的性质,而且对初值的要求比较高.因此, 我们采用近年来得到广泛应用的模型函数法来确定正则化参数, 模型函数法的好处是不但计算量大大减少, 而且还保证了收敛性.模型函数有各种选择的方法[38,39], 为简单起见采用线性模型函数mk(α,β):=Tk+αCk+βDk表示第k次迭代的结果以逼近 F (α,β) , 其中

图2 双参数模型函数法流程图Fig.2.Flow chart of two-parameter model function method.

将(9)式代入(8)式可得

确定正则化双参数的计算流程如图2所示.因为在正则化方法中物理量往往是无量纲的, 为了使正则化中的量纲一致, 将 α ,β 做了以下处理: 将 α 扩大500×std(b)倍, 将 β 扩大 s td(b) 倍, 这里 std(·) 表示取标准差.

3 数值实验

本节比较了三种确定边界层顶高度的方法, 分别是模型函数法、差分法和L曲线方法.所用数据是 2011年1, 4, 7, 10月份的 COSMIC 弯角资料.在处理之前先对数据进行质量控制, 然后借助spline插值法将资料插值到间隔为10 m的网格点上.双参数模型函数法的两个初始参数 α ,β 都取为 1, 根据弯角的量级为 10−2, 故取 δ =10−2, 收敛参数 ε =10−4.由差分法、L曲线法以及模型函数法计算出来的边界层顶高度分别记为 Hd,HL,HM.

掩星探测资料不可避免地会带有误差, 为了比较在存在噪声情况下几种方法的表现, 本节将[−δ,+δ]上的随机误差添加到弯角资料, 由于在边界层中弯角的量级为 10−2, 故令 δ ∈[0.0001,0.01],在带有噪声的弯角资料的基础上, 通过与差分法和L曲线法相比较, 检验模型函数法在确定边界层顶高度上是否具有准确性.

图3 由差分法和模型函数法得到的弯角梯度廓线(BA表示弯角) (a)弯角廓线添加随机误差 δ =0.0025 , 边界层顶高度HM=1.2km,Hd=1.4km ; (b)弯角廓线添加随机误差 δ =0.005 , 边界层顶高度 HM=1.1km,Hd=1.4km ; (c)弯角廓线添加随机误 差 δ =0.0075 , 边界层 顶高 度 HM=1.2km,Hd=1.8km ; (d)弯角 廓线 添加随 机误差 δ =0.01 , 边界 层顶高度 HM =1.2 km, Hd = 3.7 kmFig.3.Angle gradient profile obtained by the difference method and the model function method using the bending angle gradient profile (BA represents the bending angle): (a) Bending angle profile with uniform random error δ =0.0025 , boundary layer top height HM=1.2km,Hd=1.4km ; (b) bending angle profile with uniform random error δ =0.005 , boundary layer top height HM=1.1km,Hd=1.4km ; (c) bending angle profile with uniform random error δ =0.0075 , boundary layer top height HM=1.2km,Hd=1.8km ; (d) bending angle profile with uniform random error δ =0.01 , boundary layer top height HM=1.2km,Hd=3.7km.

图3选取了廓线atmPrf_C001.2011.182.05.19.G18_2013.3520_nc, 并且分别添加随机误差δ=0.0025,0.005,0.0075,0.01, 比较了差分法和数值微分结合模型函数方法得到的边界层顶高度.从结果可以看出, 图3(a)是由带有误差 δ =0.0025 的弯角得到的弯角梯度, 差分法得到的弯角梯度廓线在一定水平范围内呈密集的锯齿状, 其最小梯度在1.4 km左右, 数值微分方法结合模型函数法得到的弯角廓线比较光滑, 其最小梯度在1.2 km左右, 由图3(b)—(d)可知, 随着弯角资料误差的增大, HM始终稳定在 1.2 km 左右, 而 Hd由 1.4 km(图3(a)和图3(b))变为 3.7 km(图3(d)), 变化幅度超过2 km, 说明当掩星资料中存在观测误差时,模型函数法能够很好地过滤噪声, 保留主要信息,得到的边界层顶高度的结果更加稳定.

图4 由差分法、L曲线法和模型函数法得到的边界层顶高度 (a) 三种方法基于廓线1得到的边界层顶高度, Htrue =3.15 km, std(HM) = 0.013, std(HL) = 0.44, std(HM) = 0.61;(b) 三种方法基于廓线2 得到的边界层顶高度, Htrue = 4.55 km,std(HM) = 0.020, std(HL) = 0.89, std(HM) = 1.19Fig.4.Height of the boundary layer obtained by the three methods: (a) Three methods to get the height of the boundary layer top based on profile 1, Htrue = 3.15 km, std(HM) =0.013, std(HL) = 0.44, std(HM) = 0.61; (b) three methods to get the height of the boundary layer top based on profile 2,Htrue = 4.55 km, std(HM) = 0.020, std(HL) = 0.89, std(HM) =1.19.

图4给出了三种方法(模型函数法、差分法、L曲线法)分别在两条添加误差后的弯角廓线(atmprf_C001.2001.182.00.22.G23_2013.3520_n c, atmPrf_C001.2011.182.00.39.G20_2013.3520_nc)上的表现, 以结果的标准差衡量方法的稳定度.由图4(a)可以看出, HM比较稳定, 保持在 3.15 km左右, 而 Hd和 HL随着误差的增大, 波动越来越大,其中, HM,Hd,HL的稳定度分别为 0.013, 0.44 和0.61; 同样由图4(b) 也可以看出, HM依然是比较稳定, 保持在 4.55 km 左右, 而 Hd和 HL随着误差的增大, 波动也越来越大, 其中 HM,Hd,HL的稳定度分别为0.02, 0.89和1.19.由以上分析可以看出,弯角资料的噪声越大, Hd振荡程度就越强, HL振荡的程度次之, 而 HM的稳定性则比较好, 说明对于含有噪声的弯角资料, 模型函数法也能稳定地得到边界层顶高度, 所以模型函数法在确定边界层顶高度时, 准确性能得到保证.

4 全球边界层顶高度的季节特征

本节应用模型函数法, 对海洋上边界层顶高度的季节特征进行分析.图5展示了2007—2011年的 1, 4, 7, 10 四个月份年平均的边界层顶高度, 可以看出, 1月份南半球的边界层顶达到最高, 尤其是中高纬度海洋上的边界层顶高度达到2 km以上;4月份南半球的边界层顶高度开始降低, 南美洲与南极大陆南设得兰群岛之间的德雷克海峡附近降低的程度最为明显, 而北极附近海域的边界层顶高度也降低, 原因在于北极的海冰覆盖面积在3月份达到高峰, 4月份海冰的面积比1月份多[40], 所以4月份下垫面温度更低, 相应的边界层顶高度也在4月份比1月份低; 七月份南半球中高纬度的边界层顶高度在四个月份中最低, 北半球的低纬度地区边界层顶高度明显增高, 如墨西哥湾地区、阿拉伯地区以及日本海地区, 而受西太平洋副热带高压下沉气流影响, 我国东南沿海、孟加拉湾和菲律宾附近海域边界层顶高度明显降低; 10月份南半球中高纬边界层顶高度略有回升, 北极海域边界层顶高度开始降低, 而此时我国东南沿海、孟加拉湾和菲律宾附近海域边界层顶高度明显升高.

从以上现象分析可以得知, 随着太阳直射点的南北移动, 海洋上的边界层顶高度呈现明显的季节变化特征, 边界层顶高度与海面的温度密切相关,当海面的温度相对周围平均温度越低时, 边界层顶高度越低, 反之, 边界层顶高度就越高.但是四个月份共同的特征是, 在寒流流经的区域, 边界层顶高度普遍偏低, 比如巴西寒流和加利福尼亚寒流,暖流流经的区域边界层顶普遍较高, 比如赤道暖流.图6展示了利用掩星自带的zbalmax的资料求得的边界层顶高度的季节分布, 通过与图5相比较, 可以看出模型函数法得到的边界层顶高度略高, 但两者边界层顶高度的时空分布十分一致, 说明模型函数法得到的季节特征具有准确性.

图5 用模型函数法得到的海洋5年平均的边界层顶高度(所用资料为2007—2011年1, 4, 7, 10四个月份的掩星弯角的资料)(a) 1月份平均边界层顶高度; (b) 4月份平均边界层顶高度; (c) 7月份平均边界层顶高度; (d) 10月份平均边界层顶高度Fig.5.The 5-year average boundary layer height of the ocean obtained by the model function method, the data used is the bending angle profile of the four months of 2007−2011 in January, April, July, and October: (a) The average height of the boundary layer in January; (b) the average height of the boundary layer in April; (c) the average height of the boundary layer in July; (d) the average height of the boundary layer in October.

图6 海洋 5年平均的边界层顶高度 (所用资料为 2007—2011年1, 4, 7, 10 四个月份的掩星自带 zbalmax 的资料) (a) 1月份平均边界层顶高度; (b) 4月份平均边界层顶高度; (c) 7月份平均边界层顶高度; (d) 10月份平均边界层顶高度Fig.6.The 5-year average boundary layer height of the ocean obtained by the model function method, the data used is the zbalmax provided by CDAAC of the four months of 2007−2011 in January, April, July, and October: (a) The average height of the boundary layer in January; (b) the average height of the boundary layer in April; (c) the average height of the boundary layer in July; (d) the average height of the boundary layer in October.

5 总 结

本文提出了一种确定边界层顶高度的新的数值微分方法, 通过模型函数法确定正则化项中的双参数, 从而得到弯角梯度, 进而利用最大梯度法确定边界层顶高度.文中展示了两个数值实验的结果, 从中可知, 此方法最主要的优点在于, 尽管掩星资料中带有误差, 但其基本不受噪声的影响, 能够很好地提取主要的信息, 具有很好的稳定性和准确性.

应用本文提出的方法, 基于2007—2011年的COSMIC掩星弯角资料, 研究了全球海洋边界层顶高度的季节变化特征, 着重分析了 1, 4, 7, 10 四个月份的年平均边界层顶高度, 结果显示, 随着太阳直射点的南北移动, 海洋上的边界层顶高度具有明显的季节变化特征, 其与海温密切相关, 具体来说, 相对于周围海域的海温越低, 边界层顶高度越低, 反之, 边界层顶高度越高, 延伸这一结论可以得到, 冷洋流对应着低边界层顶, 暖洋流对应着高边界层顶.为了验证结论的准确性, 将其与利用掩星资料自带的zbalmax得出来的边界层顶高度的季节特征作对比, 结果发现, 两者边界层顶高度的时空分布特征十分一致, 说明通过双参数模型函数法得到的季节特征具有准确性.

总之, 经过数值实验中与差分法和L曲线法的比较, 发现双参数模型函数法在求边界层顶高度时具有不可替代的优势, 而且经过掩星资料自带的边界层顶高度数据zbalmax的验证, 表明文中表述的海上边界层顶高度的季节特征具有准确性.我们下一步需要继续做的工作是进一步改进方法, 使其适用于复杂天气条件下的大气波导的研究, 以期得到更好的结果.

本文利用的数据来自CDAAC提供的GPS掩星资料,在此谨表感谢! 感谢国防科技大学气象海洋学院周则明教授对文章行文结构的建议, 感谢国防科技大学气象海洋学院杜华栋副教授和闫申博士在编程方面的帮助, 感谢国防科技大学气象海洋学院陈毓敏、徐海波和北京师范大学地理学院孙铭扬对文章提供的修改建议.