引言

全球无线电探空资料最早可以追溯到1938年，是一种可以用于气候研究的长序列资料。近年来，很多学者开始逐渐注意到探空资料中人为因素造成的不均一性对气候变化研究的影响。我国高空测站建站时间较早，最早出现在1951年，大部分探空站从1957年陆续开始建站。探空资料非均一性研究从20世纪90年代开始进行，但只有少数台站进行了试验性研究(翟盘茂，1997)。近期对全国多数探空温度序列的非均一性分析证实了辐射订正方法改变、仪器换型等是造成我国探空温度序列非均一性的主要原因(郭艳君，2008；郭艳君等，2008)。然而，通过无线电探空得到的位势高度资料也是一个重要的站点探空要素。根据《常规高空气象探测规范》(中国气象局监测网络司，2002)中高空探测计算公式的规定，位势高度由两个相邻等压面间的厚度由低到高逐渐累加得到，其中温度要素也参与了计算，因此温度值偏差将导致位势高度值的偏差，且由于位势高度采用逐层累加的计算方法，这一偏差会逐层累积，因此无线电探空过程中由于温度探空造成的资料不均一性会影响到位势高度要素(马颖等，2010)。另外，马颖等(2010)的研究也指出：由于我国采用简化的订正方法等原因，导致59型探空观测数据与预报场的相对系统偏差较大，到20世纪90年代100 hPa位势高度的偏差达到60 m左右，被WMO点名的台站超过50%。因此，中国气象局于2000—2001年分两批对全国所有台站59型探空仪的误差订正方法进行了一次全面的修正。此后，59型探空仪系统偏差大为改善。

从上面的分析可以看到，由于探空温度要素引入的辐射订正方法改变、仪器型号更换，以及由于探空业务系统中位势高度算法参数以及本身订正方法的改变都可能造成位势高度要素的不均一。究竟上述因素在多大程度上影响我国的探空位势高度资料？这方面的研究工作还非常少(马颖等，2010)。本文将利用加拿大环境部气候研究中心研发的PMFT(Wang, 2008)非均一性检验方法，结合国家气象信息中心存档的中国区域各探空台站详细的元数据信息对中国区域1951—2008年123个探空台站各标准等压面层的月平均位势高度资料进行非均一性检验和分析。

1 采用的资料、元数据以及均一化方法 1.1 资料

我国探空站点总体上分布比较均一，但在青藏高原地区站点稀少。本工作进行非均一性检验与订正的1951—2008年月平均位势高度数据来自国家气象信息中心经过质量控制、质量状况较好的全国123个探空站。探空台站每天分两个时次进行探测，分别为世界时的00和12时(下同)。涉及到的标准等压面层包括100、200、300、400、500、700和850 hPa。其中，经过质量控制认为错误的每天两次探测的位势高度资料不参加月平均值的统计，这部分错误资料占总资料的0.2015%。

1.2 元数据

本部分所用的元数据信息来源于国家气象信息中心制作的中国探空台站历史沿革，其中包括台站仪器变更信息和台站迁移信息。辐射订正方法、位势高度计算参数及订正方法的变化(系统升级)信息来源于中国气象局气象探测中心。

1.3 非均一性检验与订正方法

用于气候资料非均一性检验的方法很多(刘佳等，2012)，本文采用的惩罚最大F检验(PMFT)方法是加拿大环境部气候研究中心建立的一种均一化方法(Wang，2008)。该方法经验性的考虑了时间序列的滞后一阶自相关，并嵌入多元线性回归算法，能够用于检验、订正包含一阶自回归误差数据序列的多个间断点(平均突变)。该方法不需要使用参考序列，因此本工作直接对原始探空位势高度资料进行非均一性检验。该方法首先通过回归检验算法来检验出多个间断点，检验过程简述如下：

假设ε_t代表变量，该变量的均值为0，方差为σ²，对于存在的线性趋势β，其时间序列为{X_t}，要检查在t=k时刻是否存在一个平均突变，

假设：

$ {H_0}:{X_t} = \mu + \beta {\rm{ }}t + {\varepsilon _t},{\rm{ }}1,2, \ldots ,N $

(1)

备选假设：

$ {H_a}\left\{ \begin{array}{l} {X_t} = {\mu _1} + \beta t + {\varepsilon _t},{\rm{ }}t \le k\\ {X_t} = {\mu _2} + \beta t + {\varepsilon _t},{\rm{ }}k - 1 \le t \le N \end{array} \right. $

(2)

μ₁≠μ₂，当H_a为真时，t=k时的点被称作间断点。Δ=|μ₁-μ₂|被称作平均突变的大小，最可能的间断点服从以下分布：

$ P{F_{\max }} = \mathop {\max }\limits_{1 \le k \le N - 1} \left[ {P\left( k \right){F_c}\left( k \right)} \right] $

(3)

其中，P(k)是建立的经验性的惩罚因子，其建立方法可参考(Wang, 2008)。

$ \begin{array}{l} \begin{array}{*{20}{l}} {{F_c}\left( k \right) = \frac{{\left( {SS{E_0} - SS{E_A}} \right)}}{{SS{E_A}/\left( {N - 3} \right)}}}\\ {SS{E_A} = \sum\limits_{t = 1}^k {{{\left( {{X_t} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \mu } }_1} - \beta t} \right)}^2}} + } \end{array}\\ \sum\limits_{t = k + 1}^N {{{({X_t} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \mu } }_2} - \beta t)}^2}} \\ SS{E_0} = \sum\limits_{t = 1}^N {{{\left( {{X_t} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \mu } }_0} - {\beta _0}t} \right)}^2}} \end{array} $

(4)

其中$ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \mu } _0}$和β₀是在μ₁=μ₂=μ时的估计值。

通过回归检验算法来检验出多个间断点，首先在t={N_min, N_min+1, …, (N-N_min)}时段找出最可能的间断点C₀。类似的，分别在时段t={N_min, …, (C-N_min)}和t={(C₀+1+N_min), …, (N-N_min)}找出次可能的间断点C₀₁和C₀₃，接下来在C₀₁和C₀₃之间的时段找出最可能的间断点C₀₂，然后在C₀₁、C₀₂和C₀₃之间寻找最可能的间断点，即依次分段找出序列中各段最可能的间断点，计算所有间断点的统计量PF_max，估计自相关和P的大小，找出最大的PF_max值对应的间断点，如果是显著的，该突变就被认为是找到的第一个间断点，该间断点被记为N_c=1。寻找该间断点位置之后每段最可能的间断点，估计其显著性，找出下一个可能的间断点，重复该过程，将间断点按照显著性由大到小排列，形成间断点列表。

该方法可以使用元数据，当利用回归检验检测出多个间断点后，输入元数据，再次判断输入的时间点是否存在间断点，即依次分段找出序列中各段最可能的间断点，计算所有间断点的统计量，确定第一个间断点，寻找该间断点位置之后每段最可能的其他间断点，估计其显著性，找出下一个可能的间断点，重复该过程，分步找出所有的间断点。将间断点按照显著性由大到小排列，形成间断点列表，判断最小的间断点是否显著，当不显著时剔除该间断点，再次评估剩余间断点的显著性，最终保留的统计显著的间断点即为序列检验得到的间断点。该方法能够检验的最短序列长度是20(Wang, 2008；曹丽娟等，2010；曾红玲等，2010)。

另外，对于元数据的应用方面，由于通过几年的元数据上报和核查工作，国家气象信息中心已经较全面地掌握了我国探空台站的元数据信息，所以本文主要采用PMFT方法和元数据信息相结合的方法判断资料断点，最终保留的断点都是有元数据信息支持的断点，做到了绝不没有根据地订正资料。

在最终确定断点后，主要采用了平均值订正法对断点进行订正，即以序列的最新资料为标准, 把断点前后的平均值差作为偏差, 对断点前序列进行订正。

统计最终采纳的断点数占PMFT方法检测的所有显著断点数的百分比，为84.3%，也就是说有15.7%的显著断点没有被采纳，没有采纳这些断点的主要原因之一是这些没有被采纳的显著断点中大部分是由于气候本身突变造成的资料不连续。由于适合我国探空资料的参考序列比较难找到，所以本工作中并没有使用参考序列，只是利用PMFT方法这种不需要参考序列的断点检测方法对原始资料进行断点的检测，最后结合详细的元数据进行断点的确定，这种方法在断点的检测过程中会有一些由于气候突变造成的资料不连续点被检测出，比如：1982年前后和1997前后的ENSO事件，以及1991年前后Pinatubo火山爆发等，虽然这些事件对探空温度资料的影响较大，但是探空位势高度资料是计算量，计算过程中需要用到温度这个要素，所以这些影响温度的气候突变事件也对探空位势高度要素造成了一定的影响。在断点检测过程中，尤其是300 hPa以上的位势高度序列中可以发现这些时间点前后一年内的显著断点，但是在没有元数据支持的情况下，这些断点没有在最后的结果中被采纳。另外，有2.47%左右的显著断点没有元数据支持，通过与台站反复核查这些时间点上是否有元数据的变化，并没有找到这方面的记录，本着宁可不订正也不错误地订正断点的原则，最终这些断点也没有保留。

曹丽娟等(2010)利用PMFT方法对中国地面风速资料进行了非均一性检验，其结果也表明PMFT方法可以有效地检测出中国区域753个地面台站中存在356个风速资料的间断点，其中88.5%的间断点是有元数据支持的。在结果的讨论过程中，曹丽娟等(2010)也指出其结果还需要利用更加详细的元数据信息进行讨论和分析。在PFMT方法的使用方面，由于本工作前期已经对元数据信息进行了逐台站的收集、整理，并进行了大量详细的元数据核查工作，因此，本文的断点判断方法与曹丽娟等(2010)所不同的是，本工作最终保留的间断点都是有元数据信息支持的间断点。

2 结果分析 2.1 个例分析

图 1给出邢台站订正前后年平均位势高度序列。按照1.3介绍的PMFT方法，首先利用回归检验的方法检测出该序列存在1988年4月和2000年4月两个间断点。该站的元数据信息显示，邢台站在1978年6月进行过仪器换型，在2000年7月进行过系统升级。因此按照元数据信息，加入一个可能的间断点1978年6月，另外，由于允许检测出的断点和元数据信息有前后12个月的偏差，因此对于回归检验检测出的2000年4月的间断点，把这个时间点修改成为2000年7月。重新评估这些间断点的显著性，1978年6月的间断点不显著，删除该间断点；再次进行检验、评估剩余间断点的显著性，1988年4月和2000年7月的间断点都是显著的。联系台站并查阅大量的历史材料，但是都没有找到1988年4月该站发生过元数据信息变化的记载，最终这个间断点没有保留。邢台站最终保留的一个间断点出现在2000年7月，该站在这个时间点进行了探测系统的升级，升级后的位势高度值明显降低。图 2为西沙站00时次100 hPa位势高度序列。首先利用回归检验初步检测出该序列存在1958年9月、1999年2月、2001年1月和2007年3月这4个间断点。核查元数据，西沙站在1959年1月、1966年1月和1999年进行过辐射误差订正方法改变、2001年1月进行了系统升级、2006年3月进行了仪器的换型。因此，加入1966年1月和2006年3月这两个间断点，把1958年9月这个间断点改为了1959年1月，进一步重新评估这些间断点的显著性。结果得到1966年1月、2006年3月和2007年3月这3个断点不显著。首先删除显著性最小的2006年3月这个间断点，重新评估剩余间断点的显著性。得到1966年1月的间断点不显著且显著性最小，删除该时间点，继续重新评估，得到2007年3月这个间断点不显著，继续删除该时间点，最终判断该序列存在3个显著的间断点，分别位于1959年1月、1999年2月和2001年1月。可以看出，采用该方法能够有效地检测出这3次由于人为因素造成的资料不连续。

2.2 月平均探空位势高度资料订正比例

图 3是两个时次各个标准等压面上均一的台站数的变化情况，从图中可以看出，00和12时次均在850 hPa上位势高度序列均一的台站最多，位势高度序列均一的台站数分别占所有被检验台站的49.59%和52.02%；00和12时次均在100 hPa位势高度序列均一的台站数最少，分别为26和32个，分别占被检验台站的21.1%和26.0%。Lanzante等(2003)指出探空仪的误差分布规律与大气密度有关，空气越稀薄探空仪误差越大。另外，高层仪器受到太阳辐射的影响也更为明显，所以辐射订正方法改变对高层探空位势高度资料的影响也是造成高层探空位势高度资料均一的台站数较低层少的一个原因。从图 3中可以看出，300 hPa以上12时次均一的台站数量明显多于00时次，这主要是由于12时次不受到太阳辐射的影响，因此也不存在由于辐射订正方法改变造成的资料断点，相应位势高度序列均一的台站数也会增加。

图 4为各标准等压面层各探空位势高度序列中存在1~7个断点所占的比例。从图中可以看出，从850~200 hPa各序列主要以1~2个断点为主，而100 hPa各位势高度序列则以3个断点为主，这与100 hPa受仪器误差累计效应、辐射订正方法改变等的影响较大密切相关。

2.3 探空位势高度资料的订正幅度

这里定义某探空站在某标准等压面上的订正值为：这个台站订正后多年平均的位势高度值减去订正前多年平均的位势高度值。某标准等压面上平均订正值的计算规定为：该标准等压面上不均一的位势高度序列订正值的算术平均。把订正值小于零的订正简称为负订正，订正值大于零的订正简称为正订正。图 5为2个时次各个标准等压面上平均订正量的变化图。从图中可以看出，2个时次位势高度的订正值随气压的降低呈负增加的趋势。100 hPa 2个时次123个台站平均订正量可以达到-40.9 gpm。中国气象局分别在2000年7月(50站)和2001年1月(其余站)进行了系统升级，全面修正59型探空仪记录的误差订正方法。修正前，与欧洲中心初估场相比，100 hPa高度上的位势高度系统偏高40~60 gpm，这也从一方面说明了本文的订正结果较为合理。另外，从图 5中可以看出，850~100 hPa上所有标准等压面的订正量均为负值，显示了我国探空位势高度序列存在着系统性偏高的问题。

表 1统计的是订正前后2个时次位势高度的差值(12时次值减00时次值)：订正前100~500 hPa 12时次平均位势高度值要大于00时次的平均位势高度值，订正后这种趋势没有改变，但是两个时次的位势高度差明显减小。而700和850 hPa订正前00时次平均位势高度值大于12时次平均位势高度值，订正后这种趋势没有改变，但是两个时次的位势高度差值有所增加。

2.4 订正前后月平均位势高度探空资料趋势分析

由于我国大部分台站从1957年开始建站，因此图 6中参加统计的时间段为1958—2008年。从图 6中各标准等压面位势高度资料的变化趋势上看，123个探空台站在对流层顶和对流层高层(300~100 hPa)平均位势高度序列订正前位势高度呈降低趋势，平均为-4.03 gpm·(10 a)^-1。而订正后呈位势高度增加趋势，为2.66 gpm·(10 a)^-1。在对流层中低层(850~400 hPa)原始位势高度序列呈现增高趋势，为1.24 gpm·(10 a)^-1，订正后增高趋势更加明显，为2.68 gpm·(10 a)^-1。

2.5 月平均探空位势高度资料非均一问题产生原因

我国于1966和1999年分别进行了辐射订正方法的更改，1957年4月探空观测时间从世界时的03和15时变更为00和12时，2000和2001年进行了探测业务系统的升级(包括算法参数的变化和误差订正方法的变化)，并且不定期进行探测仪器的换型等，这些原因都是导致位势高度序列不均一的主要原因。图 7为4种因素在全国123个探空台站各标准等压面上导致位势高度序列间断点个数的变化图，由于辐射误差只在300 hPa以上层次才做订正，所以图中由于辐射订正方法改变造成的断点只存在于100、200和300 hPa。从图中可以看出仪器换型和探测系统升级导致间断点的数量基本上随着高度的增加而趋于增加。00时次3种因素均在100 hPa导致的间断点数量最多，辐射订正、仪器更换和探测系统升级产生的间断点个数分别为170、184和137个。而由于观测时间改变造成的间断点个数非常少，这是因为我国探空台站的建站时间集中在1957年以后，而观测时间改变发生在1957年4月，所以受到观测时间变化影响的台站相应也非常少。仪器更换、探测系统升级以及辐射订正方法改变是造成中国区域探空位势高度资料不均一的主要原因。

3 结论

通过本文的统计和分析，可以得到如下结论：

(1) 仪器更换、探测系统升级以及辐射订正方法改变是造成中国区域月平均探空位势高度资料不均一的主要原因。

(2) 世界时00和12时次位势高度序列不均一的台站数随高度的增加而增加。

(3) 各个标准等压面上位势高度序列不均一的台站平均订正幅度随着位势高度的升高而增大。并且两个观测时次各个标准等压面上位势高度序列均为负订正，表明我国探空台站各标准等压面位势高度资料存在系统性偏高的问题。

(4) 123个探空台站在对流层顶和对流层高层的位势高度序列订正前平均为位势高度降低趋势，订正后位势高度序列呈增高趋势。对流层中低层原始位势高度序列呈位势高度增高趋势，订正后这种增高趋势更加明显。