2022, Vol. 38
2. 成都信息工程大学，四川 成都 610225;
3. 中国气象局广州热带海洋气象研究所，广东 广州 510641;
4. 中国铁道科学研究院集团有限公司通信信号研究所，北京 100081
2. Chengdu University of Information Technology, Chengdu 610225, China;
3. Guangzhou Institute of Tropical and Marine Meteorology, CMA, Guangzhou 510641, China;
4. Institute of Communication and Signaling, China Academy of Railway Sciences Group Co., Ltd., Beijing 100081, China
闪电是一种产生于大气中的高能量、长距离放电现象，其出现往往伴有强烈的闪光[1]。闪电的光学观测可以帮助人们直观地了解云外闪电放电通道的发展情况，一直是研究闪电发展过程的重要手段。
随着科学技术的发展，各式各样的光学观测设备不断被应用于闪电观测，如胶片相机、磁带式录像机、数码相机、数字摄像机、条纹相机[2-3]、高速摄像机[4]等，为研究人员分析闪电特征提供了丰富的观测资料。为保证闪电光学观测有足够大的视野，提升捕捉闪电通道的概率，经常会采用短焦距的广角镜头，甚至是鱼眼镜头。尽管广角镜头有成像视角范围大的优点，但它获取的图像通常存在着畸变，目前闪电光学观测资料的分析大多未考虑镜头成像畸变的校正，镜头成像畸变对闪电特征分析结果的影响如何迄今仍缺乏准确的认识，在一定程度上会影响闪电通道发展特征定量分析的准确性，迫切需要对光学观测设备获取的图像进行校正以消除畸变的影响。
摄像机标定方法[5]目前主要分为传统标定法[6-8] 和自标定法[9] 两类。传统标定方法需要预先准备高精度的标定物作为空间参照物，建立三维空间中标定物上的特征点与获取到二维图像中的对应点之间的映射关系，以此计算摄像机的各项参数。这种方法精度高，但标定过程复杂；自标定方法不必利用标定物，只利用二维图像信息来恢复摄像机内外参数，虽然不需要特定的标定物，然而随之带来的问题则是标定稳定性不足精准度较差。张正友等提出的平面标定法[10-13]，这种方法介于传统标定法和自标定法之间，既克服了传统标定法需要的高精度三维标定物的缺点，又解决了自标定法鲁棒性差的难题，并且只需拍摄若干幅不同角度或位置下的平面标定板图像，因此操作起来相对较为简便，在实际中得到了很广泛的应用。
本研究基于张正友摄像机平面标定方法制作了便于现场标定实验的移动棋盘格平面靶标，并针对TOLOG的多套光学观测设备拍摄的靶标图像，使用特征点提取算法获取特征点的坐标信息，并通过坐标之间的位置信息计算出多种参数，最终完成设备的标定。利用设备标定得到的参数对图像进行校正，重点研究不同设备获取的图像受畸变影响的程度，从而实际应用到闪电光学观测资料当中，消除设备对闪电通道特征分析的影响。
2 标定原理 2.1 镜头畸变由于镜头成像过程中产生的光学畸变，导致获取图像的非线性失真。畸变的产生使得原本沿直线投影的光路发生偏移，物点、光心以及像点三点不共线。通常摄像机镜头的畸变分为径向畸变、切向畸变[14-16] 等。
图像中常见的失真主要是桶形畸变与枕形畸变造成的，他们都属于径向畸变。径向畸变的成因是镜头透镜近光轴区域和远光轴区域放大倍率不同所导致，可以用下式表达：
| $ x_{c}=x_{u}+x_{u} \cdot\left(k_{1} r^{2}+k_{2} r^{4}+k_{3} r^{6}+\cdots \cdots\right) $ | (1) |
| $ y_{c}=y_{u}+y_{u} \cdot\left(k_{1} r^{2}+k_{2} r^{4}+k_{3} r^{6}+\cdots \cdots\right) $ | (2) |
式中，假设图像中一像点坐标为xu, yu，其在理想条件下无畸变的坐标为xc, yc，由此可以用k1，k2，k3，……表示径向畸变的系数，而像点与光心的距离用
切向畸变的出现主要是由于镜头模组中各镜片与像平面不严格平行造成的，这类缺陷可以表示为：
| $ \begin{aligned} x_{c}=x_{u}+\left[p_{1}\left(r^{2}+2 x_{u}{ }^{2}\right)+2 p_{2} x_{u} y_{u}\right] \end{aligned} $ | (3) |
| $ y_{c}=y_{u}+\left[2 p_{1} x_{u} y_{u}+p_{2}\left(r^{2}+2 x_{u}{ }^{2}\right)\right] $ | (4) |
切向畸变的系数用p1和p2来表示。
2.2 靶标的设计与特征点的提取为了实现摄像机的标定，寻求世界坐标系中特征点与图像中对应点的映射关系，就需要制作供摄像机拍摄提取特征点的靶标，用于标定摄像机的各项参数。靶标在保证制作精度的同时还须满足以下两个基本条件：首先靶标特征点之间相对位置关系已知，其次图像特征点的坐标容易求取。实验中选取了制作简便且精度较高的平面靶标[17]，并将其设计成棋盘格样式，靶标采用高精度打印技术并配合高清照相纸，如图 1所示，表面上均匀分布着黑白相间的正方形棋盘网格，其中每块小方格尺寸为50 mm×50 mm。靶标上每行每列均为12个小方格，总尺寸为600 mm×600 mm。
|
图 1 平面棋盘格靶标及特征点(角点)选取 |
选取特征点时应遵循坐标点稳定且容易提取的原则。在研究中选用了相邻四个棋盘格的十字中心作为特征点，称之为角点。这样12×12的小方格能组成11×11共121个角点(最外围的一圈小方格只能与其内侧组成角点)。靶标每个角点之间的设计尺寸也是已知的，为了便于移动携带试验，进而将一块65 cm×65 cm的钢化玻璃作为基底并把靶标张贴其上，通过空间中靶标与获取图像中角点之间的映射可以推算出特征点的相对关系。将拍摄到的图进行片二值化处理，采用Harris角点检测算子[18] 得到角点的初始坐标值[19]，再将这些坐标值进行亚像素级定位进一步减小噪声对标定结果的影响，这样就得到了这些特征点在图像中的实际坐标。根据坐标即可求出相机的内部参数，利用最小二乘法便可求解出畸变参数k1，k2。
3 闪电光学观测系统标定本文用以比对分析的闪电通道光学观测资料是由位于广州的高建筑物雷电观测站[20-26] (TallObject Lightning Observatory in Guangzhou, 简称TOLOG)的数字化高速摄像机以及单反相机所获取。广州高建筑物雷电观测站目前由1个主站和多个光学观测子站组成。本研究主要对主站的6套光学设备逐一进行标定。
(1) 3套高速摄像系统：摄像机均为Photron公司的FASTCAM系列高速摄像机：1台型号为SAZ，搭载焦距14 mm的Nikon f/2.8 D ED AF镜头，为叙述方便，本文采用编号SAZ代表此套高速摄像系统；1台SA3，搭载Kowa 8 mm f/1.4镜头，编号SA3；1台SA5，搭载Nikon20 mm f/1.6镜头，编号SA5。
(2) 3套单反相机：两台Nikon D7100机身，分别搭配Sigma 8 mm f/3.5 EX DG Circular Fisheye镜头(编号：D7100 - 8)和Nikon AF Fisheye - Nikkor 16 mm f/2.8D镜头(编号：D7100-16)；一台机身为D7000，搭配产自Nikon AF Nikkor 20 mm f/2.8 D镜头(编号：D7000-20)。
如图 2所示，用于TOLOG主站光学设备标定的平面棋盘格标定靶标(图 2a)以及设备的拍摄视野(图 2b)，利用张正友校正方法移动拍摄不同角度的多张(不少于20张)靶标照片，使得校正板充分分布在图片画面视野的各个位置使靶标在图片中占据视野的不同位置这样能更好地提高标定精度，鲁棒性更强。标定后得到各类参数见表 1。
|
图 2 用于相机标定的靶标照片(a)及广州高建筑物雷电观测站光学观测设备的视野(b) |
|
表 1 TOLOG闪电光学观测设备的标定参数 |
表 1中fx和fy为图像像素分布的比例因子，并且畸变系数只有k1和k2两项，由此可见实验设备获取的图像只受到了径向畸变的影响。使用这些参数对靶标图像进行校正，进而提取畸变校正后其横、纵各行列的角点坐标，并用坐标信息拟合出二阶多项式，多项式高次项系数可以达到10-6数量级，由此可见算法的标定误差很低。同时，使用表 1中的标定数据分别对6套闪电光学观测设备获取的图像进行畸变校正，计算出畸变校正前后图像中像点与光心在不同视角角度处的变化差异，并绘制出散点图进而分别拟合成曲线。如图 3所示，图像在光轴与视线夹角小于5 °时基本没有受到畸变的影响，但随着夹角的增大，畸变校正前后图像像点与光心的角度也在不断增大(畸变影响变大)。在这6套设备中两套搭配鱼眼镜头的设备产生的畸变导致图像的形变以及失真最大，除此以外其他设备获取的图像受畸变影响的程度随着搭配镜头焦距的增长相对应的变小。
|
图 3 闪电光学观测设备图像校正前后像点与光心的角度受畸变影响的相对大小 |
2018年7月6日16时41分TOLOG观测到一次广州塔上行闪电个例，该个例分别被多台光学仪器所捕获。图 4a为主观测站D7100-8的单反相机拍摄，上行闪电由广州塔顶端始发，向上发展到最高处(在D7100-8获取的图像中约为1 500像素的高度)后向左拐转向视野的左下方继续发展。图中较为完整的记录了此次上行闪电的通道，能够较为全面地分析畸变对不同区域内闪电通道的影响。图 4b是结合表 1中的参数经过校正得到的去除畸变的图像。根据镜头的成像原理图像光心位置畸变为零[20]，为了精确对比闪电通道图像畸变校正前后的差异，在原始图像的正中心位置添加了红色十字标识(图 4a)，在此基础上对图像进行畸变校正得到校正后图像(图 4b)，将校正前后的两幅图像以光心为基准(十字像素点重合)叠加，由此能够得到闪电通道校正前后的叠加对比图。由图 4c可见，图像在叠加后进行对比度增强及灰度化处理以便于通道的辨认，之后又分别把图像校正前后的闪电通道标记成蓝色和红色能够清晰地看出校正前后闪电通道差异。方框选中的区域是广州塔的塔尖，从中可以看到广州塔塔尖在图像校正前后的偏移大概在11个像素，经过换算在空间中偏移量能达到17.5 m。
|
图 4 TOLOG主站单反相机(D7100-8)于2018年7月6日拍摄的一次广州塔上行闪电个例 (a)为资料原图，(b)是将(a)经过摄像机标定后去除畸变的照片。(c)为(a)(b)以光心为基点的叠加图，并分别将校正前后闪电通道标记成蓝色和红色(图c进行了对比度增强)。 |
畸变校正前后闪电通道长度也有一定的差异，由图 4c可见，校正后的闪电通道(红色)比校正前(蓝色)的要长。通过计算，校正前闪电通道长度约为3 654个像素，校正后通道长度约为4 125个像素，两者像素个数相差471，校正后闪电通道长度较原图增加了12.9%。
为了探求闪电通道位于图像中不同位置受到畸变影响的强弱，如图 5所示，以图像光心为圆心，均匀绘制20个同心圆覆盖全部闪电通道，将其分为37段并保证每段闪电通道的位置关系已知。计算每段闪电通道校正前后的长度差值，除以校正后的长度，得到相对长度差，结合各同心圆区域到图像光心的夹角，绘制出散点图 6。根据散点拟合出一条曲线，随着闪电通道与图像光心的距离越来越远，径向畸变对校正前后闪电通道之间的长度影响也越来越大，通过拟合曲线能够看出畸变呈现出非线性的变化，并且畸变对图像的影响是随着图像从光心往边缘区域处的夹角增大而增大的，并在视野边缘处达到最大。
|
图 5 以图像光心为圆心的闪电通道同心圆分布图 |
|
图 6 校正前后区域内闪电通道相对长度随区域与图像光心的夹角的关系 |
图 7为观测站高速摄像机SA3拍摄的此次闪电事件的同步资料闪电通道图像，其中图 7a为高速摄像机拍摄的原始图像，经过畸变校正后得到图 7b，再将校正前后图像以光心为基准叠加，得到校正前后闪电通道叠加对比图(图 7c)。图 7c中以红色和蓝色分别填充了校正后的闪电通道和原始图像中的闪电通道。经过计算，SA3获取图像畸变校正前后闪电通道长度相差4.5%。
|
图 7 TOLOG主站高速摄像机(SA3)捕捉到的广州塔上行闪电光学同步观测资料 (a)、(b)分别为资料原图和畸变校正后的图像。(c)为(a)(b)以光心为基点的叠加图，并分别将校正前后闪电通道标记成蓝色和红色(图c进行了对比度增强、灰度化处理)。 |
为了定量分析这两台搭配不同种类镜头的光学设备获取图像中闪电通道长度的差异，将D7100-8图像通道末端清晰处的拐点作为特征点(如图 8左图红圈所示)，同时高速摄像机SA3图像中选取闪电通道临近末尾处相同位置的拐点作为特征点(见图 8右图红圈所示)，通过像素点间距离[22]换算分别计算两套设备从广州塔塔尖始发到特征点处闪电通道长度，可以得出：SA3原始图像中这段通道长度为5 971 m与D7100-8原始图像的5 425 m相差10.1%，而SA3图像校正后的通道长度为6 236 m与D7100-8所获取原始图像的6 124 m相差1.8%。由此可见不同光学设备畸变校正后的图像一致性更好，所以对光学设备的标定校正是很有必要的。
|
图 8 在不同设备获取的同步观测资料图像中选取闪电通道的相同拐点作为特征点 |
本文基于张正友平面标定法运用便携式棋盘格平面靶标，实现了一种既适用于实验室研究，又可以应用在雷电观测站外场摄像机标定的灵活方法[10-11]。在实际运用中通过移动平面靶标校正了高建筑物雷电观测站的各个光学观测仪器逐一标定，结合闪电光学资料数据，分析了一次广州塔上行闪电畸变校正前后图像闪电通道长度差异。具体结论如下。
(1) 通过算法标定得到观测站各光学系统的参数，分别对6套闪电光学观测设备获取的图像进行畸变校正，通过对校正前后图像比对后发现，图像光心处可看做0畸变的基准点。除此之外图像受畸变影响的程度跟拍摄设备搭载镜头的种类和焦距有着密切联系，对于同一种类的镜头，焦距越短，图像中像点与光心角度的差异也就越大，畸变也随着图像视场角的增大而愈加强烈。
(2) 对一次广州塔上行闪电个例图像进行畸变校正，D7100-8获取的原始图像闪电通道长度为：3 654个像素，校正后通道长度为：4 125个像素。校正前后闪电通道长度相差12.9%；通过对高速摄像机SA3获取的同一闪电事件的同步光学观测资料进行畸变校正，校正前后闪电通道的程度相差4.5%；并对同一个个例不同光学系统获取的闪电通道长度进行定量对比后发现，经过畸变校正后不同设备获取的闪电通道图像比原图拥有更好的一致性，由此可见校正前后闪电通道长度的差异可能会对闪电二维速度的计算带来一定影响。
(3) TOLOG六个观测站点对发生在广州珠江新城地区的闪电事件进行同步观测，利用其中两个站点拍摄的闪电同步光学观测资料，可以对闪电通道进行三维重建。目前，在闪电通道三维重建之前并没有考虑镜头畸变对图像的影响[27]。广角镜头以及鱼眼镜头获取的闪电通道二维图像畸变校正前后差异明显，可以预计，利用未经畸变校正的图像直接进行三维重建得到的闪电通道三维坐标会存在更大的偏差，这还需要进一步深入的研究。
| [1] |
郄秀书, 张其林, 袁铁, 等. 雷电物理学[M]. 北京: 科学出版社, 2013.
|
| [2] |
SCHONLAND, B F J., MALAN, D J J, COLLENS, H. 1935 Progressive lightning Ⅱ[M], Proc. R. Soc. London, Ser. A, 152, 595 625, doi: 10.1098/rspa.1935.0210.
|
| [3] |
BERGER, K, VOGELSANGER E. Photographische Blitzuntersuchungen der Jahre 1955-1965 a uf dem Monte San Salvatore[J]. Bull Schweiz Elektrotech Ver, 57, 599 - 620.
|
| [4] |
LU W, CHEN L, MA Y, et al. Lightning attachment process involving connection of the downward negative leader to the lateral surface of the upward connecting leader[J]. Geophys Res Lett, 2013, 40(20): 5 531-5 535. DOI:10.1002/2013GL058060 |
| [5] |
马颂德, 张正友. 计算机视觉——理论与算法基础[M]. 北京: 科学出版社, 1998.
|
| [6] |
ABDEL-AZIZ Y I, KARARA H M. Direct Linear Transformation from Comparator Coordinates into Object Space Coordinates in CloseRange Photogrammetry[J]. Photogrammetric Engineering & Remote Sensing, 81(2): 103-107.
|
| [7] |
STURM P F, MAYBANK S J. On Plane-Based Camera Calibration: A General Algorithm, Singularities, Applications[C]. IEEE Computer Society Conference on Computer Vision & Pattern Recognition, 1999.
|
| [8] |
TSAI R. A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses [J]. 1987, RA3(4): 323-344.
|
| [9] |
孟晓桥, 胡占义. 摄像机自标定方法的研究与进展[J]. 自动化学报, 2003, 29(1): 110-124. |
| [10] |
ZHANG Z Y. A flexible new technique for camera calibration[R]. Microsoft Corporation, NSR-TR-98-71, 1998.
|
| [11] |
ZHANG Z Y. Flexible Camera Calibration by Viewing a Plane from Unknown Orientations[C]//Proceedings of IEEE International Conference on Computer Vision, Corfu, Greece, Sep 20-27, 1999, 1: 1-8.
|
| [12] |
吴福朝. 计算机视觉中的数学方法[M]. 北京: 科学出版社, 2008: 46-52.
|
| [13] |
杨必武, 郭晓松. 摄像机镜头非线性畸变校正方法综述[J]. 中国图象图形学报, 2005, 10(3): 269-274. DOI:10.3969/j.issn.1006-8961.2005.03.001 |
| [14] |
BROWN D C. Close-range camera calibration[J]. Photogrammetric Engineering, 1971, 37(8): 855-866. |
| [15] |
张广军. 视觉测量[M]. 北京: 科学出版社, 2008: 34-42.
|
| [16] |
WEI G, MA S. Implicit and explicit camera calibration: Theory and experiments[J]. Transactions on Pattern Analysis and Machine Intelligence, 1994, 16(5): 469-480. DOI:10.1109/34.291450 |
| [17] |
余明浪, 魏振忠, 孙军华. 基于柔性平面靶标的摄像机现场标定方法[J]. 北京航空航天大学学报, 2009(3): 84-87. |
| [18] |
HARRIS C, STEPHENS M. A Combined Corner and Edge Detection[C]//InProc. 4th Alvey Vision Conel988: 189-192.
|
| [19] |
WEI G Q, MA S D. Implicit and Explicit Camera Calibration: Theory and Experiments[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1994, 16(5): 469-480. DOI:10.1109/34.291450 |
| [20] |
吕伟涛, 陈绿文, 马颖, 等. 广州高建筑物雷电观测与研究10年进展[J]. 应用气象学报, 2020, 31(2): 129-145. DOI:10.11898/1001-7313.20200201 |
| [21] |
武斌, 吕伟涛, 齐奇, 等. 双向先导正端突然延展现象的高速摄像观测[J]. 应用气象学报, 2020, 31(2): 146-155. DOI:10.11898/1001-7313.20200202 |
| [22] |
齐奇, 吕伟涛, 武斌, 等. 广州两座高建筑物上闪击距离的二维光学观测[J]. 应用气象学报, 2020, 31(2): 156-164. DOI:10.11898/1001-7313.20200203 |
| [23] |
QI Q, LU W, MA Y, et al. High-speed video observations of the fine structure of a natural negative stepped leader at close distance[J]. Atmospheric Research, 2016, 178-179: 260-267. DOI:10.1016/j.atmosres.2016.03.027 |
| [24] |
WU B, LYU W, QI Q, et al. Synchronized Two-Station Optical and Electric Field Observations of Multiple Upward Lightning Flashes Triggered by a 310-kA +CG Flash[J]. Journal of Geophysical Research: Atmospheres, 2019, 124(2): 1 050-1 063. |
| [25] |
QI Q, LYU W, MA Y, et al. High-speed video observations of natural lightning attachment process with framing rates up to half a million frames per second[J]. Geophysical Research Letters, 2019, 46(12): 12 580-12 587. |
| [26] |
迟泽英, 陈文建. 应用光学与光学设计基础[M]. 北京: 高等教育出版社, 2013.
|
| [27] |
WILLETT J, LE VINE D, IDONE V. Lightning return stroke current waveforms aloft from measured field change, current, and channel geometry[J]. Journal of Geophysical Research: Atmospheres, 2008, 113(D7): D07305. |