喀什站高空温度资料均一性检验与订正

2019-03-17胡义成张正阳王秋香

沙漠与绿洲气象 2019年6期

胡义成，张正阳，王秋香*

（1.中国气象局乌鲁木齐沙漠气象研究所，新疆乌鲁木齐830002；2.新疆气象信息中心，新疆乌鲁木齐830002）

目前存档的探空资料中蕴含了气候变化的重要信息，但是由于探空资料中人为因素造成的不均一性和气候因素造成的变化混淆，使得研究者很难利用探空资料准确地判断气候变化趋势。为了使探空资料很好地运用到气候研究中，很多学者致力于检验和订正探空资料中时间序列的不均一性。和地面气象要素均一性检验相类似的是[1-9]，对于探空资料不均一性的检验需要识别由于人为因素造成的间断点，估计间断点偏离大小并对间断点进行订正[10]。我国高空观测站建站时间较早，最早出现在1950 年，大部分高空站从1960 年陆续开始建站，站点遍布全国。均一性研究在20 世纪90 年代开始进行，但只有少数台站进行了试验性研究[10]。近期对全国多数探空温度序列的均一性分析证实了辐射订正、仪器换型等方法的修改造成了序列的非均一性，这种非均一性在各站探空资料中不同程度的存在，均一化过程也存在着相当的不确定性[11，12]。

喀什高空气象观测站（以下简称“喀什站”）经历了多次仪器变更、台站位置变动、观测规范变更、观测时间变化以及辐射订正方法变化（以下简称：辐射订正）。这些活动是否对高空资料的均一性产生了影响，高空资料是否可以连续使用，都是值得去研究和探讨的问题。本文利用加拿大环境部气候研究中心研发的RHtest[13，14]均一化检验方法，对喀什站1958—2016 年各个标准等压面层的年平均温度资料进行均一性检验与订正，对喀什站高空资料的完整性、连续性和均一性进行了分析，有利于帮助台站提高单站高空资料的质量，也可以为新疆高空气象观测业务的运行提供参考。

1 资料方法与观测环境

1.1 气象观测环境变迁及现状

喀什站始建于1954 年7 月，自1957 年6 月1日承担着北京时间08、20 时常规探空和测风的同步观测任务。建站以来，喀什站址多次变动。1955 年5月1 日因不明原因发生了站址位置变动（图1）；1957年7 月1 日，站址由疏附县浩罕庄迁至喀什市西南郊；1964 年1 月1 日又因不明原因导致站址经纬度信息发生变化；因2014 年1 月1 日观测环境已远远不能满足气象观测的要求。

图1 历次迁站过程图

1.2 资料与方法

使用的基础资料来自国家气象信息中心“中国高空规定等压面定时值数据集（V2.1）”及高空台站元数据集。资料时段为喀什站1958—2016 年，包括北京时间08、20 时（世界时00 时和12 时）各规定等压面位势高度、气温、露点温度观测数据和风向风速。数据集中1951—2010 年数据来源于高空观测纸质月报表数字化资料，2010—2016 年12 月的数据来源于高空全月观测数据文件（高空月报G 文件）。该数据集在研制过程中借鉴了《无线电探空资料质量控制（QX/T 123-2011）》中相关质量控制方法及参数，以多源数据对比检测为主，对可疑资料实施人工核查与订正。不仅更正了纸质资料本身存在的观测和抄录错误数据，还更正了信息化过程导致的资料错误和无数据现象，最大程度地确保了研究数据的可靠性。

参考序列是ERA-Interim，ERA-Interim 是欧洲中期天气预报中心（ECMWF）发布的最新的全球大气数值预报再分析资料，它是继早期资料（如ERA-15，ERA-40）后的新产品，目的是对ERA-40 等资料进行完善，从而逐渐取代ERA-40。ERA-Interim所采用的同化系统是ECMWF 集成预报系统（IFS，32R2 版本的全球谱模式）。相对于ERA-40 在很多的关键部分都进行了改进，包括模式的改进、四维空间变量分析的运用、对卫星资料变量偏差的修正以及对观测系统中其他数据的处理等。

RHtestsV5 软件包是在RHtestsV4 分位数匹配订正（Quantile-Matching adjustments，QM）中增加使用参考序列的基础上，在序列检验中提供了一阶滞后自相关，还提供了数据序列的均值订正、QM 订正和序列曲线图以及结果的回归模拟。可以对日或者小时数据序列进行QM 订正。检验方法包括两种，一种是惩罚最大t 检验（PMT）方法，其检验过程中需要建立参考序列，待检序列与参考序列的差值（差值序列）是被检验的对象；另一种是惩罚最大F 检验（PMF），适用于无参考序列的检验过程。

PMT 方法检验过程如下：

假设｛Xt｝（t=1，…，N）代表一个ⅡD 正态分布序列，要检验时间序列｛Xt｝是否存在断点。

备择假设：

μ1≠μ2，并且“ ｛ Xt｝~‖DN（μ，σ2）”代表｛ Xt｝遵循正态分布，其均值为μ，方差为σ2。当Ha为真时，t=k 时的点被称作间断点，Δ=被称作平均突变的大小，最可能的间断点服从以下分布：

其中：P（κ）是建立的经验性的惩罚因子，其建立方法可参考文献[13]。

通过回归检验算法来检验出多个间断点，首先在t∈｛Nmin，Nmin+1，…，（N-Nmin）｝时段找出最可能的间断点C0，Nmin指的是时间序列中被分割片段的最小长度，也就是说，两个相邻断点之间子序列的长度或者第一（最后一）个断点到序列起点（终点）之间的长度，类似的，分别在时段t={Nmin，…，（C0-Nmin）}和t={（C0+1+Nmin），…，（N-Nmin）}找出次可能的间断点C01和C03，接下来在C01和C03之间的时段找出最可能的间断点C02，然后在C01，C02和C03之间寻找最可能的间断点，即依次分段找出序列中各段最可能的间断点，计算所有间断点的统计量PTmax，估计lag-1 自相关和p值的大小，找出最大的PTmax值对应的间断点，如果是显著的，该突变就被认为是找到的第一个间断点，该间断点被记为NC=1。寻找该间断点位置之后每段最可能的间断点，估计其显著性，找出下一个可能的间断点，重复该过程，将间断点按照显著性由大到小排列，形成间断点列表，判断最小的间断点是否显著，当不显著时剔除该间断点，需要再次评估剩余间断点的显著性，最终保留都统计显著的间断点即为序列检验得到的间断点，受分段检验的影响，该方法能够检验的最短序列长度是20，详细过程见文献[14]。

PMF 方法检验过程如下：

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

备选假设：

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

其中，P（k）是建立的经验性的惩罚因子，其建立方法可参考[14]。

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

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

RHtests 方法已被人们广泛地应用于气温、降水、风速及相对湿度等气候序列的均一性分析中。因喀什站高空各层次日温度资料有大量缺测值，且待检序列（1958—2016 年）与对应的ERA-Interim 参考序列（1979—2013 年）时间不一致，在对喀什站高空日、年温度资料进行均一性检验和订正的过程中发现，在利用RHtestsV5 软件包进行订正时，大量缺省值的代入造成了订正序列与实际观测序列的偏差量过大，导致订正序列不可信。故本文利用RHtestsV5 软件包提供的PMT 方法，采用对应的ERA-Interim 作为参考序列，对喀什站高空日温度资料进行均一性检验和订正；利用RHtestsV5 软件包提供的PMF 方法，对喀什站高空年温度资料进行均一性检验和订正。

1.3 月观测资料序列完整性

喀什站常规探空和测风的同步观测任务起始于1957 年6 月1 日，以有整年观测数据（1958 年1 月）以来的平均温度序列按照850、700、600、500、400、300、200 hPa 和100 hPa 共8 个规定等压面的月值（共708 个月）来统计缺测率及数据缺测年月分布状况，统计结果见表1 和表2，由表可见，喀什站8 个规定层1967 年9—10 月均缺测。

根据表1 和表2 给出的缺测统计结果，08 时8个层次的月平均温度除100 hPa 的缺测率>3%，完整性较差外，其他层次完整性较好；20 时100 hPa的月平均气温完整性较差，200 hPa 的完整性一般，其他6 个层次的完整性较好。综合分析认为喀什站850、700、600、500、400 hPa 和300 hPa 共6 个层次的高空温度观测资料序列完整性较好。

表1 喀什站高空温度序列月值缺测情况（08 时）

表2 喀什站高空温度序列月值缺测情况（20 时）

2 观测资料序列均一性分析

2.1 日观测资料均一性分析

利用RHtestsV5 软件包提供的PMT 方法，采用对应的ERA-Interim 资料作为参考序列，分别对喀什站08 时、20 时对应的850、700、500、400、300、200 hPa 和100 hPa 共14 条（因搜集的ERA-Interim资料中无对应的600 hPa 序列，故此处没有对喀什站600 hPa 的日气温序列进行断点检测和分析。）温度序列作置信度为0.95 的PMT 检验。保留其中的两类断点，一类是在年或月序列中检验出的显著间断点，并且有元数据的支持，如果该断点出现的时间与台站元数据记录信息相差一年以内，根据元数据记录时间替换该断点位置；另一类是没有元数据支持，但是在日和月（或年）序列中被同时检验出的相同年份的显著间断点，断点位置依据日平均温度序列被检验出的断点时间。

检验得到的断点信息显示，喀什站14 个温度序列被检验出的断点共有221 个（图2），其中有元数据支持的断点占89.14%。从检测结果和元数据分析情况看，喀什站始建于1954 年7 月，而20 世纪50 年代以来中国高空观测仪器和规范变化较大。1957 年4 月高空观测时间从北京时间11 时和23 时变更为08 时和20 时。1959 年前后、1966 年以及1999 年7 月1 日后又经历了辐射订正方法的变化，1962 年5 月31 日、1968 年10 月1 日、1970年7 月1 日和2005 年喀什站又陆续更换探空仪器。综合分析认为，喀什站温度资料非均一性主要来源于台站位置变动、观测仪器换型、观测规范变更和辐射订正。

图2 喀什站温度序列检测出的断点数量

通过PMT 方法和元数据信息相结合的方法来综合判断资料的断点，所确定的断点均通过了0.05的显著性检验。其中，有一部分比例的断点虽然有元数据支持，但调查数据时发现由于各种原因，在断点处、前或后的数据为缺测。如，14 条数据序列均在1967 年9 月12 日—1967 年11 月10 日缺测，而检测出得到的1967 年9 月12 日和1967 年11 月10日两个断点及其间缺测的数据不进行订正。另外，在对日原始数据序列进行订正时，仅对有参考序列对应的时段进行订正。通过判断，20 时700 hPa、08时500 hPa 和20 时100 hPa 的非均一性最为明显，分别检测出32 个、26 个和25 个断点（表3），经查询台站元数据信息，这3 条序列断点的元数据支持率分别为84.38%、57.69%和88%，综合考虑以上各种情况，这3 条序列最终确定需要订正的断点分别为2 个、6 个和2 个，占所检测出的断点比例仅分别为6.25%、23.08%和8%，做到了绝不无根据地订正资料。

表3 2 个观测时次7 个标准等压面上断点数

通过检验出的各规定层气温序列断点的位置和订正幅度来看，14 条序列均检验出了断点，除了08时100 hPa 以外，其他13 条序列均检验出一个相同位置（1997-12-09）的断点，经查询台站元数据文件发现，在上世纪90 年代前后，喀什站经历了辐射订正方法的变化。由此可以认为，辐射订正方法的变化导致13 条序列产生了间断点。另外，08 时850~400、100 hPa 和20 时400 hPa 共7 条序列因台站位置的变动（1980-05-01 和1993-07-01）而产生了间断点。在利用RHtestsV5 软件包进行断点检测和订正时，序列中的缺测值需用缺省值“-999.99”来代替。在对各序列进行订正时发现，300 hPa 以上各序列的订正幅度非常大，其中20 时100 hPa 的日气温序列订正量高达-184 ℃，这么大的订正量是不可信的。通过查找原因，发现300~100 hPa 这6 条日值序列的缺测率较高（表4），而08 时和20 时各序列订正量的大小与对应的序列缺测率之间呈显著的正相关关系，相关系数分别高达0.97 和0.99，由此不难理解，大量缺省值的代入，会导致订正序列与观测序列的偏差过大。为了保证订正序列的可信度，本研究仅对缺测率小于1%的温度日值序列进行订正，即仅对08 时和20 时850~400 hPa 共8 条序列进行了订正，具体订正情况见表5。

2.2 年平均气温序列均一性分析

利用RHtestsV5 软件包中的最大惩罚F 检验（PMFT）方法对喀什站1958 年以来观测的850~100 hPa 共14 条平均气温的年值序列进行断点检验。该方法经验性的考虑了时间序列的滞后一阶自相关，并嵌入多元线性回归算法，能够用于检验、订正包含一阶自回归误差数据序列的多个间断点（平均突变），而且不需要使用参考序列。

表4 各规定等压面日气温序列缺测率

表5 日温度序列检验出的断点位置和订正量

针对PMFT 方法检测出的年序列中的显著间断点，基于喀什站详实的高空历史沿革记录，挑出仪器变更、台站位置变动、观测规范变更、观测时间变化以及辐射订正方法变化等有记录的时间点（表6），来判断在该时间是否为间断点。如果该断点出现的时间与台站元数据记录信息相差一年以内，根据元数据记录时间替换该断点位置。若检测出的间断点不显著，也没有元数据支持，则舍弃该间断点，也不进行订正，做到绝不无根据地订正资料。

表6 喀什站存在变化的时间点及原因

年平均气温资料序列均一性检验结果如图3～8所示。

图3 喀什站08 时100 hPa 年平均气温检测结果

（不连续点：1997 年，辐射订正方法变化）

图4 喀什站08 时300 hPa 年平均气温检测结果

图5 喀什站08 时500 hPa 年平均气温检测结果

图6 喀什站08 时700 hPa 年平均气温检测结果

图7 喀什站08 时850 hPa 年平均气温检测结果

（不连续点：1967 年，两个月无观测；1977 年，海拔高度变化、观测手册变化；1997 年，辐射订正方法变化；2005 年，雷达和探空仪换型、L 波段系统运行、数据文件换型、观测规范变化）

图8 喀什站20 时100 hPa 年平均气温检测结果

喀什站14 条温度序列中有6 条序列即08 时100a、300、500、700、850 hPa 和20 时100 hPa 存在不连续点，其中，08 时850 hPa 的不连续点高达4个，20 时除100 hPa 外，其他序列均无不连续点。结合元数据分析，造成喀什站观测资料不连续的原因主要有：1967 年有近2 个月停止观测，1974 年观测规范变化，1977 年海拔高度变化和观测手册变化，1997 年辐射订正方法的变化，2005 年雷达和探空仪换型、L 波段系统运行、数据文件换型、观测规范变化。其中，08 时850 hPa 年平均温度在2005 年出现的断点，很可能是由于L 波段与59-701 系统的感应元件和测温原理不同而存在的器差所致[15]。

2.3 年气温订正幅度及订正前后趋势分析

经过均一性分析和趋势分析发现，喀什站14 条高空温度序列中，08 时850、700、500、300、100 hPa和20 时100 hPa 共6 条序列不连续，需要订正，其他8 条序列连续，则不需订正。通过PMFT 方法和元数据信息相结合的方法确定断点后，主要采用了RHtestsV5 软件包中的平均值订正法对最终被采纳的断点进行了订正，即以序列的最新资料为标准，把断点前后的平均值差作为偏差，对断点前序列进行订正，订正结果见表7。最终采纳的断点数占PMFT方法检测出断点数的百分比为76.9%，也就是说有23.1%的断点没有被采纳。

通过订正，6 条气温序列订正后与订正前均呈现出相反的变化趋势，其中，08 时100、300 hPa 和500 hPa 均为副订正，而且订正前气温均为下降趋势，订正后的气温均为上升趋势；08 时700 hPa 和20 时100 hPa 则为正订正，订正后的气温均为下降趋势。可以很明显地发现订正幅度随着规定层高度的变化呈现出两头大中间小的现象，其中08 时和20 时100 hPa 的订正幅度分别为-1.13 ℃和-0.80 ℃，而08 时850 hPa 的订正幅度则高达2.56℃。08 时需要订正的序列要多于20 时，08 时7 条序列中有6 条序列做了订正，订正率达到了85.71%，而20 时仅有100 hPa 做了订正。