引言

近年来随着我国气象观测网建设的快速发展，自动气象站观测数据时空密度大大提高，并且在大力发展气象无缝隙精细化智能网格预报技术背景下和人类生产活动对基于位置的气象服务要求下，将高时空分辨率不规则、离散站点实况资料快速、准确生成规则的格点资料具有重要的社会服务和业务应用价值。

将站点资料利用数据融合分析技术生成格点资料主要有遥感反演、资料同化、客观分析和统计插值等方法，由于受卫星空间分辨率、反演算法和复杂下垫面的制约，遥感反演的精度还有待提高^[1-3]，而后三种方法在气象领域应用较广泛，如美国NOAA/ESRL(National Oceanic and Atmospheric Administration/Earth System Research Laboratory)开发的局地分析与预报系统LAPS(Local Analysis and Prediction System)，后来升级为时空多尺度分析系统STMAS(Space-Time Multiscale Analysis System)，方法也从Barnes客观分析转变为连续变分资料同化，通过对不同尺度的观测信息进行分析能够消除长波信息与短波信息之间的干扰，得到更精确、细致的格点分析场^[4-6]。张涛等^[7]通过对比分析LAPS和STMAS地面2 m气温融合效果，结果表明在东部观测密集区LAPS和STMAS都有不错的表现，在资料稀疏地区STMAS的优势更为明显。李超等^[8]以LAPS系统本地化为基础研发了一套地面、高空多要素客观分析产品，结果显示气温的融合效果最好。师春香等^[9-10]基于LAPS/STMAS系统研制完成的中国气象局陆面数据同化系统CLDAS陆面要素融合分析产品，与中、外同类产品比较在中国区域的时、空分辨率和质量更高^[11-13]，在CLDAS-V2.0基础上，高分辨率中国气象局陆面数据同化系统(HRCLDAS-V1.0)也已经投入试运行，产品分辨率提高至1 km^[14]，产品质量有了进一步提高。但是，地形高度调整的不合理性对LAPS、STMAS在山区的融合效果有所影响，而且STMAS方法中将背景误差协方差和观测误差协方差的比值简单定义为常数，较为简单。张璐等^[15]将多重网格策略引入非线性最小二乘三维变分同化方法NLS-3DVar(Non-linear Least Squares-based on Three-dimensional Variational Data Assimilation)，并应用于2 m气温数据融合，NLS-3DVar产品相比于STMAS和Cressman插值产品来说，其具有更高的精度。何斌等^[16]利用Multiquadric和Cressman插值法对浙江省自动气象站的2 m气温资料进行客观分析，结果表明当站点较密集时两种方法分析结果差异非常小。邝建新等^[17]对广东省自动气象站的2 m气温、压强、风场、雨量进行质量控制后，采用Barnes方法进行客观分析，建立一套自动气象站资料客观分析的业务运行系统。何险峰等^[18]提出KNN、Barnes重采样技术为主线，梯度提升树(GBT)回归最优化的客观分析方法，实现了对传统气象Barnes客观分析、k最邻近法、Spark机器学习框架的改进与融合，生成1 km高分辨全要素格点场客观分析产品。马诺等^[19]利用协同克里金、统一高度普通克里金、多元线性回归3种空间插值方法对天山区域的气温进行插值，统一高度普通克里金方法虽然考虑了地形偏差的影响，但是认为气温垂直递减率是一常数(6.5 ℃·km^-1)。以上研究表明在资料密集的平原地区不同客观分析方法得到的产品质量相差不大，接近站点观测值。但是在站点稀疏和地形复杂区域，由于不同客观分析方法使用的背景场和地形处理方案不同，产品质量存在明显的差异。

而Cressman方法作为气象领域应用广泛且比较成熟的客观分析方法，是将站点视为同一平面上分布的离散点，其忽略了地形高度等影响因素，在降水客观分析方面取得较为理想的分析结果，但是对于2 m气温这种受地形高度偏差影响比较大的气象要素，Cressman方法较难获得高精度的估算结果。因此本文基于山东省自动气象站分布较密集结合地形特征，采用统一高度Cressman客观分析方法、地形订正等关键技术，对山东省2 m气温进行客观分析生成逐1 h、0.01°×0.01°高分辨率格点产品，供全省业务应用和社会服务。

1 资料和质量控制 1.1 资料介绍

本文使用的资料主要有2020年5—12月山东省及周边省自动气象站(国家级气象观测站、区域气象观测站、海岛观测站、浮标观测站等，以下统称“自动站”)观测资料，包含逐小时地面2 m气温、逐1 h累计降水量和逐小时地面气压; 欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts，ECMWF)预报的海上2 m气温和海平面气压，水平分辨率为0.125°×0.125°; 1970—2019年50 a间山东省123个国家级气象观测站2 m月最高温和月最低温资料; 山东省中尺度数值预报系统(WRF-RUC)预报的垂直位温递减率和地面气压，水平分辨率为3 km; 分辨率为90 m的地形高程数据SRTM(shuttle radar topography mission)。

Cressman客观分析时由于边界观测资料较少会引起较大分析误差，因此Cressman客观分析所用的离散站点区域要大于生成格点产品的区域范围，图 1中白色d01框(34.0°~39.0°N，114.0°~123.0°E)为不规则离散站点分布，红色d02框(34.25°~38.55°N，114.65°~122.80°E)为Cressman客观分析生成的格点场，陆地上蓝色离散点是自动站位置，约有3 500个站，海洋上白色点是ECMWF预报的格点位置，约有900个格点，鲁中山区和半岛地区自动站密度分布稀疏，平原地区自动站密度分布较密集。由于海上自动站资料较匮乏，根据多模式的检验结果，ECMWF预报的地面2 m气温平均误差和均方根误差最小，因此使用ECMWF预报的2 m气温作为补充，和自动站2 m气温组成Cressman客观分析所用的离散站点，约有4 400个站。离散的站点经过Cressman客观分析后输出一个纬度方向431个格点，经度方向816个格点，分辨率为0.01°×0.01°的等经纬度格点产品。

1.2 质量控制

自动站观测资料易受到站点周边环境、人工维护、数据传输、仪器设备等多种因素的影响导致自动站观测数据出现缺失或异常等现象，严重影响了气象观测数据的质量和应用。因此Cressman客观分析前需要对观测资料进行质量控制。本文对自动站观测资料主要使用极值检查、背景场检查、时间一致性检查、区域一致性检查常规的质量控制方法来剔除异常和可疑的数据^[20-24]。

1.2.1 极值检查

极值检查主要用来剔除明显异常的气温数据，利用1970—2019年50 a间山东省123个国家级自动站观测的月最高温和月最低温作为极值检查的阈值指标(表 1)，为了避免极端天气发生时的实况资料丢失，极值检查的阈值范围比统计阈值范围略大，当观测的自动站2 m气温高于或低于这个阈值指标时剔除。

1.2.2 时间一致性检查

时间一致性检查主要是检查地面2 m气温随时间的变化是否为一个连续变化的过程，通过检查观测数据时间变化率，剔除不真实的跳跃或无变化值。根据山东省地面2 m气温1 h变温的统计结果，本文时间一致性检查分为：当自动站有降水天气现象发生时，前1 h温差绝对值大于12.0 ℃时剔除该自动站观测资料; 当自动站没有降水发生时，前1 h温差绝对值大于8.0 ℃时剔除该自动站观测资料。

1.2.3 区域一致性检查

自动站空间分布距离相近的站点的特征值比距离较远的站点相关性要大，因此利用单站资料对该站进行质量控制的基础上再结合邻近台站的观测资料进行质量控制能够进一步提高观测资料的准确率，区域一致性检查就是充分利用同时刻邻近的多个台站观测资料，对该站观测资料进行质量控制的方法^[24]。考虑气象要素随海拔高度的垂直变化、气象要素距离相关性的衰减等原因，本文在地形订正后再进行区域一致性检查，首先将地面2 m气温利用地形订正方法订正到海平面高度并转换为位温，然后采用箱线图异常值的方法进行区域一致性检查，影响半径取10 km，若订正后的位温大于上限(Q₃+1.5×I_QR)或小于下限(Q₁-1.5×I_QR)，均视为异常值，其中，Q₃为75%分位数，Q₁为25%分位数，I_QR=Q₃-Q₁。

1.2.4 背景场检查

数值预报产品能够反映三维大气的时空变化演变特征，利用数值预报产品对观测资料进行质量控制，能有效剔除不符合气候态的观测资料。将山东省中尺度数值预报系统(WRF-RUC)预报的地面2 m气温作为背景场(B)，利用反距离权重方法插值到自动站，和自动站观测的2 m气温(O)计算观测余差(O-B)，当自动站观测的降水这一天气现象和数值预报相同，且|O-B|大于8.0 ℃时剔除该自动站观测资料，当自动站降水这一天气现象和数值预报不一致时，即数值模式降水出现空报和漏报时，|O-B|大于12.0 ℃时剔除该自动站观测资料。

2 方法 2.1 地形订正方法

统计表明地面2 m气温和地形高度呈负相关并具有统计学上的依赖关系，地形每升高100 m温度降低0.65 ℃左右^[25]，但当有地面逆温、焚风效应等发生时，温度梯度不再是常数。因此本文利用山东省中尺度数值模式预报的位温递减率，借鉴WRFDA变分同化系统中地面2 m气温地形矫正方法进行地形订正，公式如下：

$ {\theta _{r'}} = {\theta _1} + \frac{{{\theta _2} - {\theta _1}}}{{{h_2} - {h_1}}}\left( {{h_r} - {h_1}} \right) $

(1)

$ {\theta _{o\prime }} = {\theta _1} + \frac{{{\theta _2} - {\theta _1}}}{{{h_2} - {h_1}}}\left( {{h_o} - {h_1}} \right) $

(2)

$ \Delta \theta = \left( {{\theta _o} - {\theta _{o\prime }}} \right)\frac{{{h_2} - {h_r}}}{{{h_2} - {h_o}}} $

(3)

$ {\theta _r} = {\theta _{r\prime }} + \Delta \theta $

(4)

其中，h_o为订正前地形高度，θ_o、θ_o′为高度h_o处位温和插值位温，h_r为订正后地形高度，θ_r、θ_r′为高度h_r处位温和插值位温，h₁、h₂为距地面气压差为100 hPa和200 hPa模式层的高度，h₁、h₂的位温分别为$ {\theta _1}、{\theta _2}, \frac{{{\theta _2} - {\theta _1}}}{{{h_2} - {h_1}}}$为位温递减率。

首先利用位温公式计算地面2 m位温θ_o、距地面气压差为100 hPa和200 hPa模式层的位温θ₁、θ₂，然后利用$ \frac{{{\theta _2} - {\theta _1}}}{{{h_2} - {h_1}}}$位温递减率计算h_r的插值位温θ_r′和地面2 m插值位温θ_o′(公式(1)、(2))，接着计算地面2 m位温θ_o与插值位温θ_o′的增量并将其转换到订正后高度(公式(3))，最后地形订正后的位温为θ_r(公式(4))，并利用位温公式得到地形订正后的2 m气温。

2.2 Cressman方法

Cressman客观分析方法由CRESSMAN^[26]于1959年首先提出，至今依然用于各种气象资料的空间插值。Cressman客观分析方法采用逐步订正方法将离散点(站点)插值到规则格点上。先假定格点上初猜场的值，本文使用扫描半径R范围内的区域平均值作为初猜场，然后用扫描半径内的观测资料与初猜场的分析增量(ΔG_ij)订正初猜场，得到新的初猜场，逐步迭代直到订正的初猜场逼近观测资料为止。通常，一次客观分析过程需要设置多个扫描半径R进行多次Cressman分析才能取得较好的分析结果，公式如下：

$ G_{ij}^n = G_{ij}^{n - 1} + \Delta {G_{ij}} $

(5)

$ \Delta {G_{ij}} = \frac{{\sum\limits_{k = 1}^k {{W_{ij, k}}} \left( {{O_k} - G_{ij}^{n - 1}} \right)}}{{\sum\limits_{k = 1}^k {{W_{ij, k}}} }} $

(6)

$ \left\{ \begin{array}{l} {W_{ij, k}} = \frac{{{R^2} - {r^2}_{ij, k}}}{{{R^2} + {r^2}_{ij, k}}}, 当{r_{ij, k}} \le R\\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;当{r_{ij, k}} > R \end{array} \right. $

(7)

$ G_{ij}^{(n, R)} = {V^R}G_{ij}^{(n, R)} + \left( {1 - {V^R}} \right)G_{ij}^{(n, R)} $

(8)

其中，R是扫描半径，r_{ij, k}是扫描半径范围内第k个观测与格点(i, j)的距离，W_{ij, k}是第k个观测与格点(i, j)的影响权重，O_k是扫描半径范围内的第k个观测值，G_ij是初猜场的值，G_ijⁿ为第n-1次迭代的Cressman分析值。V^R为扫描半径R的Cressman分析权重，V和R的维数相同，G_ij^{(n, R)}为扫描半径R的Cressman分析结果，G_ij^{(n, R-1)}扫描半径R-1的分析结果。本文经过多组对比优化试验，分析权重和扫描半径为V=(0.9, 0.9, 0.9，0.8, 0.8, 0.7, 0.7)，R=(0.20, 0.16, 0.12，0.09, 0.07, 0.05, 0.03)。

2.3 统一高度Cressman方法

山东省自动站观测网分布较密集，但山区和地形复杂区域尤其是山顶观测资料仍然较少且分布不均匀，如鲁中山区和山东半岛区域(图 1)，这些区域站点的时空代表性较差，因此需要对这些区域进行地形订正减少地形对气象要素插值引起的误差。另外，由于常规Cressman方法客观分析时没有考虑地形对地面2 m气温的影响，造成气温等值线和地形不匹配这显然是不合理的，因此本文采用统一高度Cressman方法对山东省地面2 m气温进行客观分析。步骤(图 2红色框)如下：首先利用反距离权重插值方法将WRF-RUC预报的垂直位温递减率插值到离散站点位置，然后利用基于位温递减率的地形订正方法将地面2 m气温转换为海平面高度上的位温，接着对位温进行Cressman客观分析生成0.01°×0.01°海平面高度位温格点分析场，最后再利用基于位温递减率地形订正方法将海平面位温订正到真实地形高度(SRTM)生成地面2 m气温格点分析场。

2.4 检验办法

本文检验使用了平均误差、平均绝对误差、均方根误差、预报准确率检验指标，公式如下：

$ 平均误差:ME = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {{F_i} - {O_i}} \right)} $

(9)

$ 平均绝对误差:MAE = \frac{1}{N}\sum\limits_{i = 1}^N {\left| {{F_i} - {O_i}} \right|} $

(10)

$ 均方根误差:RMSE = \sqrt {\frac{1}{N}{{\sum\limits_{i = 1}^N {\left( {{F_i} - {O_i}} \right)} }^2}} $

(11)

$ 预报准确率:PC = \frac{{NA}}{{NB}} \times 100\% $

(12)

其中，O_i为第i站(次)实况温度，F_i为第i站(次)客观分析温度，NA为 F_i-O_i ≤2 ℃或 F_i-O_i ≤1 ℃的站(次)数，NB为客观分析或实况的总站(次)数。PC温度预报准确率指温度客观分析绝对误差≤2 ℃或1 ℃的百分比。

3 结果分析 3.1 地形订正影响

两种客观分析方法的主要区别是统一高度Cressman方法考虑了地形偏差对2 m气温的影响。图 3为2020年5月14日14时山东省鲁中山区2 m气温格点产品等值线图，图 3a为常规Cressman方法客观分析结果，图 3b为统一高度Cressman方法客观分析结果，检验时两种产品都没有进行临近站点替换。通过对比，利用统一高度Cressman方法客观分析后，地面2 m气温等值线和地形结构走向一致，气温随高度上升温度降低，地形梯度越大的地方温度梯度也越大。而常规Cressman客观分析的地面2 m气温等值线在鲁中山区较平滑。如徂徕山(图 3a、b红色圆圈)海拔高度为1 028 m，两种Cressman方法客观分析时都没有使用该站观测资料，显然统一高度Cressman方法客观分析结果更合理。

表 2为2020年5月山东省2 m气温Cressman客观分析格点产品检验结果，Test1为常规Cressman方法，Test2为统一高度Cressman方法。Test2、Test1的平均误差分别为-0.008 6 ℃、-0.162 9 ℃，降低了0.154 3 ℃，平均绝对误差分别为0.159 7 ℃、0.417 2 ℃，降低了0.257 5 ℃，均方根误差分别为0.353 7 ℃、0.614 7 ℃，降低了0.261 0 ℃，2 ℃以内准确率分别为99.60%、98.14%，提高了1.46%，1 ℃以内准确率分别为98.09%、96.51%，提高了1.58%。结果表明统一高度Cressman方法相比常规Cressman方法提高了客观分析产品的质量，尤其在没有资料和地形复杂区域，客观分析产品的结果也更加合理。

3.2 检验与分析

利用中国气象局陆面数据同化系统CLDAS和HRCLDAS地面2 m气温格点实况产品，和统一高度Cressman客观分析格点产品(没有进行临近站点替换)进行对比，HRCLDAS产品2020年7月底业务试运行。HRCLDAS和统一高度Cressman产品水平分辨率为1 km，自动站和格点距离较近，因此取临近格点和自动站实况资料进行匹配，CLDAS产品水平分辨率为5 km，自动站和格点距离较远，因此利用双线性插值到山东省1 635个自动站上和实况资料进行匹配。通过分析2020年5—12月CLDAS、HRCLDAS和统一高度Cressman格点产品地面2 m气温的检验结果(表 3)，三种产品平均误差分别为-0.025 6 ℃、-0.020 0 ℃和-0.003 9 ℃，均为负偏差，这主要是由于插值对极值平滑而引起的，平均绝对误差分别为0.516 0 ℃、0.174 0 ℃和0.146 9 ℃，均方根误差分别为0.872 3 ℃、0.407 7 ℃和0.359 7 ℃，2 ℃以内准确率分别为97.24%、99.62%和99.64%，1 ℃以内准确率分别为86.40%、97.43%和98.24%。总体上CLDAS产品不如HRCLDAS和统一高度Cressman产品，这主要是因为CLDAS产品分辨率较低，CLDAS产品临近站点的格点和自动站距离较远，而HRCLDAS和统一高度Cressman产品临近站点的格点和自动站点距离较近甚至重合; HRCLDAS产品8月误差较大，各项检验指标均不如统一高度Cressman产品，9月统一高度Cressman产品质量略优于HRCLDAS产品，10、11和12月HRCLDAS产品质量除了平均误差外其他各项检验指标均略优于统一高度Cressman产品，以上分析结果表明了通过统一高度Cressman方法客观分析得到的产品质量接近HRCLDAS产品质量。

图 4为2020年5—12月统一高度Cressman方法客观分析的山东省地面2 m气温月平均误差空间分布，通过对每个月的对比分析，月平均误差基本在±1 ℃以内，鲁西北、鲁西南、鲁东南、山东半岛等平原地区月平均误差较小，多数在±0.5 ℃以内。鲁中山区地形高度较高的站月平均误差较大，多数在±0.5~±1.0 ℃之间，泰山、崂山月平均误差最大，主要是由于在地形复杂区域，自动站空间分布不均匀，站点高度不一致，温度梯度较大引起的。气温偏低的10、11、12月温度准确率均略低于5、6、7、8、9月，这可能是由于冬季空气的相对湿度偏低，热容量小，湍流交换较弱，大气层结较稳定，温度受局地辐射和地形的影响使冬季温度局地性相比夏季偏强引起的，因此统一高度Cressman影响半径和影响权重系数应根据不同的季节优化不同的参数。

4 业务应用

利用统一高度Cressman方法生成的地面2 m气温格点分析产品，临近自动站的格点和自动站较近甚至重合，如果自动站实况和临近格点的偏差较大，说明客观分析在此格点可能不合理，例如从地面2 m气温月平均误差分布(图 4)来看，仍然存在误差较大格点，因此利用经过质量控制后的自动站2 m气温临近替换格点可以矫正分析不合理的格点，最终生成临近替换后的山东省逐1 h、0.01°×0.01°高时空分辨率地面2 m气温格点产品并业务应用。图 5为2020年5月20日14时经过临近替换后的客观分析格点产品在Micaps软件上的显示效果，格点产品在山东省不同区域，如鲁中山区、海岸线、海岛、平原地区，地面2 m气温等值线和地形匹配合理，并且和自动站观测无偏差。图 6为稀疏化后格点图像产品，该产品在山东省灾害性预警平台业务应用，自业务试运行以来，系统运行稳定、产品无缺测，系统收到实况资料后在2 min内运算完成并生成格点图像和数据产品。

5 结论与讨论

本文将地形订正方法和常规Cressman方法相结合形成了统一高度Cressman客观分析方法，这种处理方式解决了常规Cressman方法在地形复杂区域客观分析的格点产品和地形不匹配问题。将统一高度Cressman方法应用到山东省地面2 m气温客观分析中，生成了逐1 h、0.01°×0.01°高分辨率地面2 m气温格点产品。经检验统一高度Cressman客观分析格点产品质量接近HRCLDAS格点实况产品，取得了较理想的分析结果。

1) 利用山东省中尺度数值模式预报的位温递减率，借鉴WRFDA变分同化系统中地面2 m气温地形矫正方法进行地形订正，通过和常规Cressman方法进行对比分析，表明统一高度Cressman方法相比常规Cressman方法提高了客观分析产品的质量，并且在没有资料和地形复杂区域，客观分析产品的结果更加合理。

2) 将改进后的统一高度Cressman方法应用到山东省地面2 m气温客观分析中，进行多级尺度Cressman客观分析，该方案经多次试验后，确定了较佳的Cressman客观分析影响半径和权重系数，获得较为合理的温度梯度分布并使临近格点逼近自动站观测值。

3) 评估了自2020年5月业务试运行以来统一高度Cressman客观分析格点产品的质量，全省平均误差基本在±1 ℃以内，平均误差分布鲁中山区地形高度较高的区域略大于鲁西北、鲁西南、鲁东南、山东半岛等平原地区。5—12月1 ℃以内准确率为98.24%，温度偏低的10、11、12月温度准确率均低于5、6、7、8、9月。