引言

气象上，通常需要借助温度场来分析对流层低层暖空气和中高层冷空气及其之间的配比，需要在对流层低层分析暖脊，在中层分析冷槽(张杰，2006；朱乾根，2007)。暖脊是指从高温区中延伸出来的较狭长区域，通常用暖脊线来表征，即从暖中心出发，沿系列等温度线曲率极大值点勾勒而成。在我国地域，暖脊线通常呈南到北、西到东或西南到东北的走向(张小玲等，2010)。在我国的气象数据资料分析处理中，暖脊线的识别和绘制是通过人工或借助计算机软件(如我国的MICAPS-气象信息综合分析处理系统; 李月安等，2010)提供的辅助工具人工完成的。其绘制质量难免受到人为因素(如经验不足)的影响，本文试图通过构建合理的算法，实现暖脊点的自动提取及暖脊线的自动绘制，并为基于天气图的自动分析、挖掘出更多的客观天气规律奠定基础。

在地质领域，从数字化地形资料中自动提取山脊线和山谷线的研究相对成熟，常见的有适于数字化等高线资料的二维算法和适于高程网格数据的三维算法，分别称为等高线曲率极值法和地形断面极值法(黄培之等，2005)，两者均是从数据中提取脊线或谷线上的候选点，再根据一定的规则对这些候选点进行分析和归类。如果将对流层等压面上的温度值等同于地形高度，暖脊线与山脊线便具有了类似的几何性质。

1 基于断面极值法的暖脊特征点分析

断面极值法从求取数据格点的梯度入手，令梯度方向为过数据点的垂直断面的法线方向，那么，等值线最大曲率点就是此断面上格点数据的局部极值点。

设曲面函数为

$ Z = f\left({x, y} \right) $

(1)

如果我们关注函数值的变化最大方向，可求梯度：

$ \nabla = \mathit{\boldsymbol{n}}\left({{f_x}, {f_y}} \right) $

(2)

式中，$ {f_x} = \partial f\left({x, y} \right)/\partial x, {f_y} = \partial f\left({x, y} \right)/\partial y $。此梯度矢量的方向即为函数值的变化最大方向：

$ {\theta _0} = \arctan \left({{f_y}/{f_x}} \right) $

(3)

对于离散的格点数据，可用差分代替微分。设格点(i, j)上的数据为Z_{i, j}，其八邻域的格点坐标及其取值如图 1所示，式(3) 中的f_y和f_x的三阶反距离平方权差分计算方法(陈楠等，2007)如下：

$ \left\{ \begin{array}{l} {f_y} = \frac{{\left({{Z_7} - {Z_1}} \right) + 2\left({{Z_8} - {Z_2}} \right) + \left({{Z_9} - {Z_3}} \right)}}{{8{\Delta _y}}}\\ {f_x} = \frac{{\left({{Z_1} - {Z_3}} \right) + 2\left({{Z_4} - {Z_6}} \right) + \left({{Z_7} - {Z_9}} \right)}}{{8{\Delta _x}}} \end{array} \right. $

(4)

将θ₀的计算范围(－90°, 90°)映射到(0°, 360°)，并将其离散成8个方向，八邻域近似断面如表 1所示。

图 2a是根据MICAPS提供的850 hPa下的温度格点数据绘制的等温线和利用断面极值法得到的暖脊候选特征点标识，根据我国温度场中暖脊的走向，只需关注0°~90°的特征点，○表示0°特征点，●表示45°特征点，▲表示90°特征点。图 2b展现了天津市气象台预报员卢焕珍在同一张等温线图上手工绘制的暖脊线。通过对比，不难发现：

(1) 在真实存在的四条暖脊线附近有一些暖脊候选特征点被提取出来，特别称这些点为有效点。

(2) 有较多的候选特征点不能形成对暖脊线的支撑，称这些点为无效点。就图 2a而言，无效点占全部候选特征点的比例高达52.3%(46/88)。

将图 2a中的有效点和无效点进行对比发现(参见点A和点B)，有效点附近等温线两侧下滑迅速，无效点附近等温线下降平缓且等温线密集。在此，特别称前者为陡峰点，后者为坡峰点。

(3) 在有效点中，有些点的方向与暖脊线走向不一致(差45°)，如图 2a中的点C，对其做局部水平剖线，将剖线上的各点温度值如图 3，可以发现，这类点两侧温度值变化一缓一急，称这类点为坡-陡峰点。面对坡-陡峰点，断面极值法给出的特征点对暖脊而言出现了两点偏差，即方向偏差和位置偏差。

需要强调的是，上述定义的陡峰点和坡峰点是针对特征点处等值线形态而言的，而坡-陡峰点则指的是特征点处水平或垂直方向温度值沿其两侧下降的特点。

2 暖脊特征点识别 2.1 对坡-陡峰点进行纠偏

针对上述方向偏差和位置偏差，本文将单侧梯度扩展到双侧梯度，同时扩大计算中心极值点两侧梯度的范围。首先需要将传统断面极值法划分的梯度方向范围进行扩充，为叙述方便，首先给出如下定义：

定义1脊宽：在断面极值点两侧或一侧连续存在略低于该极值点的系列点形成准极值区。

定义2显著性脊：在脊的准极值区两侧温度值迅速下降的脊称为显著性脊。

定义3伸展性特征点：在特征点两侧，等值线呈迅速下滑态势的特征点称为伸展性暖脊特征点。显然，前述定义的陡峰点即为伸展性暖脊特征点，由多个伸展性暖脊特征点构成的暖脊线一定是沿距离缓慢下降的，而由多个坡峰点构成的“脊”线则沿距离下降迅速，反映为一组密集的等值线，如图 4所示。

下面以90°暖脊线上左坡右陡型坡-陡峰点为例叙述纠偏算法。

(1) 求准极值区确定特征点处的脊宽

① 扫描第i行数据，对梯度值落入[45°, 135°]范围内的点做过该点的长度为n的水平剖线。在剖线范围内找到极大值，将极大值点(i, j)标记为水平坐标轴的原点O，令原点右侧为正方向；

② 分别将(i, j+1) 和(i, j-1) 的温度插值到另一侧的等数值处，获得插值距离d_R和d_L，令r＝int{max(d_R, d_L)}。例如，图 3中点C的水平断面数据如图 5所示，(i, j+1)点向左的插值距离d_R＝2.947，(i, j-1)点向右的插值距离d_L=0.029，显然，r＝2；

③ 设等温线分辨率为ΔT，若Z_{i, j}－Z_{i, j-r}≥ΔT，则r－1→r，回到③，否则继续；

④ 将(i, j－r)点的温度向右插值求出插值距离d′_L，若d′_L≥0.5，则c＝1，否则c＝0；令集合S∈[(i, j－r), (i, j+c)]为准极值区；

⑤ 为确保区域S上的各点温度值相近：找到i+Δ₁和i+Δ₂，使Z_{i+Δ1, j-r}＝Z_{i+Δ2, j+c}＝Z_{i, j}，若max{Δ₁, Δ₂}＞ρ_T(默认值取0.5·ΔT)，则令r－1→r，c－1→c，回到⑤，否则计算此处的脊宽：s＝r+c+1。

(2) 显著性判断

在满足单调递减条件的情况下，将暖脊两侧各外扩d个单位，若：

$ \max \left\{ {{Z_{i, j}} - {Z_{i, j - r - d}}, {Z_{i, j}} - {Z_{i, j + c + d}}} \right\} > {\alpha _T} $

则宽度为s的脊是显著性的，保留显著性特征点、舍弃非显著性特征点(本文设定α_T＝1.5℃)。

(3) 将脊宽的中心确定为显著性脊的特征点(i, j₀)，其中，j₀＝j+(c－r)/2。

图 6展现了本文纠偏算法的效果，对比可见，图 2a中传统断面极值法给出的19个坡-陡型峰点的角度及位置偏差全部得以纠正。同时，非坡-陡型峰点依然被正确地提取出来，这是因为，面对坡-坡型或陡-陡型的峰点，算法结果是j＝0，c＝0，s＝1，j₀＝j。另外，此算法还在一定程度上降低了无效点的数量(52.3%降低到41.4%)。

2.2 滤除无效的坡峰点

由暖脊的定义可以看出，在温度值满足极大值阈值的情况下，等值线水平投影需同时满足曲率足够大的条件。由图 4可知，由坡峰点组成的脊伸展力度不够，应该滤除。在等温线图上，坡峰点附近的等温线曲率小、等温线密集(图 7)，温度下降坡度大(伸展度小，见图 4)，于是，在坡峰点所在脊宽的外部对满足单调递减条件的k个格点数据进行如下运算(以90°坡峰点为例)：

求$ \Delta {i_{Lk}}和\Delta {i_{Rk}}, k = 1, 2, 3 $，使$ {Z_{i + \Delta {i_{Lk}}, {j_0} - r - k}} = {Z_{i, j}}, {Z_{i + \Delta {i_{Rk}}, {j_0} + c + k}} = {Z_{i, j}} $

令$ \Delta {i_L} = \max \left\{ {\Delta {i_{Lk}}} \right\}, \Delta {i_R} = max\left\{ {\Delta {i_{Lk}}} \right\} $

若$ \min \left\{ {\Delta {i_L}, \Delta {i_R}} \right\} < 1.5\Delta T且\Delta {i_L} + \Delta {i_R} < 5\Delta T $，则(i, j₀)为无效的坡峰点，将其滤除。

2.3 传统断面极值法所得特征点的有条件保留

前两小节不论是显著性特征点的识别，还是坡峰点的滤除，均是针对90°和0°脊线而言的。考虑到45°暖脊线不是主流方向，首先借助传统断面极值法将它们检出，再附加判断条件，将错误的坡-陡点滤除，即只保留满足附加条件的45°特征点。

错误的45°特征点是其两侧梯度变化不均匀所致(即坡-陡峰点)，为此，本文将待考核点视为均值，利用方差概念构建约束量β

$ \beta = \arctan \left({{{f'}_y}/{{f'}_x}} \right) $

(5)

式中，f′_y、f′_x分别是八邻域格点数据纵向和横向方差，即

$ \left\{ \begin{array}{l} {{f'}_y} = \frac{{\sum\nolimits_{k = 2, 5, 8} {\left[ {{{\left({{Z_k} - {Z_{k - 1}}} \right)}^2} + {{\left({{Z_k} - {Z_{k + 1}}} \right)}^2}} \right]} }}{6}\\ {{f'}_x} = \frac{{\sum\nolimits_{k = 4, 5, 6} {\left[ {{{\left({{Z_k} - {Z_{k - 3}}} \right)}^2} + {{\left({{Z_k} - {Z_{k + 3}}} \right)}^2}} \right]} }}{6} \end{array} \right. $

(6)

显然，候选点纵横向方差相同时，β＝45°，两者相差越大，β越趋向0°或越趋向90°。将式(6) 作为保留45°特征点的附加条件，即若θ′₀(i, j)为45°且Z_{i, j}为135°方向断面上的极值，则进一步计算β，若β∈[30°, 60°]，则认为(i, j)点的45°方向为真，保留此点，否则舍弃之。

图 8是经过以上各步的暖脊特征点提取结果。与图 2a相比，错误的45°点全部被纠正、有效点由47.7%(42/88) 增加到91.9%(34/37)，无效点所占比例由52.3%下降到8.1%。

2.4 暖脊特征点识别流程

设某高度的温度格点数据规模为N×M，对格点数据Z_{i, j}(i＝2, …, N, j＝1, …, M)进入图 9所示流程。

3 暖脊线的自动绘制

在识别出暖脊特征点的基础上，绘制暖脊线的工作分为“定向搜索-形成折线”和“光滑中轴-拟合暖脊线”两个步骤。

3.1 定向搜索形成折线

鉴于暖脊线的方向限制在[0°，90°]范围内，按照从右到左、从上到下的顺序搜索特征点时，首先遇到的特征点一定是暖脊线的低值(尾端)。搜索算法设计如下：

(1) 在暖脊特征点集内找到未被标记过的最低点，将该点作为当前点并作标记(黄培之，2001)；

(2) 从当前点出发，根据其梯度方向确定搜索范围(3+Δx)×(3+Δy)如图 10所示，需要根据暖脊线的方向确定Δx和Δy的取值；

(3) 重复步骤(2)，直到搜索范围内未标记点空；

(4) 重复步骤(1)~(3)，将连续搜索到的暖脊特征点顺序连线，并记录起点P₁和终点P₂，计算暖脊长度L＝‖P₁－P₂‖；

(5) 遍历暖脊线，若某条暖脊线的终点和另一条暖脊线的起点距离小于等于Δ，且两者走向基本一致，则首尾相接，合并成一条新的暖脊线。

(6) 遍历各条暖脊线，若L≤3，则表明此区域暖脊强度较弱，删除此暖脊线。

对图 8执行上述算法的4条折线如图 11所示。

本文以暖脊线上多数特征点方向作为暖脊方向，以暖脊线的长度表示其强度。

3.2 光滑中轴拟合算法

人们通常使用最小二乘法用已知数据点拟合出一条光滑的曲线(倪慧等，2011)，这时，需要设定曲线的函数形式，而暖脊线因其多变性和不定性，试图建立统一有效的数学模型是不客观的，本文基于“找到某条光滑曲线，使其即反映数据点预示的变化趋势，又使光滑曲线对数据点拟合的误差平方和最小”的最小二乘设计思想，构建光滑中轴算法，其中，中轴概念如图 12所示。

图中，设A、B、…、I为定向跟踪算法得到的属于一条暖脊线的暖脊特征点，虚线为辅助线，b、c分别为AC和BD的中点，以此类推，顺势连接Bb中点、Cc中点、Dd中点……并将端点适当延长，即可得到均匀地穿过离散数据点的中轴折线，再用二次B样条对中轴折线做进一步平滑处理，得到光滑的中轴。经此算法处理得到的平滑缓变的暖脊线如图 13所示。

4 测试及分析

为了验证本文算法的正确性和有效性，随机采集2005－2013年间30组温度场数据(850 hPa)进行测试。图 14是2005年7月31日08时的数据样本。图 14a中的四条暖脊线由气象专家在MICAPS系统上手工绘制而成，图 14b的背景是以三角网格法得到的等温线(戴泽军等，2003；张利等，2009)，基于本文算法的全部暖脊线识别结果自动绘于其上。

表 2是对30组样本检验后的统计结果。可以看出，针对测试样本，本文算法的暖脊线击中率高达96.91%，且未出现误识，凡击中的暖脊线，在一定程度上比人工分析更精细，且方便提供量化参数；漏识原因分析：本文算法基于局部思想而设计，识别暖脊线时欠缺顾全整体，在特征点稀疏处，容易出现一个暖脊系统被分裂为两个及以上的情况，如果分暖脊线长度较短，会因其强度不够被删除。

5 结论和讨论

针对天气图分析的实际需求和人工绘制暖脊线效率低、质量参差不齐的问题，本文进行了暖脊线自动识别方法的研究。将地形中山脊线的识别方法进一步改进用于暖脊识别，改进措施是根据对流层温度场的特点进行了改进，提出的“纠偏算法”将极值点放宽到准极值区，有效解决了面对坡-陡峰特征点时传统断面极值法出现的方向偏差和位置偏差问题。在函数模型未知的情况下，根据最小平方误差的思想，提出基于离散数据点进行曲线拟合的“光滑中轴算法”，使自动绘制的暖脊线符合气象业务规范与气象专家手工绘制的暖脊线吻合较高。对于测试样本，在未出现误识的前提下获得了接近97%的击中率。

暖脊自动识别的实现一则可以协助业务员快速准确找出格点数据所反映出来的具有气象意义的特征，二则特别适合于从多年的气象资料中，挖掘出与暖脊相关的客观气象信息。再就是，将该方法做适当改进，可以实现高度场上脊线的自动识别，以及副高、南亚高压等系统脊线的自动识别。