参考文献

ACTA Scientiarum Naturalium Universitatis Pekinensis - - Contents -

全光相机是一种深度相机, 又称为光场相机,能够通过单次拍摄来提供场景的三维信息。随着Lytro 公司和 Raytrix 公司研制成功手持式相机, 利用全光影像进行三维重建(深度信息)已成为研究热点[1]。虽然 Lytro 公司和 Raytrix 公司的相机结构分别属于全光1.0相机和全光2.0相机[2], 但这两款相机均有记录光线强度和方向信息的能力, 这也是利用全光相机进行三维重建的基础[3]。

对于Lytro相机, 通常进行深度估计的方式是将

微透镜影像分解成图像阵列的形式[4–5], 图像阵列的行列数目等同于单个微透镜下覆盖的探测器单元数目, 而每一幅影像的分辨率等同于微透镜的总数目。为了获取深度图, 一些研究[6–10]利用传统的立体视觉特征点匹配的方式进行视差计算, 将图像阵列的中心影像作为参考影像。另外一些研究利用核线平面影像[11–15]检测同名点构成的直线斜率, 利用直线斜率与深度值一一对应的关系, 计算得到深度

[16]图 。这两种方式计算得到的深度影像均可以利

用图割Graph Cuts[17–18]和置信度传播bp (belief propagation)[19–20]等消除误匹配。但是,由于受光学结构的限制, 计算出的深度图与参考影像的分辨率相等(只有原始影像的1%)。

通常, Raytrix相机微透镜的数目更少, 使用上述方法进行深度计算会产生分辨率更低的深度图,因此需要基于微透镜阵列直接进行深度计算。在基于微透镜直接计算深度的研究中, Georgiev等[21]提出一种利用归一化交叉相关算法, 沿着横向和纵向的方式, 进行同名点匹配搜索计算视差, 最终采用WTA (winner-take-all)方法得到深度图。之后, Zeller等[2223]采用同样的思路, 但是通过强度差异计算深度, 同时改进了视差搜索方式, 主要沿着核线方向搜索同名点, 最终通过多次观测去除异常深度信息,

[24]得到稀疏的点云信息。Fleischmann等 提出一种不同的利用传统视觉的方法, 首先以单个微透镜为基础, 构建代价立方体CV (cost volume), 然后采用

[25]半全局方式 对单个微透镜的CV进行优化, 消除误匹配。由于微透镜的边缘限定了优化的范围, 因此这种方式有待改进。本文提出一种全局优化深度图的方法, 深度图的构建不会受到微透镜数目的限制。首先, 基于微透镜进行代价计算[22–24], CV的高度由虚拟深度的范围决定[26–28],同时仍然保留高精度的WTA深度结果; 然后, 将CV通过获取的初步深度信息投影到像空间, 同时构建新的CV, 并利用全局BP优化产生平滑且可靠的深度图。这里获取的深度图是离散的, 因此我们会将高精度且连续的WTA深度图结果一并投影, 并利用此结果对上一步获取的离散深度图进行优化, 通过马尔科夫随机场传播(Markov Random Fields Propagation, MRFSP)[12,29]获取既连续也平滑的深度图。为验证方法的效果, 通过实验与文献[24]的方法进行对比。MRFSP优化深度图的流程如图1所示。

1 光场三维信息恢复方法 1.1 光学结构分析

全光相机包含3个主要部件: 主透镜、微透镜阵列以及探测器阵列。主透镜将物方空间成像至像方空间, 并且物方景深被压缩到一个合理的像方景深, 如图 2所示。可以将微透镜阵列和探测器阵列视为一个多相机阵列系统, 对像方空间进行观测,这时就可以利用多视几何相关理论进行三维信息计算, 并获得深度a的值。探测器与微透镜之间的距离 b 通常是一个固定值, 因此将虚拟深度v作为深度计算的最终结果[26]:

v  a / b 。(1)在实际计算过程中, 微透镜被视为一个小孔成像模型, 因此在微透镜–探测器组成的观测系统中,深度计算类似于双目结构, 如图3所示。

虚拟深度v由下式计算得到: l l

v  , (2) xp1  xp2 x

其中, xp1 和 xp2 是点 p1 和 p2相对于微透镜的局部像点坐标, x为视差。根据式(2), 可以得到虚拟深度的变化与视差变化的关系:

v2  v   ( x)。 (3)

l式(3)表明, 为视差计算精度 ( x)在一定的范围时,越长的基线能够保障越高的虚拟深度精度, 这时选择基线长的虚拟深度值有利于提高深度图的精度。

1.2 代价计算策略

为了获得视差信息, 首先需要进行同名点的匹配来构建 CV, 匹配的具体过程如图4(a)所示, 计算时深度的离散值数目为DL , 如图5所示。本文采用经典的 SAD (sum of absolute difference)进行同名点

对的强度一致性度量:

ci ( m , n ,  x )  min(| I ( m , n )  I ( m m, SAD i

n  n )|,  1) (4)其中, 1是一个限定代价大小的常数, Ii 是待计算微透镜周围的第i个微透镜影像, 如图4(a)所示。(∆m, ∆n)是沿着基线方向的两个微透镜的偏移量。

2式(2)中 ∆x等于偏移量的欧式距离  m  n2 , 如图4(b)所示。此外, 考虑到鲁棒性, 视差计算采用自适应窗口[30]: ci ( m , n ,  x)  SAD   w ( p , q ) w ( p , q , I ) w ( p , q , Ii ) distance color color  p q

 w ci ( m , n ,  x)   ( p , q )w ( p , q , I ) SAD  distance color  pq w ( p , q, Ii ) , (5) color

p 2  q2 rd 在式(5)中, w ( p , q )  e和w ( p , q , I ) distance color e  I  m , n  I  m  p , n  q  / rc 分别是基于距离信息和颜色信息的权重, rd 和 rc 是常量。点( p , q)是相邻像素点的坐标。由于微透镜影像中的同一个点对于不同的基线会有不同的视差, 所以CV应该建立在虚拟深度上,因此最终的代价函数应该为1 c ( m , n , v )   ci ( m , n , x), (6) SAD SAD N i其中, N 是计算中涉及的微透镜的数目, v可以利用视差x通过式(2)计算得到。考虑到式(3)的影响, 一般利用式(6)计算代价信息会降低虚拟深度的精度。由于全光相机是一个多基线的系统, 且视差值与虚拟深度之间并不是简单的线性关系, 若要保持精度只能提高CV中的深度方向的数量。为了平衡深度信息的精度以及深度方向的数量, 我们在利用式(6)进行深度计算的同时, 保留高精度的WTA深度结果(该结果直接采用视差进行深度计算), 最后的深度图由式(2)转化至虚拟深度, 并用于后续的优化。深度计算的流程见算法1(其中符号 round() 表示取整)。算法1 基于微透镜代价及WTA深度图构建算法。1 输入: 全光数据2 针对每一幅微透镜影像执行步骤3~5 3 对于微透镜中每一个点( m , n) , 搜索与之最近的N 幅微透镜影像, 执行步骤4~5 1 4 计算代价 c ( m , n , v )   ci ( m , n , x) , 其SAD N SAD i中, v  round( li / x) 5 如果 ci ( m , n , x) 是所有可能的视差中最小的, SAD则将 WTA 赋值深度 deplenslet (m , n ) li / x 6 输出: CSAD 和 WTA 深度图dep lenslet

1.3 深度与代价投影

由1.2节计算得到的CSAD只能在单个微透镜影像内部进行局部优化[24], 因此本节对CV和WTA深度信息进行投影, 建立新的基于投影影像的CV和 深度图。投影之前, 首先利用CV信息进行深度可靠性的估计, 建立基于微透镜影像的可信度图[20]。可信度计算公式如下:

c second ( m , n )  c First ( m , n) conflenslet ( m , n) SAD SAD , (7) c second ( m , n) SAD

c First ( m , n) c second ( m , n)为点( m , n) 对应的最小代价值, SAD SAD为仅次于最小代价值的代价值, 如图5(b)所示。可信度 conflenslet ( m , n)的取值范围为[0,1) 。可以利用可信度图选择合适的代价曲线进行投影。参考图2和 3, 像方空间点 P′ 和图像点 pi (i=1, 2)以及透镜的中心满足共线条件, 因此投影到像方空间点的坐标通过下式确定:

 xp   xp   xc , s ' v   v 1    (8)    y '  y   y  P p c

其中,( x c, yc)是微透镜的中心坐标, s为决定投影影像的分辨率控制因子。显然, 投影获得的坐标( xp , yp ) 不是整数, 因

此本文采用插值的手段获取投影代价曲线和投影深度图。

  iw , y  y (11)  ( x x ) , distance P P P P   i i c xpi ypi v) )分( , ,、depproj ( xpi ypi , )以及confproj ( xpi ypi , proj 别代表投影代价值、深度值和可信度值, xp 和y的P i i取值在实数域连续。值得注意的是 X (xp , yp ) proj i i X ( xp , yp )( X 代表式(9)~(11)中{ c, conf , dep} )。lenslet由于同一个点有多条代价曲线以及深度值, 因此进行插值时, 只选取可信度最大的前N 条代价曲线, proj以此来提升鲁棒性。

1.4 深度图优化

经过投影步骤, 可以得到基于投影影像的CV,本文采用由粗到精的3层BP对其进行优化[1920]。由于BP不能保证聚合, 故设定多次迭代得到最后的结果。代价值的每一次迭代使用式(12)和(13)。

c( x , y , v )  c ( x , y , v )   mt ( v) , (12) new P  P  P  P  q p qn p mq t ( v )=min V ( v , vi )  c ( vi )  ms t 1 ( vi ), (13)  p i  q s  Nq , s p

t q p其中, m 表示传播的信息, 代表像点 邻域的

q  p像素(本文中, 相邻像素点之间的深度约束V ( v , vi )使用线性模型)。BP优化后的深度影像和可信度图计算结果如下:

depop ( xp , yp )  arc min( c ( xp , yp , v)) , (14)   v new   c second ( x , y ) c First ( x , y ) conf ( xp , yp )  new P P new P P 。 (15) OP  c second ( x , y ) new P P

利用式(14)得到的深度图depop 是离散的, 利用式(10)计算的深度图 depwta虽然包含误匹配信息但是连续的。为了结合两者的优势, 我们利用MRFSP模型[12,29]对深度图进行融合优化, 最终的深度图计算如下:

dep  arc min  conf (,)|(,) xy dxy d i i {WTA,OP} x ,y depi ( x , y )|  d  d , (16) flat 1 smooth 1  ( ) 代表梯度算子,  ( )代表拉普拉斯算子。详细的计算过程见算法2。算法 2 MRFSP 优化深度图算法。

1 输入: c , deplenslet 和 conflenslet SAD 2 利用式(7)~(11)进行投影, 计算得到 depwta 和confwta conf3利用式(12)~(14)进行全局BP优化, 获得depop和OP 4 利用式(16)的MRFSP 获取融合两类深度图5 输出: 优化后的深度图dep

2 实验

本文的数据集由Raytrix 42相机提供, 每个微透镜覆盖约31  31个像素单元。计算深度图前, 先用白光影像进行微透镜中心点的检测[4], 结果如图 6所示。

获取微透镜中心点信息后, 按照图 1的步骤构建深度图, 计算中使用的参数如表1所示。

此外, 使用式(5)时, 自适应窗口的大小为5 5像素, 而利用式(9)~(11)进行插值计算的窗口大小为 2  2。本实验分别从深度图连续性以及深度估计的鲁棒性两个方面与现有算法的对比进行分析。

2.1 深度图分析

3 类深度图分别建立在DL 为30, 50, 70的CV上,对比结果如图7所示。图7(B)中的结果为WTA深度图, 可以看到其中包含很多噪点, 但深度分布是连续的。用全局优化得到的深度图(图7(c), (e)和(g))可以消除这些噪点, 但是深度信息明显离散, 具有一定的阶梯性, 从放大的图像(图7(i))可以明显看到这种离散性。经过MRFSP优化后的结果(图7(D), (f)和(h))既可以消除噪点信息, 也能保持深度图的连续性。

为了进一步对比, 本文以DL=70的深度图作为参考影像, 对图(d), (f)和(h)进行相减操作, 并将差值归一化到[0, 1], 差值图像如图8所示。可以看出,经过MRFSP优化后, 深度图的主要差异只体现在边缘地区。利用DL=50计算得到的深度图已经趋近于DL =70的深度图, 因此在深度计算中选择DL =50

已经可以满足要求。

2.2 深度图对比

本节使用 4组数据进行深度计算, 并与文献[24]的基于微透镜的优化方法进行对比, 结果如图9~12所示。利用基于微透镜进行局部优化的方法得到的深度图虽然比较平滑, 但是存在错误的深度信

息。得益于MRFSP优化, 本文提出的方法能够增加深度图结果的鲁棒性, 且误匹配的深度信息被有效地抑制。图9(f)~12(f)中4个放大的图像中依次是WTA深度图、全局优化深度图、MRFSP优化的深度图以及利用基于微透镜半全局优化的深度图。

图9是对一个机械加工的钢结构件的三维信息

提取结果, 对比图9(b)~(e), 可以明显看出, 利用WTA 方法得到的深度图(图9(b))存在较多错误估计的深度信息, 图9(c)和(e)分别利用全局算法和基于微透镜的半全局算法, 有效地减少了深度误估计的区域。从图9(f)的区域放大图可知, 利用半全局优化的深度图存在误估计的区域, 与全局优化结果 (图9(c))相比效果较差。利用 MRFSP 融合全局优化以及 WTA 结果的深度图(图9(d)), 既保留了深度信息的连续性, 也能消除绝大部分的误匹配区域。

图10是对人脸进行三维信息提取的结果。可以看出, 图10(b)存在大量的噪点; 图10(c)依赖全局优化可以有效地去除噪点, 但有明显的阶梯感; 利用

基于微透镜的半全局优化算法得到的深度图(图10 (e))依然存在深度的误匹配点, 特别是对鼻子部位的深度估计。通过对比图10(f)的细节, 可以看出,利用 MRFSP 方法融合得到的深度图可以保持与全局优化算法一样的效果, 且阶梯感明显降低。

图11是对水果进行三维信息提取的结果。从图11(f)可以看出, 利用基于微透镜的半句全局优化与利用 WTA结果计算得到的深度图均包含明显的误估计点, 可知在此场景下, 基于微透镜的半全局优化不能有效地降低深度误估计, 而利用全局优化则有效地优化了误估计点。对比图11(f)的细节可知,利用 MRFSP 方法得到的深度图继承了全局优化的优势, 同时避免了全局优化中代价的离散性。

图12是对墙角进行三维信息提取的结果。由于该场景下墙角的纹理信息丰富, 基于微透镜的半全局优化(图12(e))和利用全局优化(图12(c))均能很好地避免深度误估计, 但是WTA方法得到的深度图依然有很强的噪声。对比图12(f)的细节可知, 利用MRFSP方法融合全局优化后, 误估计区域得到极大的抑制, 同时全局优化的离散性得到改善。

3 结论

本文探索了一种全局优化并获得连续深度图的方法。由于全光相机的结构与传统双目立体系统不同, 所以在进行深度计算时, 对初始的CV, 需用虚拟深度决定CV的高度离散量DL。为了提高全局优化的鲁棒性, 本文将代价信息投影到像方空间构建新的CV。此外, 为了解决全局优化只能获得离散深度图的问题, 在基于微透镜进行深度计算的同时保留WTA深度图, 并将该深度图同样投影到像方空间, 最终采用MRFSP优化获得既连续又鲁棒性好的深度图。

为了验证本文方法, 我们对不同代价离散量DL下获取的全局优化结果进行分析, 发现利用MRFSP进行优化的深度图可以用较小的DL 获得较好的连

[24]续性。与基于微透镜的半全局优化 方法对比发现, 本文提出的方法在钢结构件、人脸、水果和墙角等 4种常见场景下均能取得鲁棒性较好的深度结果。在场景纹理信息较弱的情况下, 基于微透镜的半全局优化方法获得的深度图包含大量的误估计点。值得注意的是, 全光相机目前适用于近距离摄影测量, 当距离增大时, 深度信息估计的精度会降低[27],因此本文选择的场景均属于近景摄影测量。

参考文献

[1] Yang P, Wang Z, Yan Y, et al. Close-range photogrammetry with light field camera: from disparity map to absolute distance. Appl Opt, 2016, 55(27): 7477–7486 [2] Georgiev T G, Lumsdaine A. Superresolution with plenoptic 2.0 cameras // Proceedings of the Signal Recovery and Synthesis, Optical Society of America. California, 2009: STUA6 [3] Ng R. Digital light field photography [D]. Palo Alto: Stanford University, 2006 [4] Dansereau D, Pizarro O, Williams S. Decoding, calibration and rectification for lenselet-based plenoptic cameras. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2013, 9(4): 1027–1034 [5] Bok Y, Jeon H G, Kweon I S. Geometric calibration of micro-lens-based light-field cameras using line features. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2017, 39(2): 287–300 [6] Bishop T E, Favaro P. Full-resolution depth map estimation from an aliased plenoptic light field // Proceedings of the 10th Asian conference on Computer vision: Volume Part II. Queenstown: 2010: 186–200 [7] Jeon H G, Park J, Choe G, et al. Accurate depth map estimation from a lenslet light field camera // Proceedings of the Computer Vision and Pattern Recognition (CVPR). Boston, 2015: 1547–1555 [8] Navarro J, Buades A. Reliable light field multiwindow disparity estimation // Proceedings of the 2016 IEEE International Conference on Image Processing (ICIP). Phoenix, 2016: 1449–1453 [9] Roberts W, Thurow B S. Correlation-based depth estimation with a plenoptic camera. AIAA Journal, 2017, 55(2): 1–11 [10] Yu Z, Guo X, Lin H, et al. Line assisted light field triangulation and stereo matching // Proceedings of the IEEE International Conference on Computer Vision. Darling Harbour, 2013: 2792–2799 [11] Criminisi A, Kang S B, Swaminathan R, et al. Extracting layers and analyzing their specular properties using epipolar-plane-image analysis. Computer Vision and Image Understanding, 2005, 97(1): 51–85 [12] Tao M, Hadap S, Malik J, et al. Depth from combining defocus and correspondence using light-

field cameras // Proceedings of the IEEE International Conference on Computer Vision. Darling Harbour, 2013: 673–680 [13] Wang T C, Efros A A, Ramamoorthi R. Depth estimation with occlusion modeling using light-field cameras. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2016, 38(1): 2170–2181 [14] Wanner S, Goldluecke B. Globally consistent depth labeling of 4D light fields. Computer Vision & Pattern Recognition, 2012, 23 (10): 41–48 [15] Zhang S, Sheng H, Li C, et al. Robust depth estimation for light field via spinning parallelogram operator. Computer Vision and Image Understanding, 2016, 145 (C): 148–159 [16] Zhang Y, Lv H, Liu Y, et al. Light-field depth estimation via epipolar plane image analysis and locally linear embedding. IEEE Transactions on Circuits & Systems for Video Technology, 2017, 27(4): 739–747 [17] Kolmogorov V, Zabih R. Multi-camera scene reconstruction via graph cuts // Proceedings of the 7th European Conference on Computer Vision: Part III. London, 2002: 82–96 [18] Kolmogorov V, Zabin R. What energy functions can be minimized via graph cuts?. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2004, 26(2): 147–59 [19] Felzenszwalb P F, Huttenlocher D P. Efficient belief propagation for early vision. International Journal of Computer Vision, 2006, 70(1): 41–54 [20] Yang Q, Wang L, Yang R. Stereo matching with colorweighted correlation, hierarchical belief propagation, and occlusion handling. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2009, 31(3): 492– 504 [21] Georgiev T, Lumsdaine A. Focused plenoptic camera and rendering. Journal of Electronic Imaging, 2010, 19(2): 021106 [22] Zeller N, Quint F, Stilla U. Calibration and accuracy analysis of a focused plenoptic camera. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2014, 2(3): 205–212 [23] Zeller N, Quint F, Stilla U. Depth estimation and camera calibration of a focused plenoptic camera for visual odometry. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 118: 83–100 [24] Fleischmann O, Koch R. Lens-based depth estimation for multi-focus plenoptic cameras // Proceedings of the German Conference on Pattern Recognition. Münster, 2014: 410–420 [25] Hirschm H. Accurate and efficient stereo processing by semi-global matching and mutual information. Proceedings of the Computer Vision and Pattern Recognition, 2005, 2(2): 807–814 [26] Perwass C, Wietzke L. Single lens 3D-camera with extended depth-of-field // Human Vision and Electronic Imaging XVII, International Society for Optics and Photonics. San Diego, 2012: 829108 [27] Sardemann H, Maas H G. On the accuracy potential of focused plenoptic camera range determination in long distance operation. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 114: 1–9 [28] Strobl K H, Lingenauber M. Stepwise calibration of focused plenoptic cameras. Computer Vision and Image Understanding, 2016, 145(C): 140–147 [29] Janoch A, Karayev S, Jia Y, et al. A category-level 3D object dataset: putting the kinect to work // Proceedings of the IEEE International Conference on Computer Vision. Washington, DC, 2011: 1168–1174 [30] Kanade T, Okutomi M. A stereo matching algorithm with an adaptive window: theory and experiment. IEEE Transactions on Pattern Analysis & Machine Intelligence, 1994, 16(9): 920–932

根据相态, 大气中的云可以大致分为 3 种类型:水云(液态)、冰云(固态)以及混合云(混合相态)。温度是云相态判识的重要依据: 云顶温度在 0℃以上的云为水云, 云底温度在−40℃以下的云为冰云, 0~−40℃之间的云可能为水云、冰云或混合云[1]。0℃以下的水云称为过冷水云, 云粒子因为缺乏冰核或其他原因仍然呈现为液态, 状态不稳定。云相态识别是精确评估云辐射强迫的前提, 也是遥感反演云的其他微物理特性参数的必要条件。对过冷水云区的识别, 在提高人工增雨的成功率以及保障航

空飞行器飞行安全等方面具有重要意义。

目前对过冷水云的探测, 主要有机载粒子测量系统实地探测、探空资料反演以及微波辐射计、激光雷达、毫米波雷达等主动和被动传感器探测方法[2]。Hogan 等[3]利用航天飞机携带的激光雷达估计层状云中的过冷水在全球的分布, 发现在−10 ~ −15℃之间的云中有 20%左右含过冷水。Hu 等[4]依据 CALIOP/IIR/MODIS 测量结果, 发现在高纬度地区 0~−40℃之间的低层云中, 过冷水云的比例达到

[5] 95%。Westbrook 等 根据激光雷达及毫米波雷达

探测资料, 利用阈值法识别云顶过冷水层。Hogan

[67]等 利用地基激光雷达, 并结合毫米波雷达得到的光学厚度来识别过冷水存在区域, 发现在−10~ −15℃之间有 20%的云中含过冷水。Shupe 等[89]利用多种仪器的探测资料, 对混合云的特性进行研究, 证明液态云的比例随着温度的升高而增大, 总结出识别过冷水的阈值。Luke 等[10]利用毫米波雷达多普勒速度谱的特征来识别过冷水, 结果与云高仪测得的结果一致。国内对过冷水的研究一般采用机载粒子测量系统或结合星载雷达, 如严卫等[11]根据 Cloudsat 携带的激光雷达和 94 GHZ 测云雷达资料, 利用支持向量机法识别出混合相云。但是, 机载仪器探测受时空及飞机飞行的限制, 当星载雷达与云层之间距离过远时, 信号容易受到干扰。在上述对过冷水云的研究中, 有的所用仪器众多, 方法复杂, 有的只研究云顶的过冷水, 有的因时间受限而不能作为长期观测手段。利用激光雷达探测资料进行云检测, 在星载激

[12] [13]光雷达 、机载双波长偏振激光雷达 和地基激光雷达[14]等方面均有相关研究, 检测方法各有优缺点, 但都采用阈值筛选法, 无论一个阈值还是两个阈值, 检测结果都存在不确定性。随着激光雷达的发展, 利用激光雷达线性体积退偏振比(linear volume depolarization ratio, 简称退偏比)的观测来识别云的相态, 具有时间分辨率高, 探测手段简单易行的优势[1516]。本文利用地基微脉冲激光雷达识别过冷水云层, 技术上主要分两步。首先利用激光雷达后向散射信号区分云层和气溶胶层; 然后对检测出来的云层进行相态识别, 再结合探空数据的温度廓线, 鉴别出云层中的过冷水。本文对前人的云识别分类算法做了改进, 并利用改进后的算法对2016 年 3 月20日至2017年3月 19日北京大学观测点上空出现的云进行识别、分类和统计分析, 取得令人信服的结果。

1 设备和原理

微脉冲激光雷达(Micro Pulse Lidar, MPL)是 20世纪 90 年代美国国家宇航局(NASA)研制的一种人眼安全的探测云和气溶胶的新型激光雷达, 具有激光脉冲能量稳定、重复频率高、信噪比高等优势。本研究采用 Sigmaspace 公司生产的 MPL, 观测点位于北京大学物理楼六层楼顶。该激光雷达的光源

为 ND:YV04 532 nm 激光, 脉冲宽度为10.3 ns, 脉冲能量为 6~8 μj, 脉冲频率为 2500 Hz。该激光雷达出射线偏振脉冲信号, 同时以光子计数的方式接收与出射激光偏振方向平行和垂直的弹性散射信号,最大探测高度夜间约为 45 km, 白天无云、清洁条件下可达 15 km, 垂直分辨率为 15 m。该激光雷达接收视场角较小, 接收的后向散射信号以单次散射为主, 因而降低了处理多次散射问题的复杂性, 并在一定程度上减小了太阳光背景噪音的影响。激光雷达接收到的信号强度可以用式(1)表示: p ( z)

C  E  O ( z )  ( z)  exp  2  z  ( z  )d z   N  A ( z) z2  0  b

, (1)

D  p ( z) 

式中, p ( z)是探测器接收到的信号强度, C是系统常数, E 是脉冲能量, O ( z)是重叠区订正因子,  ( z )是后向散射系数, z 是探测高度,  (z) 是消光系数, Nb 是背景噪音, A ( z) 是残余脉冲, D [ p ( z )]是延时订正函数。对原始信号 p(z)进行探测器延时订正、残余脉冲订正、背景噪声订正、距离订正和重叠区订正, 并对发射能量进行归一化处理, 可以得到激光雷达归一化后向散射(normalized relative backscatter)信号NRB(Z): [ p ( z )  D [ p ( z )] N  A ( z )]  z2 NRB(Z)  b 。 (2) E  O ( z)

激光雷达发出的是线偏振光, 遇到各种粒子发生后向散射后, 偏振状态会发生不同程度的改变。激光雷达退偏比 可由下式求得: Cp NRBS( z)

 ( z)  , (3) Cs NRBP( z)其中, 下角标 s 表示与发射光偏振方向垂直的探测通道, 下角标 p表示与发射光偏振方向平行的探测通道。

利用激光雷达退偏比, 可以识别云的相态。原理是: 云层中的水滴都是球形或近球形的轴对称粒子, 激光雷达发射的线偏振光经过球形粒子和轴对称粒子散射后, 返回的后向散射光偏振状态不发生改变, 即 NRBS( z) 为零, 因此水云的退偏比很小;冰晶粒子多为不规则的形状, 返回的后向散射光的偏振状态会发生显著的改变, 即 NRBS( z) 会增加,

NRBP( z) 会减小, 因此冰云的退偏比较大。

2 云识别算法

背景噪音、仪器观测的随机误差以及低层高浓度气溶胶会导致高层信号信噪比低, 使得很多利用激光雷达探测信号进行云识别的算法都存在错判和误判的情况。Zhao 等[17]提出一种云和气溶胶识别算法, 主要分为 3 个步骤: 1) 由于对雷达信号进行距离订正会放大高处的噪音信号, 因此在进行距离订正前, 先对雷达信号进行三点滑动平均以及半离散化处理, 消除随机噪音; 2) 用直方图均衡化方法,判断云层和气溶胶层; 3) 选择合适的阈值, 利用定

[17]义的特征函数区分云和气溶胶。我们在 Zhao 等算法的基础上, 优化了区分云与气溶胶的方法。下面详细介绍改进后的云识别算法。

2.1 去除噪音

未进行距离订正的激光雷达后向散射信号可由式(4)得到:

p ( z )  D [ p ( z )]  N  A ( z) P ( z)  b 。 (4) C  E  O ( z)

P ( z)虽然去除了背景噪音, 但是信号仍然受测量的随机噪音和大气湍流影响。

计算 15 km 以上 P ( z) 的标准差, 将 5 作为随机噪音的阈值。首先对 P(z)做三点滑动平均, 得到 Ps(z)。然后进行半离散化处理: Ps 按照高度递增的顺序, 即 Ps(zi)(i=2, 3, 4, …, N), 如果 Ps(zi)与Ps(zi−1)之差的绝对值小于随机噪音阈值5 , 则Ps(zi)=ps(zi−1), 否则 Ps(zi)保持不变, 将新获得的Ps(z)记为 PD1(Z); Ps 按照高度递减的顺序, 即 Ps(zi) (i=n−1, N−2, N−3, …, 1), 如果 Ps(zi)与 Ps(zi+1)之差的绝对值小于随机噪音阈值5 , 则 Ps(zi)=ps(zi+1),否则 Ps(zi)保持不变, 新获得的 Ps(z)记为 PD2(Z)。将 PD1(Z)和 PD2(Z)平均, 得到 PD(Z)。我们认为,当信号之间的差异大于噪音水平时, 信号才是有效的。这样, 随机噪音和湍流对信号的影响已经被减到最小, 可以使用 PD(Z)作为信号进行云层和气溶胶层的判断。

2.2 判断云层和气溶胶层

使用直方图均衡化方法处理信号 PD(Z), 判断信号中的云层和气溶胶层。假设每一组信号有N个数据, 则 PD(Z)可表示为 Pd(i)(i=1, 2, 3, …, N)。

1) 将 PD(Z)进行升序排列, 得到 Rs(i)(i=1, 2, 3, …, N)。排列后, Rs(i)在原序列 PD中的位置记为Is(i), 即 Rs(i)=pd(is(i))。最大值 RS(N)记为 MA,最小值 Rs(1)记为 MI。

2) 定义 Pe(i)=i/n(i=1, 2, 3, …, N), 如果 Rs(i)= Rs(i−1), 则 PE(I)= PE(I−1)。

3) 重新建立一个升序序列 Y(I)=PE(I)·(MA− MI)+MI (i=1, 2, 3, …, N)。再将 y(i)按照高度排序,得到直方图均衡化后的信号 PN(Z), 即 PN(IS(I))= y(i)(i=1, 2, 3, …, N)。如果整个垂直方向没有云层或气溶胶层出现, 那么 PN应该随高度 z的增加而减小。

4) 定义基准线 B(z), 两个端点分别为为(z1, MA)和(zn, MI)。理论上, 如果信号中不包含云层或气溶胶层, 信号 P(z)应该随高度的增加而减小, PN(Z)与 B(z)应该相等。如果信号中包含云层或气溶胶层, 就会出现高处的信号大于低处信号的情况, 那么就可以通过 PN(Z)与 B(z)的差异区分出云层和气溶胶层。按照高度增加的顺序, 寻找 PN(Z)开始大于 B(z)的高度作为一层的起点, 记为 base(ilayer),将之后 PN(Z)小于 B(z)的第一个点作为这一层的终点, 记为 top(ilayer), 其中 ilayer 是垂直高度上分出来的层数。同时, 为了降低随机信号的干扰, 要求每一层的厚度要大于 45 m。这样, 就得到一组信号中云层和气溶胶层的位置。

2.3 区分云和气溶胶

判断出云层和气溶胶层后, 需要区分云与气溶胶。一般情况下, 云层比气溶胶层出现在更高的位置, 并且云粒子的后向散射比气溶胶层强。这个特征应该体现在云层内的 PN 与基线 B 的差异更大。定义 Area(ilayer)为云层内 PN 与 B 围成的归一化面积:

 top(ilayer)

[PN( z )  B( z )]dz Area(ilayer)  base(ilayer) (MA  Mi)[top(ilayer)  base(ilayer)]。(5)

考虑到底层大气污染的严重程度也会影响阈值的选取, 我们用归一化后向散射信号的积分FMX来表征底层大气的污染情况:

 base(ilayer)

Fmx(ilayer)  NRB( z )d z。 (6) 0

最终, 确定新的阈值函数F:

F (ilayer)  Area(ilayer)  Fmx(ilayer) 。(7)经过比对大量数据中阈值函数的输出结果, 考虑到不同污染情况及不同高度的差别, 结合人工判断, 我们发现, 可以选取 1000 作为统一的阈值来区分云与气溶胶, 即当 F>1000 时为云层, 否则为气溶胶层。以 2016 年 3 月 20 日为例, 24 小时的激光雷达归一化后向散射信号如图1所示。对 2016年 3 月20日的雷达信号进行处理, 云识别结果如图2 所示。可以看出, 改进后的算法避免了高空噪音信号的误判以及近地面密度较大气溶胶层的误判,效果很好, 与人工在图像上对云的识别结果基本上一致。

3 过冷水云的识别

识别过冷水云的难点主要是云相态的识别。利用退偏比来区分云的相态, 必须要求后向散射信号以单次散射为主。本文采用的激光雷达接收视场较小, 满足收到的信号以单次散射为主, 可以忽略水云多次散射带来的偏振现象。即使考虑了接收到部分多次散射信号, 水云底部退偏比接近 0, 并随厚度增加呈线性增长, 在云顶也不会超过 0.1。下面根据北京市南郊观象台探空数据得到的温度廓线,找到典型的水云和冰云, 分析其退偏比特征。

图 3 为 2016 年 5 月 14 日 05:00 的温度廓线、退偏比及归一化后向散射信号廓线, 1.6 km 左右的云层所在高度的温度为零上, 可以看出是典型的水

云。在云层中观察退偏比的特征, 发现从云底到云顶, 退偏比经历了从大到小、趋近 0、再从小变大的过程。分析其原因, 云底是气溶胶与云接触的位置, 因为气溶胶的存在, 退偏比稍大。在达到最小值后, 随着云层厚度的增加, 激光雷达接收到的信号中有部分来自水云的多次散射, 因此退偏比随高度呈线性地增加, 但也小于 0.1, 一般在 0.05 左右。

再来看一下冰云退偏比的特征。从图 4 中温度廓线看, 7~9 km 的云层对应的温度都在−30℃以下,是典型的冰云, 可以看出在云层所在位置, 退偏比有显著的变化, 普遍大于 0.3。

冰云中有一种特殊的云, 云中的冰晶是固定取向的水平片状冰晶。水平导向片状冰晶的后向散射信号很强, 而且偏振状态不发生改变, 情况与水云

相似。与水云不同, 水平冰晶的退偏比虽然很小,但没有趋近零的趋势, 是不均匀的, 也非单调的,不会随着云层厚度的增加而呈线性变化(图 5)。

因此, 结合退偏比以及温度廓线, 可以判断云的类型。温度在 0℃以下, 退偏比大于 0.3 的是颗粒状冰云; 温度在 0℃以下, 退偏比在 0.05~0.3 之间的是混合相态云; 退偏比小于 0.05, 从云底到云顶, 退偏比先减小后线性增加的是水云, 水云中温度在 0 ℃以下的则为过冷水云; 温度在 0 ℃以下,退偏比小于 0.05, 从云底往上, 退偏比没有线性增长趋势的是水平导向的片状冰晶云。

以 2016 年 5 月 4日的激光雷达数据为例展示分类结果。结合激光雷达的归一化后向散射信号(图 6)、退偏比(图 7)和云分类结果(图 8)可以看出, 10:30 开始出现后向散射信号比较强的厚云层, 退偏比数值较小, 分析其变化特征后, 确定为水平导向的片状冰晶。13:30 左右退偏比发生变化, 其间夹杂过冷水、混和云及冰云, 14:00— 15:00 的云层主要是冰云。从 15:00 到 20:30 左右, 云层不稳定,各种云都有, 18:30 左右有短暂的降水。从 20:30 直到 24:00, 有一层后向散射信号较强的云层, 从退偏比变化特征来看, 主要是过冷水云层。

结合退偏比信号特征和温度廓线, 可以将已经得到的云检测图划分为五类: 0℃以上的水云、随机导向颗粒状的冰云、混合云、水平导向的片状冰晶 以及过冷水云(图 8)。只有在比较复杂的天气情况下, 五类云才会出现在一张云图里, 出现各类云夹杂的情况。当然, 也可能存在误判的情况, 主要是水平导向的片状冰晶与过冷水云相似度太高, 分类算法不够完善所致。

4 统计结果

北京大学物理楼六层楼顶的微脉冲激光雷达于2016 年 3月 20日调试完毕, 正式启用。本文统计了 2016 年3月20日至 2017 年3月 19日为期一年的数据, 有效数据总时间为344 天5小时12分钟。其中, 有云的时间为151 天 12小时6分钟; 有0℃以上水云的时间为41 天 2小时 6分钟, 平均高度为2.72 km; 有随机导向颗粒状的冰云的时间为 89 天21 小时 24 分钟, 平均高度为 8.03 km; 有混合云的时间为 116 天 12 小时 24 分钟, 平均高度为 6.52 km; 有水平导向的片状冰晶云(水平冰)的时间为 21天 6 小时, 平均高度为 4.62 km; 有过冷水云的时间为 14 天 21 小时 48 分钟, 平均高度为 4.90 km。过冷水云出现的时间占有云时间的比例为 9.84%; 对于全部 0~−40℃的云, 过冷水云出现的时间占比为11.99%。各种类型云层的月平均高度变化趋势如图 9 所示。3 月的数据是 2016 年 3 月下旬的数据与 2017年 3 月上中旬的数据拼到一起的, 因此在变化趋势

上不具代表性。结果表明, 水云的平均高度是所有类型云层中最低的。冰云的月平均高度最高, 然后是混合云, 水平冰和过冷水云的月平均高度位于中间且不分上下。各云层平均高度的年变化趋势都与温度正相关。

各种类型云层出现时间的占比情况如图 10 所示。可以看出, 1 月、2 月和 12 月没有出现水云,其他月份水云占有效数据时间的比例和水云在云中所占比例都与温度正相关; 冰云和混合云在云中所占比例与温度负相关; 水平冰和过冷水没有明显的变化规律。各月份过冷水云在 0~−40℃的云中的占比如图 11 所示, 过冷水云的占比在 8 月、9 月和 1月出现峰值。

5 结语

本文利用地基微脉冲激光雷达, 对过冷水云层进行识别, 方法简单易行, 时间分辨率高, 具有其他方法无可比拟的优势。对已有的云识别算法进行改进, 不仅避免了使用分段阈值的不便, 而且结果更准确。

利用长达一年的资料进行分析, 发现观测站点上空过冷水云出现的时间占有云时间的 9.84%, 对于全部 0~−40℃的云, 过冷水云出现的时间占比为11.99%, 过冷水云在 8 月、9 月和 1月出现峰值。

在对云的相态进行判别时, 主要依靠不同类型的云滴粒子在退偏比方面的差异。水云、冰云和混

合云都很容易区分, 难点在于区分过冷水云与水平导向的片状冰晶。这不仅因为两者的退偏比及激光雷达后向散射信号都极其相似, 而且因为两者之间也存在相互转化和夹杂的情况。因此, 区分过冷水和水平导向的片状冰晶的算法难以做到尽善尽美。

另外, 由于条件限制, 本文采用的温度廓线不是实时获得的, 而是用北京市南郊观象台早晚的探空数据进行平均和插值得到, 因此针对整天进行判断会存在一定的误差。如果条件许可, 今后可以采用飞机探测的云微物理观测结果, 与地基激光雷达进行云相态判别的比对, 也期待与云雷达的观测数据进行比对。

参考文献

[1] Pruppacher H. A new look at homegeneous ice nucleation

1924‒1933 in supercooled water drops. Journal of Atmosphere Science, 1995, 52(11): [2] Wu Juxiu, Wei Ming, Wang Yilin. Retrieval of the supercooled water in stratiform clouds based on

227‒235 Milimeter-wave Cloud Radar. Journal of Arid Meteorology, 2015, 33(2): [3] Hogan R J, Behera M D, O’connor E J, et al. Estimate of the global distribution of stratiform supercooled

325‒341 liquid water clouds using the LITE lidar. Geophys Res Lett, 2004, 31(5): [4] Hu Y, Rodier S, Xu K, et al. Occurrence, liquid water content, and fraction of supercooled water clouds from combined CALLOP/IIR/MODIS measurements. J Geophys Res, 2010, 115: D00H34 [5] Westbrook C, Illingworth A. Evidence that ice forms primarily in supercooled liquid clouds at temperatures >−27℃. Geophys Res Lett, 2011, 38: L14808 [6] Hogan R J, Francis P N, Flentje H, et al. Characteristics of mixed-phase cloud: 1. Lidar, radar and

2089‒2116 aircraft observations from CLARE’98. Q J R Meteorol Soc , 2003, 129: [7] Hogan R J, Illingworth A J, O’connor E J, et al. Characteristics of mixed-phase cloud: 2. A climatelogy

2117‒2134 from ground-based lidar. Q J R Meteorol Soc, 2003, 129: [8] Shupe M D, Matrosov S Y, Uttal T. Arctic mixedphase cloud properties derived from surface-based

697‒711 sensors at SHEBA. Journal of the Atmospheric Sciences, 2006, 63(2): [9] Shupe M D. A ground-based multiple remote-sensor cloud phase classifier. Geophys Res Lett, 2007, 34: L22809 [10] Luke E P, Kollias P, Shupe M D. Detection of supercooled liquid in mix-phase clouds using radar Doppler spectra. Journal of Geophysical Research, 2010, 115: D19 [11] 严卫, 任建奇, 陆文, 等. 联合星载毫米波雷达和

68‒73激光雷达资料的云相态识别技术. 红外与毫米波学报, 2011, 130(1): [12] 卜令兵, 单坤玲, 吕晶晶, 等. 云和气溶胶探测星62‒64载激光雷达及其应用. 光学与光电技术, 2009, 7(4): [13] 伯广宇, 刘东, 王邦新, 等. 探测云和气溶胶的机203‒208载双波长偏振激光雷达. 中国激光, 2012, 39(10): [14] 邱金桓, 郑斯平, 黄其荣, 等. 北京地区对流层中

1‒7上部云和气溶胶的激光雷达探测. 大气科学, 2003, 27(1): [15] Hu Y, Winker D, Vaughan M, et al. CALIPSO/

2293‒2309 CALIOP cloud phase discrimination algorithm. Atmos Ocean Tech, 2009, 26(11): [16] Schotland R M, Sassen K, Stone R. Observations by

1011‒ lidar of linear depolarization ratios for hydrometers. Journal of Applied Meteorology, 1971, 10(5): 1017 [17] Zhao Chuanfeng, Wang Yuzhao, Wang Qianqian, et al. A new cloud and aerosol layer detection method based

614‒616 on micropulse lidar measurements. Journal of Geophysical Research, 2014, 119(11):

中国夏季降水受到东亚季风环流系统的影响,表现出区域差异以及年际和年代际变化的特征[1–2]。中国东部夏季降水到底有多少种空间上的异常分布类型, 需要利用客观的方法加以识别。前人对中国东部夏季降水异常分布提出不同的分型或分区。廖

[3]荃荪等 根据夏季总降水量距平百分率及500 hpa环流形式, 将中国东部地区夏季大范围降水异常分布概括为主要多雨带分别位于黄河流域及其以北、黄河至长江之间以及长江流域或江南3个类型。魏

凤英等[4]用经验正交函数(empirical or-thogonal function, EOF)中前3个特征向量的时间系数, 将中国东

[5]部夏季雨带划分为3种主要类型。王绍武等 利用

[6] EOF, 划分6种雨型。孙林海等 利用EOF和主成分分析等方法, 将中国东部季风区降水划分为二类四型。黄荣辉等[7]揭示出中国东部降水从南到北的经向三极子型和偶极子型两种分布模态, 并发现存在准两年的周期振荡和年代际变化。龚振淞等[8]将阻塞高压和西太平洋副高的配置与中国夏季旱涝分

布特点结合起来, 将中国夏季旱涝分布划分为 8 种雨带类型。可见, 前人的分型(或雨带)从 3 个到 8个不等。

前人对中国夏季降水的客观分型多以 EOF 方法为主, 同时考虑大气环流异常[3,5,8–13]。由于 EOF分解需要满足模态之间的正交性, 其结果不能完全反映中国夏季降水与大气环流和外强迫异常之间的关系。正如Dommenget等[14]指出的, EOF分解得到的模态缺少物理含义, 即模态可能并不存在。作为一种聚类和模式识别方法, 自组织映射(self-organizing maps, SOM)方法在 20 世纪 90 年代后期引入气象学和气候学研究中[15–19]。人们发现, SOM方法在不同方向和不同时空尺度的研究中都是一个有效的工具, 例如极端气候与降水模态分析、云型分类、气候变化分析等方面的研究[20]。Tadross 等[21]利用 SOM 方法, 提取南非和津巴布韦

[22]降水的特征模态。Nishiyama 等 使用SOM方法,分析降水和850 hpa风场的联合分布。Lin等[23]通过比较, 发现SOM方法的聚类结果优于K-means 和

[24] Ward’s方法。Hsu等 使用SOM和小波分析法, 分析了台湾地区22年降水数据的时空变化特征。

前人的研究表明, 与 EOF 方法相比, SOM 方法具有更多的优势。由于模态之间有内在的联系, 将SOM分型结果以矩阵的形式排列, 是可视化方面

[25–26]的一个突破。Liu 等 发现, SOM 分型往往比

[27] EOF的前几个模态更准确和直观。Liu 等 发现,在较复杂的模态集合中, SOM方法能成功地提取所有模态, EOF方法却不能。Reusch等[28]发现, EOF方法有时会将多个模态混合成单个模态, 或者不能正确地分解不同部分的方差, 而 SOM 方法能够鲁棒性更强地提取不同模态及对应的方差。有研究认为, 当取特定参数时, SOM 方法即退化为 K-means方法[29], 并且更灵活[30]。

本文基于 EOF 方法在中国东部夏季降水分型中的广泛应用以及 SOM 方法在其他地区降水分型中的应用, 使用这两种方法分析中国东部夏季降水异常, 藉此比较 EOF 和 SOM 方法的差异, 并考察SOM分型的天气和气候环流异常背景。

1资料与方法

1.1资料

中国东部地区降水数据来自中国气象局提供的1961—2010 年中国日降水格点数据集, 空间分辨率

为 0.5°×0.5°。大气风和水汽资料来自NCEP Reanalysis 1, 空间分辨率为 2.5°×2.5°。位势高度场资料来自 ERA Interim, 空间分辨率为 0.75°×0.75°。我们用总场变量减去对应日(30 年)的瞬时气候场, 得到逐日扰动场[31]。

1.2 EOF 分析

EOF 分析, 也称为主成分分析(principal component analysis, PCA), 最早由 Person[32]提出。20 世纪 50 年代, Lorenz[33]率先尝试将 EOF方法应用于天气预报。通过EOF分析, 可以得到气象要素场变化的正交模态空间分布特征及其系数随时间的演变。它的前几个模态往往能够解释较多的时空变化方差, 因此该方法常用于气象要素场时空变化特征的研究。

EOF方法的基本原理是将由m个空间点n次观测值(样本时间长度)构成的变量矩阵X 分解为m n空间特征向量矩阵V以及对应的时间系数矩阵T两部分:  v v  v  t t  t1m   11 12 1m  11 12  v v  v t t  t2m  21 22 2m  21 22  X  。 (1)  m n              vm1 vm2  vmm  tm1 tm2  tmm 

对于第 i 个站点第 j 次观测值 xij, EOF 可将其分解为空间函数和时间函数两部分。其中, 特征向量表示某区域要素场的变率分布结构, 可以代表该要素场的主要分布特征; 所对应的时间系数则反映实际分布与之相似或相反, 且绝对值越大, 表示该空间分布特征越显著。

EOF的优势在于可以反映时空变率最大的模态, 但它的劣势也很明显: 1) 分解得到的模态不稳定, 受选取的时间序列长度和空间区域大小影响; 2) 选取的 EOF 前几个模态通常无法解释所有方差; 3) EOF 分解要求各个模态是正交的, 这可能导致分解得到的模态没有物理意义。

1.3 SOM 方法

SOM 方法是一种基于竞争学习的无监督式的神经网络方法[34–36],能够将高维空间的输入样本非线性地映射到一组二维格点(神经网络节点)中。每一个节点都有一个权重向量 mi, 它可以被随机地初始化, i 的取值从 1 到 M, 其中 M为 SOM 矩阵的大小。具体的训练方法为: 每次都选取其中一个输入样本(向量) Xk, 计算其到每个神经元节点的激活函

数。通常将欧氏距离作为激活函数, 激活函数最小(与被选取的输入样本距离最小)的神经元将被激活。该过程可被表示为

ck = argmin||xk − mi||, (2)其中, ck即为激活神经元的序号, argmin||xk − mi||表示||Xk − mi||最小时 i 的取值。我们可将激活神经元的权重向量移向被选取的输入样本。同样地, 激活神经元周围的邻近节点的权重向量也被调整并向输入样本移动。与此同时,移动的距离与临近函数有关。该过程可被表示为

mi(t +1) = mi(t) + α(t)hci(t)[x(t)− mi(t)], (3)其中, t 表示当前的迭代次数, X 表示当前选取的输入样本, α(t)为随时间衰减的学习参数, hci(t)为邻近函数。

不断重复上述迭代过程, 直到收敛, 即为自组织映射算法。SOM 节点的输出权重向量被改造为与输入数据有相同的特征。这种学习过程导致对输入数据有次序的映射, 相似的模态映射到临近区域,而差异较大的模态则映射到相距较远的区域。

与传统方法相比, SOM 方法主要有以下优势: 1) SOM 分型结果为有内在联系的连续态, 而非离散的模态, 这能够帮助我们理解关键模态之间的过渡态, 并对关键模态之间转换的预报提供帮助; 2) SOM 方法是可视化方面的一个突破, 通过将 SOM分型结果以矩阵的形式排列, 使我们更容易理解这些分类结果。

本文使用降水距平资料进行 SOM 分型, 没有 进行归一化处理, 选取Ep 函数(Epanechnikov function)作为邻近函数(neighborhood function), 即

hci(t) = max{0, 1 − (σt − dci)2}, (4)其中, σt 为 t 时刻的邻近半径, dci为激活神经元到邻近神经元的距离。选取最大训练步长为 100 步。

与所有动态聚类问题一样, SOM方法需要事先指定聚类数目。然而, 最优聚类数目难以客观选取,并且尚无定论[37–39]。在此, 可以用 Chattopadhyay[40]的链接网络概念来确定分型数目范围。在一个长度为(L+1)的时间序列中, 共有 L 种演变形态, 而每一种演变形态又对应 SOM 中一个节点到另一个节点的转换。那么, 有

N(N − 1)/2 = L, (5) N即为神经元节点个数, L为演变形态数目。

对季节平均降水, 我们关注的是降水的年际及年代际变化, 因此选取 1961—2010 年共 50 年的数据, L 可取为 50, 解得 N ≈ 10。对日平均降水, 我们关注的是夏季降水的日变化。由于夏季的时间尺度为 3~4 个月(90~120 天), 即L = 90~120, 解得 N= 14~16。

2 中国东部夏季季节平均降水距平EOF 与 SOM模态的比较 2.1 空间形态分析

我们先用夏季 3 个月(6—8 月)的季节降水距平资料做 EOF 分析。图 1 是我国东部夏季降水距平EOF 分解得到的前 3 个模态的空间分布, 它们的解

释方差分别为 16.8%, 15.9%和 9.4%。EOF 第一模态最大的异常中心在华南, 并表现为以长江为界的南北反向振荡(图 1(a)); 第二模态的异常带沿长江分布(图 1(b)); 第三模态表现为以环渤海为中心的北方异常以及东南地区的反向变化(图 1(c))。上述EOF 特征与前人的分析结果[13,41]相似。前人的研究表明, 中国东部季节尺度的夏季降

[4–5]水主要有3 个 EOF模态 。据图 1, 用6 个 SOM分型已经足以解释前3个EOF模态, 因此选取6个分型是比较合适的。

中国东部夏季降水距平SOM分型如图2所示。这 6 个(3×2)分布类型可以分为两两相反的 3组。A组包括1型和6型, 这组降水型主要反映以长江为界, 南、北方降水反向变化的特征, 但三北地区(西北、华北和东北)又与南方相同, 表现为南方偏旱(涝)时, 黄淮流域降水偏多(少), 三北地区偏旱(涝)。B组包括2型和5型, 反映西北与其他地区的相反变化。C组包括3型和4型, 与A组有相似之处, 但异常中心的位置发生变化。SOM的 6个分

型中, 中心强度和分布范围是连续变化的, 有其内在的规律性。从1→4→5→6型看, 华南的干旱中心逐步北移到淮河; 从6→3→2→1型看, 南方的湿中心在北移。SOM的6个分型中, 干(湿)中心的强度和位置是连续变化的, 这在 EOF分析结果中是看不到的。在50年的夏季降水中, 有24%的年份降水分布与4型相似, 即长江及其以南降水偏少。各有18%的年份, 江淮降水偏多(少)而华南降水偏少(多), 类似1型和6型。只有10%的年份, 江淮夏季降水偏少(5 型)。

进一步对比 EOF 分解结果和 SOM 分型可见,1型对应 EOF1 的正相位,6型对应 EOF1 的负相位;3型对应 EOF2 的正相位,4 型和 5 型对应 EOF2 的负相位; EOF-3 在 SOM分型中不明显, 即 EOF的第三模态不是独立存在的。

2.2 时间序列分析

图3 是3个 EOF模态对应的主分量的时间序列。正如黄荣辉等[7]指出的, 中国东部夏季降水EOF的前3个模态均呈现一定程度的周期性震荡。

EOF第一模态表现出明显的年代变化, 如 1980— 1993年期间以正值为主, 其前后则以负值为主, 存在年代际变化; EOF第二模态以几年尺度和年际变化为主; EOF第三模态也存在年代际和年际变化。但是, EOF主分量序列的时间变化特征不够直观。

对每年夏季的降水空间分布, 我们可以通过计算其与图2 所示 6 个 SOM分型的欧氏距离, 得到最佳匹配神经元(best matching unit, BMU), 即与之 最相似的降水分型。图 4 给出 1961—2010 年夏季BMU 出现的时间序列, 其中, 1, 4, 5 型为南方偏旱型, 用实心符号表示; 2, 3, 6 型为南方偏湿型, 用空心符号表示。A 组(1, 6 型)用圆形符号表示,B组(2, 5 型)用三角形符号表示, C 组(3, 4 型)用正方形符号表示。1, 6 型是华南地区两个相反的干湿型,对应的降水偏多(偏少)信号较强。从图4可以清楚地看到, 20 世纪 80 年代至 90 年代初期,1 型居多

(华南华北干, 降水偏少), 而 6 型在 90 年代中后期

[1]偏多(华南湿, 降水偏多), 这与黄荣辉等 的结论一致。在 1995 年前后, 3型连续出现, 长江流域降水偏多。此外, 在图4中也展现出不同区域的年际变化。 2.3 合成的环流和水汽扰动根据图4中各个夏季降水异常分型对应的年份, 分别合成相应的850 hpa 扰动风和700 hpa 水汽扰动。比较图2与图5可以看出, 这些降水分型的空间特征与 850 hpa 扰动风的辐合以及 700 hpa 水汽扰动中心有密切的联系。例如, 华南沿海降水偏多的 6 型是华南地区850 hpa扰动风的辐合与700 hpa水汽扰动中心共同作用的结果。同样, 其他降水型的干湿分布都可用850 hpa扰动风和 700 hpa

水汽扰动来解释。

3 中国东部夏季日平均降水距平EOF与 SOM模态的比较3.1 空间形态分析

图 6 给出夏季日平均降水距平 EOF 分解得到的前 6 个模态的空间分布。EOF 第一模态反映长江南北相反的逐日降水异常分布, 第二模态反映沿长江降水与华南降水的相反分布, 第三模态反映长江下游与黄河中下游及华南的相反降水分布, 第四模态反映环渤海和东北与黄河中游的相反降水分布, 第五模态反映黄河和江南降水与淮河和华南的相反降水分布, 第六模态反映南方东、西部相反的降水分布。解释方差从第一模态的 9.1%逐渐下降到 3.5%。前 6个模态的总方差贡献为 34.4%。

对比图1与图2可以看出, 对季节平均降水来说, 前3个 EOF模态相当于 SOM展开的 6个型。对日平均降水来说, EOF的第五和第六模态中的正负中心最多达到5个。逐日降水SOM分型的 16个型主要反映EOF 前 3个模态的变化。在图7的 1 型中, 大降水带在长江流域与淮河流域之间, 该大降水带以南和以北降水偏少。降水带的地理位置变化是有规律的, 按图 7 中灰色虚线箭头方向, 降水偏多的中心带从华北向华南沿海靠近。

为了得到各个降水型之间的转变规律, 我们对30 年中的所有天数做统计分析, 得到各个降水型后一天出现的降水型概率的估计(图略)。比如, 在出现 1 型的后一天, 出现 2 型的概率是 23.1%。但是,出现14型的第2天再次出现14型的概率是 33.2%。根据每个降水型后一天的最大转移概率, 我们在图7 中画出黑色实线的箭头环线, 以此反映中国东部逐日降水演变的最大概率, 也反映天气系统的转换规律。

3.2 降水异常型对应的高度异常

预报员通常关注在怎样的气压场分布下会出现图 7 中各降水型。我们很容易根据出现 1 型的那些天数的 850 hpa位势高度扰动, 合成出对应的高度异常场(图8 中 1 型)。显而易见, 图7 和 8 中 1 型的江淮流域降水偏多带是850 hpa 的一个负扰动中心在黄海向西南延伸的扰动槽的产物。这个高度扰

动中心和扰动槽在图8中的指向(灰色虚线箭头)与多雨带的指向一致。图8反映降水偏多区域对应的低层大气高度扰动负异常与2012 年 7月 21日北京特大暴雨时出现的瞬时高度扰动垂直剖面低层负值 分布[42]是一致的。3.3 高度异常型对应的降水异常从图 7 到图 8, 反映从降水异常到环流异常的正向对应关系, 即降水异常可以用气压(或环流)异

常来解释。我们再考察一个反向问题, 即由高度扰动判断对应的降水异常。这相当于, 如果数值模式已预报未来时刻的高度异常, 是否有对应的降水异常。

图 9给出中国东部地区 850 hpa 高度扰动场 SOM的 16个(44)分型结果。在左边一列从上向下1→2→3→4 型的顺序中, 中国大陆受到的正高度扰动影响在减弱, 而每一行从左到右的负高度扰动在增强, 右边一列的负高度扰动最强。对各型高度场

进行降水距平的合成分析(图 10), 发现 13 型和 16型的合成降水偏多, 这两个降水型均位于右边一列,有较强的负高度扰动。然而, 其他负高度扰动型并不对应强降水, 再次说明发生强降水时总是可以找 到负高度扰动的信号(图 7 和 8)。但是, 用单一变量的高度扰动不能解释所有降水事件(图 9 和 10),原因是高度扰动只是产生强降水的条件之一, 其他条件包括水汽扰动、温度扰动和风扰动等[43–44]。

结果和讨论

本文基于 1961—2010年中国东部夏季降水格点资料, 分析SOM方法在我国降水分型中的应用,并与传统的EOF分解得到的结果做比较。在此基础上, 探讨SOM分型中降水型的年际变化及转换关系。

对中国东部夏季季节平均降水距平, 使用 SOM分型可得到3组两两相反的降水型, 它们的排序在空间分布上有较好的规律和连续性, 易于理解。与SOM不同, EOF 分解得到的前3个模态只能解释约42%的方差, 不能全面地反映中国东部地区降水的真实情况。通过比较SOM分型与 EOF分解结果,发现SOM分型能体现EOF的前两个模态, 其时间

序列也验证了前人的结论。EOF第三模态在降水分型中是不存在的。SOM方法得到的夏季降水分型能够很好地用850 hpa的风扰动和700 hpa 的水汽比湿扰动解释。作为倒向问题, 强降水的发生可以从850 hpa的高度扰动负值中心得到解释。但是,作为正向问题, 数值模式预报得到的高度负值扰动不一定对应强降水。这是由于强降水的影响因素除高度扰动外, 还有水汽和温度扰动等。

本文初步揭示了 SOM 方法在中国东部夏季降水分型中的优势。在对降水型的类型转换分析方面, 将来可以深入研究转换的时间尺度, 并解释转换的物理意义。在本文探讨的高度场扰动指示降水异常之外, 还可以寻找与预报相关的其他先期信号,如风场异常、温度扰动和水汽扰动等。 [1] 黄荣辉, 徐予红, 周连童. 我国夏季降水的年代际465‒476变化及华北干旱化趋势. 高原气象, 1999, 18(4): [2] 蔡榕硕, 谭红建, 黄荣辉. 中国东部夏季降水年际

35‒46变化与东中国海及邻近海域海温异常的关系. 大气科学, 2012, 36(1):

1‒9 [3] 廖荃荪, 赵振国. 我国东部夏季降水分布的季度预报方法. 应用气象学报, 1992, 3(增刊 1):

15‒19 [4] 魏凤英, 张先恭. 我国东部夏季雨带类型的划分及预报. 气象, 1988, 14(8):

329‒341 [5] 王绍武, 赵宗慈. 近五百年我国旱涝史料的分析.地理学报, 1979, 34(4): [6] 孙林海, 赵振国, 许力, 等. 中国东部季风区夏季

56‒62雨型的划分及其环流成因分析. 应用气象学报, 2005, 16(增刊 1): [7] 黄荣辉, 陈际龙, 黄刚, 等. 中国东部夏季降水的545‒560准两年周期振荡及其成因. 大气科学, 2006, 30(4):

46‒50 [8] 龚振淞, 杨义文. 中国夏季旱涝气候预测相似模型.气象, 2010, 36(5):

331‒340 [9] 李栋梁, 王文. 中国西北夏季降水特征及其异常研究. 大气科学, 1997, 21(3): [10] 宋正山, 杨辉. 夏季东亚季风区 500 hpa 月环流异

401‒404常及与我国降水关系的向量 EOF 分析. 大气科学, 2001, 25(3): [11] 黄山江, 王谦谦, 刘星燕. 西北地区春季和夏季降336‒346水异常特征分析. 大气科学学报, 2004, 27(3): [12] Zhou T J, Yu R C. Atmospheric water vapor transport

associated with typical anomalous summer rainfall

211‒211 patterns in China. Journal of Geophysical Research Atmospheres, 2005, 110(8): [13] 庞轶舒, 祝从文, 刘凯. 中国夏季降水异常 EOF 1137‒1146模态的时间稳定性分析. 大气科学, 2014, 38(6): [14] Dommenget D, Latif M. A cautionary note on the interpretation 216‒225 of EOFS. Journal of Climate, 2002, 15(2): [15] Hewitson B C. Neural nets: applications in geography (Vol. 29). New York: Springer Science & Business Media, 1994 [16] Hewitson B C, Crane R G. Self-organizing maps: applications

13‒26 to synoptic climatology. Climate Research, 2002, 22(1): [17] Cavazos T. Using self-organizing maps to investigate extreme climate events: an application to wintertime

1718‒1732 precipitation in the Balkans. Journal of Climate, 2000, 13(10): [18] Malmgren B A, Winter A. Climate zonation in Puerto Rico based on principal components analysis and an

977‒985 artificial neural network. Journal of Climate, 1999, 12(4): [19] Ambroise C, Sèze G, Badran F, et al. Hierarchical

47‒52 clustering of self-organizing maps for cloud classification. Neurocomputing, 2000, 30(1): [20] Liu Y, Weisberg R H. A review of self‐ organizing map applications in meteorology and oceanography // Mwasiagi J I. Self-organizing maps: applications and novel algorithm design. Rijeka: Intech, 2011: 253– 272 [21] Tadross M A, Hewitson B C, Usman M T. The interannual variability of the onset of the maize

3356‒3372 growing season over South Africa and Zimbabwe. Journal of Climate, 2005, 18(16): [22] Nishiyama K, Endo S, Jinno K, et al. Identification of typical synoptic patterns causing heavy rainfall in the

185‒200 rainy season in Japan by a self-organizing map. Atmospheric Research, 2007, 83(2): [23] Lin Gwo-fong, Chen Lu-hsien. Identification of homogeneous regions for regional frequency analysis

1‒9 using the self-organizing map. Journal of Hydrology, 2006, 324(1): [24] Hsu Kuo-chin, Li Sheng-tun. Clustering spatialtemporal precipitation data using wavelet transform

190‒200 and self-organizing map neural network. Advances in Water Resources, 2010, 33(2):

[25] Liu Y G, Weisberg R H. Patterns of ocean current variability on the West Florida Shelf using the selforganizing map. Journal of Geophysical Research: Oceans, 2005, 110: C06003 [26] Liu Y G, Weisberg R H. Ocean currents and sea surface heights estimated across the West Florida 1697‒1713 Shelf. Journal of Physical Oceanography, 2007, 37(6): [27] Liu Y G, Weisberg R H, Mooers C N. Performance evaluation of the self-organizing map for feature extraction. Journal of Geophysical Research: Oceans, 2006, 111: C05018 [28] Reusch D B, Alley R B, Hewitson B S. Relative performance of self-organizing maps and principal component analysis in pattern extraction from synthetic 188‒212 climatological data. Polar Geography, 2005, 29(3):

19‒36 [29] Lobo V J. Information fusion and geographic information systems. Berlin: Springer, 2009: [30] Solidoro C, Bandelj V, Barbieri P, et al. Understanding dynamic of biogeochemical properties in the northern Adriatic Sea by using self-organizing maps and k-means clustering. Journal of Geophysical Research, 2007, 112: C07S90.1–C07S90.13 [31] Qian Weihong, Shan Xiaolong, Liang Haoyuan, et al. A generalized beta-advection model to improve unusual typhoon track prediction by decomposing total flow into climatic and anomalous flows. Journal of Geophysical Research: Atmospheres, 2014, 199(3): 1097–1117 [32] Pearson K. On lines and planes of closest fit to systems of points in space. The London, Edinburgh

559‒572 and Dublin Philosophical Magazine and Journal of Science, 1901, 2(11): [33] Lorenz E N. Empirical orthogonal functions and statistical weather prediction [R]. Science Report 1. Cambridge, Massachusetts: Statistical Forecasting Department of Meteorology, MIT, 1956: NTIS AD 110268 [34] Kohonen T. Self-organized formation of topologically

59‒69 correct feature maps. Biological Cybernetics, 1982, 43(1):

501‒505 [35] Kohonen T. Self-organizing maps. 3rd ed. New York: Springer, 2001: [36] Vesanto J, Himberg J, Alhoniemi E, et al. SOM toolbox for Matlab 5 [R]. Finland: Helsinki Univ of Technol, 2000 [37] Christiansen B. Atmospheric circulation regimes: can

2229‒2250 cluster analysis provide the number?. Journal of Climate, 2007, 20(10): [38] Michelangeli P A, Vautard R, Legras B. Weather

1237‒1256 regimes: recurrence and quasi stationarity. Journal of the atmospheric sciences, 1995, 52(8): [39] Riddle E E, Stoner M B, Johnson N C, et al. The impact of the MJO on clusters of wintertime circulation anomalies over the North American region. Climate Dynamics, 2013, 40: 1749–1766 [40] Chattopadhyay R, Vintzileos A, Zhang C. A description of the Madden-julian Oscillation based on a 1716‒1732 self-organizing map. Journal of Climate, 2013, 26(5):

289‒295 [41] 邓爱军, 陶诗言, 陈烈庭. 我国汛期降水的 EOF 分析. 大气科学, 1989, 13(3): [42] Jiang Ning, Qian Weihong, Du Jun, et al. A comprehensive approach from the raw and normalized anomalies to the analysis and prediction of the Beijing

1551‒1567 extreme rainfall on July 21, 2012. Natural Hazards, 2016, 84(3): [43] Qian Weihong. Temporal Climatology and anomalous weather analysis. Singapore: Springer, 2017 [44] Qian Weihong, Jiang Ning, Du Jun. Anomaly based weather analysis vs. traditional total-field based weather analysis for depicting regional heavy rain events. Weather and Forecasting, 2016, 31(1): 71–93

图 1深度图计算流程Fig. 1 Depth estimation workflow

虚线部分表示不考虑微透镜和成像平面时光线的汇聚图 2全光相机结构Fig. 2 Structure of the plenoptic camera

(a) 绿色点代表相邻影像的中心点, 蓝色十字代表同名点的位置, (b) 同名点一般沿着核线方向进行搜索图 4多基线同名点匹配Fig. 4 Multi-baseline corresponding points matching

图 3微透镜–探测器系统Fig. 3 Microlens-detector system

图 5 深度图代价计算Fig. 5 Cost calculation for depth map

(a) 彩色全聚焦影像; (b) 投影的 WTA 深度图; (c) DL=50优化的深度图; (d) MRFSP 优化的深度图; (e) 文献[24]得到的深度图; (f) 图(b)~(e)中黑色矩形内放大的深度图 图 9钢结构件的深度图Fig. 9 Depth map for a workpiece

(a) 彩色全聚焦影像; (b) 投影的 WTA 深度图; (c) DL=50优化的深度图; (d) MRFSP 优化的深度图; (e) 文献[24]得到的深度图; (f) 图(b)~(e)中黑色矩形内放大的深度图 图 10人脸的深度图Fig. 10 Depth map for a woman face

(a) 彩色全聚焦影像; (b) 投影的 WTA 深度图; (c) DL=50优化的深度图; (d) MRFSP 优化的深度图; (e) 文献[24]得到的深度图; (f) 图(b)~(e)中黑色矩形内放大的深度图 图 11水果的深度图Fig. 11 Depth map for a fruit pile

图 1 2016 年 3 月 20日激光雷达归一化后向散射信号Fig. 1 Lidar normalized relative backscatter on March 20th, 2016

图 2 2016 年 3 月 20 日云检测识别结果(蓝色为云) Fig. 2 Cloud identification results on March 20th, 2016 (cloud in blue)

图 4 典型冰云的退偏比变化特征(2016 年 3 月 20 日 04:00) Fig. 4 Typical depolarization raito of ice cloud (20160320 04:00)

图 3典型水云的退偏比变化特征(2016 年 5 月 14 日 05:00) Fig. 3 Typical depolarization raito of water cloud (20160514 05:00)

图 5 水平导向的片状冰晶退偏比的变化特征(2016 年 3 月 22 日 16:00) Fig. 5 Depolarization raito of horizontal oriented ice flakes (20160322 16:00)

图 9各种类型云层月平均高度变化趋势Fig. 9 Monthly-averaged cloud layer height of various types of clouds

图 10各种类型的云月均出现时间占比Fig. 10 Monthly-averaged percentages of occurrence of various types of clouds

图 11 各月份过冷水云在 0~−40℃云中的月均占比Fig. 11 Monthly-averaged percentages of occurrence of supercooled water clouds among clouds with temperature from 0℃ to −40℃

右下角数字为该模态的解释方差; 实线表示降水偏多, 虚线表示降水偏少; 等值线范围为−3~3 mm, 间隔 0.5 mm图 1 1961—2010 年中国东部夏季(JJA)降水距平场前 3 个 EOF 模态空间分布(单位: mm) Fig. 1 Spatial distributions of the first three EOF modes of summer (JJA) precipitation anomaly over Eastern China from 1961 to 2010 (unit: mm)

右下角数字为该降水型发生的频率; 实线表示降水偏多, 虚线表示降水偏少; 等值线范围为−3~3 mm, 间隔 0.5 mm图 2 1961—2010 年中国东部夏季降水距平场 6 个(3×2) SOM 分型Fig. 2 Six (3×2) SOM modes of summer precipitation anomaly over Eastern China from 1961 to 2010

图 3中国东部夏季降水前 3 个 EOF模态相对应的主分量时间序列Fig. 3 Time series of the first three EOF modes for the summer precipitation over Eastern China

右上角数字为该降水型发生的频率; 实线代表正数, 虚线代表负数, 等值线间隔 0.16 m/s图 5 SOM各个降水型对应的 850 hpa 扰动风(m/s)和 700 hpa 水汽比湿扰动(g/kg)的合成Fig. 5 Composites of 850 hpa wind vectors anomaly (m/s) and 700 hpa specific humidity anomaly (g/kg) of the six SOM modes

图 4 6 个(3×2) SOM 分型对应的 BMU时间变化序列Fig. 4 Temporal evolution of the BMU for the six (3×2) SOM modes

右下角数字为该模态的解释方差; 实线代表正数,虚线代表负数, 等值线间隔 0.01 mm/d图 6 1981—2010 年中国东部夏季日平均降水距平场前6 个 EOF 模态空间分布Spatial distributions of the six leading EOF modes of daily summer precipitation anomaly over Eastern China from 1981 to 2010

右上角数字为该型发生的频率; 实线表示降水偏多, 虚线表示降水偏少; 等值线范围为−30~30 mm/d,头指示降水偏多的地理位置, 黑色实线箭头指示后一天降水异常型的最大概率转移路径图 7 1981—2010 年中国东部夏季日平均降水距平场 16 个(4×4) SOM 分型结果Sixteen (4×4) SOM modes of daily summer precipitation anomaly over Eastern China from 1981 to 2010 Fig. 7 间隔 5 mm/d; 灰色虚线箭

实线表示正高度扰动, 虚线表示负高度扰动; 等值线范围为−30~30 gpm, 间隔 5 gpm图 8 降水距平 SOM 分型的 850 hpa 位势高度扰动(gpm)合成分析Composites of 850 hpa geopotential height anomaly (gpm) of the SOM patterns Fig. 8

右上角数字为该型发生的频率; 实线表示正高度扰动, 虚线表示负高度扰动; 等值线范围为−30~30 gpm, 间隔 5 gpm图 9 中国东部地区 850 hpa 日平均高度场扰动(gpm)场 SOM 的 16 个(44)分型结果Sixteen (4×4) SOM modes of daily 850 hpa geopotential height anomaly (gpm) over Eastern China Fig. 9

实线表示降水偏多, 虚线表示降水偏少; 等值线范围为−30~30 mm/d, 间隔 5 mm/d图 10 中国东部地区 850 hpa 高度扰动场 SOM 分型(图 9)对应的日降水距平(mm/d)合成Fig. 10 Composites of daily precipitation anomaly (mm/d) of each daily 850 hpa geopotential height anomaly SOM modes over Eastern China in Fig. 9

Newspapers in Chinese (Simplified)

Newspapers from China

© PressReader. All rights reserved.