陈青青1,2，汤志亚1，杨玲1,2

　　（1.成都信息工程大学，四川 成都 610225；2.中国气象局 大气探测重点开放实验室，四川 成都 610225）

       摘要：采用地基镜头旋转扫描拍摄不同俯仰角的天空图像，通过全景图拼接获得全天空图像，避免了鱼眼透视畸变，保留了低仰角信息。目前全景图拼接技术主要局限于单层柱面，通过空间坐标转换及以渐入渐出加权平均融合算法和插值算法，借鉴地图学投影原理，将28张图像在垂直投影面和等积柱面投影面进行无缝拼接得到360°×180°的全天空全景图，为识别云状、等积计算云量提供了连续完整的全天空图像。

　　关键词：地图投影；全天空；云图；全景图拼接

0引言

　　图1原始图像随着大气探测技术的发展，云的观测逐渐自动化，国内外研发了大量不同波段不同视角的测云仪器［1］，全景拍摄方式简单快捷但低仰角畸变大；另外还有全景图拼接方式［2］。本文采用监控摄像机，无盲区扫描拍摄可见光波段天空图像，得到36张具有重叠部分的子图，通过还原物体在天球中的相对位置，并借鉴地图投影转换为平面图像，对重叠部分进行融合，可得到高清连续的全天空云图。

1观测平台

　　云台型号为HY900B，相机为索尼FCBEX1010P，采集器为HB18204T，图像传感器CCD为EXview HAD 1/4英寸，成像面积为3.2 mm×2.4 mm，固定焦距为3.4 mm，视角为57.8°，图像尺寸为704×576，长边视角为44.70°，宽边视角为36.58°，图像上下左右都存在重叠，无盲区。实验用机配置Intel Core i3 CPU，2 GB内存，Windows XP系统。原始图像如图1所示，拍摄俯仰角顺序如表1所示。

　　令拍摄次序t=1~36，其中1~9表示相机光轴仰角为0°的图像，10~18表示仰角30°图像……由此可见第t张图像方位角θt仰角φt与t存在如式(1)所示的关系，其中［·］表示向下取整，%表示求余。

　　相机标定以标准网格角点作为匹配点进行多项式拟合［3］，将拟合关系式应用到整幅图像并对坐标进行双线性插值，截去无效边缘，如图2所示。

　　图2相机标定

　　2图像在球面上的配准

　　图像拼接包括配准和融合两个环节。本文图像有确定的空间坐标，能实现图像在空间位置上的镶嵌，从而找到配准点。

　　以相机位置为原点，拍第一幅图像时相机光轴方向为0°方位角，观测点所在的大地水平面为xoy面，0°方位角为x轴的正向，90°方位角为y轴正向，过原点xoy面指向天空的垂线为z轴。0°方位角图像在空间xoz平面上的剖面图如图3(a)所示。其中E′表示相机光轴投影点E点在xoz面上的投影，在xoz平面上的坐标为(xE,zE),由图3(a)可见。

　　其中，L表示物距，P为图像中任意一点在xoz面上的投影，如图3(a)所示，P在x轴上的坐标为xP=xE+Δx,P在z轴上的坐标为zP=zE－Δz，而

　　由此得出两个关系式：

　　zp=Lsinφt－(j－a)cosφt(5)

　　xp=Lcosφt+(j－a)sinφt(6)

　　联立式（5）、式（6）解得：

　　zp=Lcscφt－xPcotφt

　　且

　　y∈［－b+1,－b+W］

　　x∈［Lcosφt+(1－a)sinφt,Lcosφt+(H－a)sinφt］(7)

　　0°方位角图像在空间中xoy平面的剖面图如图3(b)所示，可见

　　yp=－b+i(8)

　　数字图像坐标是离散值，直接由（j，i，t）计算（x，y，z），会导致有些空间坐标找不到与之对应的图像坐标，投影后图像出现空洞。采用反算插值的方法可避免出现这种情况，即根据投影面上每个点来反查其在原图中的位置和像素值［4］。

　　当相机光轴θt=0°，φt≠0°时，联立式（6）、式（8）反解j、i得：

　　联立式（7）、式（9）即得到θt=0°，φt≠0°时球面各点（x，y，z）所对应的第t张图的图像坐标（j，i，t）。

　　同理，当θt=0°、φt=0°时，联立公式（5）、式（8），反解出

　　联立式（6）、式（10）即得到θt=0°、φt=0°时球面各点对应的第t张图的图像坐标。

　　探测范围方位角0°~360°，仰角0°~180°，半径0~∞,该范围的物体可投影到以相机为球心，给定值为半径的半球面上。假想景物已从空间投影到了天球面上，相机拍摄将天球物体投影到CCD，因此应采用点光源在球心的球面投影方式［5］。

　　设P(xP,yP,zP)是方位角0°，仰角φt的图像上任意一点，利用直角坐标与球坐标的转换公式得到P点的球坐标为(θp,φp,rp)。

　　P点与投影点P′位于同一条直线上，二者的θ和φ相同，而所有投影点的r都等于投影球面半径R，因此P点的投影坐标为P′(θp,φp,R)。

　　云台顺时针旋转，那么球面投影坐标（θp－θt，φp，R）对应第t张图像的图像坐标（j，i）。如式（1）所示，改变θt就能通过方位角0°的图像投影点求得其余8个方位对应的图像投影坐标点，从而得到360°全景图。仰角90°时只取一张图就能覆盖拍摄范围，因此仅需用到前28张图。为避免信息冗余和损失，本文取物距与长边弦长的中值为投影球面半径。

3垂直投影面的拼接

　　垂直投影光源为一束与z轴平行的平行光，投影面为xoy面，投影公式如式（11）所示，（x，y，z）表示空间某点，（x′，y′，z′）为其投影点。

　　图4全天空云图垂直投影实验结果MATLAB仿真垂直投影结果如图4所示，衔接无误位置准确，但边界亮度跃变，这是相机自动曝光引起的。数字图像以矩阵形式存储，反算下标会取整而产生误差，应尽量减少取整次数，在最后一级投影才对下标取整和图像融合。

　　加权平均融合是指在匹配点重合的n个子图各自对应的3个颜色通道的亮度值进行加权平均，新值赋给匹配点，当权值相等时称为均值融合，均值融合结果如图5所示。

　　结果表明均值融合后图像个体完整，但边界明显，原因是没有考虑权值差异，距离匹配点近的子图贡献权值应更大。以距离比为权值，将距离近的图像在该点的值作为主元素进行融合，称为渐入渐出加权平均融合（后面简称加权融合）［6］。

　　子图垂直投影之后的图像是扇形的，由于图像是围绕圆心旋转的，权值便不应用直线距离，而是该点到边界的弧长。对30°层进行加权融合实验结果如图6所示，可见该方法解决了子图边缘突变，效果良好。

　　对60°层进行一次加权融合，仅显示其重叠部分如图7(a)所示，可见重叠部分还有重叠，再次加权融合，结果如图7(b)所示，重叠消失。最后将一级融合后的像元和无重叠像元放入，形成60°层全景如图7(c)所示。

　　先进行同仰角层内融合，再进行层间融合。层间融合用半径的变化量来计算权值，图8（a）是30°、60°层间未融合图像，图8（b）是30°、60°层间加权融合，图8（c）是融合后的全天空图像。

　　4等积柱面投影面的拼接

　　云量是按视面积来计算的，等积投影［5］为云量的计算提供了借鉴。前人证明轴射柱面投影下，任意形状的球面与其在外切圆柱面上投影面积相等［7］。设A是球面上任意一点，A的柱面坐标为(θA,ρA,zA)，A′是A的轴射柱面投影，其柱面坐标为(θA′,ρA′,zA′)。投影公式为：

　　利用柱面坐标与直角坐标的转换公式，以直角坐标为媒介，将球坐标转换为柱面坐标。以柱面底圆的弧长为横轴，z轴为纵轴，沿着0°方位角处的母线将柱面展开成投影平面。重合部位采用加权融合，先同仰角层内融合，再层间融合，60°仰角仍需二次融合。

　　全天空图加权融合实验结果如图9(a)所示，完整呈现了低仰角图像信息，云量无失真地得以保留，但顶部（仰角90°附近）图像拉伸产生了空洞。采用渐入渐出加权平均插值算法利用空洞周边像素对空洞插值。实验结果如图9(b)所示，空洞消除且图像呈连续自然变化。

5结论

　　本文通过空间坐标转换和地图投影方法，用垂直投影和等积柱面投影实现了全天空云图拼接，采用渐入渐出加权平均融合算法解决了边界融合。该方法成本低廉，能够保留子图的信息量，解决了低仰角图像信息严重损失的弊端，为识别云状、计算云量提供了参考。提高电机角速度可缩短子图时差，根据经纬度和时间计算位置，再插值可去除太阳，结合其他波段或双目相机进行联合观测可进一步优化。

　　致谢：感谢中国气象局大气探测中心马舒庆研究员、陶法博士等对本课题的支持。

参考文献

　　［1］ 高太长，刘磊，赵世军,等.全天空测云技术现状及进展［J］.应用气象学报，2010，21(1):101108.

　　［2］ 孙晓钢，孙学金，牛珍聪，等．全天空云图获取的一种方式及算法实现［J］．气象科学，2008，28(3):338341.

　　［3］ 唐烁,缪源.基于Harris角点的图像匹配算法［J］.微型机与应用,2013,32(2):4143.

　　［4］ 李晓辉,周荫清,王祖林，等.基于曲面拼接的球面全景生成算法［J］.北京航空航天大学学报,2007,33(6):668681.

　　［5］ 焦健,曾琪明.地图学［M］.北京：北京大学出版社，2005.

　　［6］ 郑永安,宋建社，曹继平，等.基于小波变换的区域特征加权自适应图像融合算法［J］.计算机工程与应用，2004，28(6):9495，127.

　　［7］ 杨道富,彭明.球面和圆柱面相等对应关系的研究［J］.商丘师范学院学报,2005,21(5):6163.