陳青青1,2，湯志亞1，楊玲1,2

　?。?.成都信息工程大學(xué)，四川 成都 610225；2.中國(guó)氣象局 大氣探測(cè)重點(diǎn)開放實(shí)驗(yàn)室，四川 成都 610225）

       摘要：采用地基鏡頭旋轉(zhuǎn)掃描拍攝不同俯仰角的天空?qǐng)D像，通過全景圖拼接獲得全天空圖像，避免了魚眼透視畸變，保留了低仰角信息。目前全景圖拼接技術(shù)主要局限于單層柱面，通過空間坐標(biāo)轉(zhuǎn)換及以漸入漸出加權(quán)平均融合算法和插值算法，借鑒地圖學(xué)投影原理，將28張圖像在垂直投影面和等積柱面投影面進(jìn)行無縫拼接得到360°×180°的全天空全景圖，為識(shí)別云狀、等積計(jì)算云量提供了連續(xù)完整的全天空?qǐng)D像。

　　關(guān)鍵詞：地圖投影；全天空；云圖；全景圖拼接

0引言

　　圖1原始圖像隨著大氣探測(cè)技術(shù)的發(fā)展，云的觀測(cè)逐漸自動(dòng)化，國(guó)內(nèi)外研發(fā)了大量不同波段不同視角的測(cè)云儀器［1］，全景拍攝方式簡(jiǎn)單快捷但低仰角畸變大；另外還有全景圖拼接方式［2］。本文采用監(jiān)控?cái)z像機(jī)，無盲區(qū)掃描拍攝可見光波段天空?qǐng)D像，得到36張具有重疊部分的子圖，通過還原物體在天球中的相對(duì)位置，并借鑒地圖投影轉(zhuǎn)換為平面圖像，對(duì)重疊部分進(jìn)行融合，可得到高清連續(xù)的全天空云圖。

1觀測(cè)平臺(tái)

　　云臺(tái)型號(hào)為HY900B，相機(jī)為索尼FCBEX1010P，采集器為HB18204T，圖像傳感器CCD為EXview HAD 1/4英寸，成像面積為3.2 mm×2.4 mm，固定焦距為3.4 mm，視角為57.8°，圖像尺寸為704×576，長(zhǎng)邊視角為44.70°，寬邊視角為36.58°，圖像上下左右都存在重疊，無盲區(qū)。實(shí)驗(yàn)用機(jī)配置Intel Core i3 CPU，2 GB內(nèi)存，Windows XP系統(tǒng)。原始圖像如圖1所示，拍攝俯仰角順序如表1所示。

　　令拍攝次序t=1~36，其中1~9表示相機(jī)光軸仰角為0°的圖像，10~18表示仰角30°圖像……由此可見第t張圖像方位角θt仰角φt與t存在如式(1)所示的關(guān)系，其中［·］表示向下取整，%表示求余。

　　相機(jī)標(biāo)定以標(biāo)準(zhǔn)網(wǎng)格角點(diǎn)作為匹配點(diǎn)進(jìn)行多項(xiàng)式擬合［3］，將擬合關(guān)系式應(yīng)用到整幅圖像并對(duì)坐標(biāo)進(jìn)行雙線性插值，截去無效邊緣，如圖2所示。

　　圖2相機(jī)標(biāo)定

　　2圖像在球面上的配準(zhǔn)

　　圖像拼接包括配準(zhǔn)和融合兩個(gè)環(huán)節(jié)。本文圖像有確定的空間坐標(biāo)，能實(shí)現(xiàn)圖像在空間位置上的鑲嵌，從而找到配準(zhǔn)點(diǎn)。

　　以相機(jī)位置為原點(diǎn)，拍第一幅圖像時(shí)相機(jī)光軸方向?yàn)?°方位角，觀測(cè)點(diǎn)所在的大地水平面為xoy面，0°方位角為x軸的正向，90°方位角為y軸正向，過原點(diǎn)xoy面指向天空的垂線為z軸。0°方位角圖像在空間xoz平面上的剖面圖如圖3(a)所示。其中E′表示相機(jī)光軸投影點(diǎn)E點(diǎn)在xoz面上的投影，在xoz平面上的坐標(biāo)為(xE,zE),由圖3(a)可見。

　　其中，L表示物距，P為圖像中任意一點(diǎn)在xoz面上的投影，如圖3(a)所示，P在x軸上的坐標(biāo)為xP=xE+Δx,P在z軸上的坐標(biāo)為zP=zE－Δz，而

　　由此得出兩個(gè)關(guān)系式：

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

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

　　聯(lián)立式（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)

　　數(shù)字圖像坐標(biāo)是離散值，直接由（j，i，t）計(jì)算（x，y，z），會(huì)導(dǎo)致有些空間坐標(biāo)找不到與之對(duì)應(yīng)的圖像坐標(biāo)，投影后圖像出現(xiàn)空洞。采用反算插值的方法可避免出現(xiàn)這種情況，即根據(jù)投影面上每個(gè)點(diǎn)來反查其在原圖中的位置和像素值［4］。

　　當(dāng)相機(jī)光軸θt=0°，φt≠0°時(shí)，聯(lián)立式（6）、式（8）反解j、i得：

　　聯(lián)立式（7）、式（9）即得到θt=0°，φt≠0°時(shí)球面各點(diǎn)（x，y，z）所對(duì)應(yīng)的第t張圖的圖像坐標(biāo)（j，i，t）。

　　同理，當(dāng)θt=0°、φt=0°時(shí)，聯(lián)立公式（5）、式（8），反解出

　　聯(lián)立式（6）、式（10）即得到θt=0°、φt=0°時(shí)球面各點(diǎn)對(duì)應(yīng)的第t張圖的圖像坐標(biāo)。

　　探測(cè)范圍方位角0°~360°，仰角0°~180°，半徑0~∞,該范圍的物體可投影到以相機(jī)為球心，給定值為半徑的半球面上。假想景物已從空間投影到了天球面上，相機(jī)拍攝將天球物體投影到CCD，因此應(yīng)采用點(diǎn)光源在球心的球面投影方式［5］。

　　設(shè)P(xP,yP,zP)是方位角0°，仰角φt的圖像上任意一點(diǎn)，利用直角坐標(biāo)與球坐標(biāo)的轉(zhuǎn)換公式得到P點(diǎn)的球坐標(biāo)為(θp,φp,rp)。

　　P點(diǎn)與投影點(diǎn)P′位于同一條直線上，二者的θ和φ相同，而所有投影點(diǎn)的r都等于投影球面半徑R，因此P點(diǎn)的投影坐標(biāo)為P′(θp,φp,R)。

　　云臺(tái)順時(shí)針旋轉(zhuǎn)，那么球面投影坐標(biāo)（θp－θt，φp，R）對(duì)應(yīng)第t張圖像的圖像坐標(biāo)（j，i）。如式（1）所示，改變?chǔ)萾就能通過方位角0°的圖像投影點(diǎn)求得其余8個(gè)方位對(duì)應(yīng)的圖像投影坐標(biāo)點(diǎn)，從而得到360°全景圖。仰角90°時(shí)只取一張圖就能覆蓋拍攝范圍，因此僅需用到前28張圖。為避免信息冗余和損失，本文取物距與長(zhǎng)邊弦長(zhǎng)的中值為投影球面半徑。

3垂直投影面的拼接

　　垂直投影光源為一束與z軸平行的平行光，投影面為xoy面，投影公式如式（11）所示，（x，y，z）表示空間某點(diǎn)，（x′，y′，z′）為其投影點(diǎn)。

　　圖4全天空云圖垂直投影實(shí)驗(yàn)結(jié)果MATLAB仿真垂直投影結(jié)果如圖4所示，銜接無誤位置準(zhǔn)確，但邊界亮度躍變，這是相機(jī)自動(dòng)曝光引起的。數(shù)字圖像以矩陣形式存儲(chǔ)，反算下標(biāo)會(huì)取整而產(chǎn)生誤差，應(yīng)盡量減少取整次數(shù)，在最后一級(jí)投影才對(duì)下標(biāo)取整和圖像融合。

　　加權(quán)平均融合是指在匹配點(diǎn)重合的n個(gè)子圖各自對(duì)應(yīng)的3個(gè)顏色通道的亮度值進(jìn)行加權(quán)平均，新值賦給匹配點(diǎn)，當(dāng)權(quán)值相等時(shí)稱為均值融合，均值融合結(jié)果如圖5所示。

　　結(jié)果表明均值融合后圖像個(gè)體完整，但邊界明顯，原因是沒有考慮權(quán)值差異，距離匹配點(diǎn)近的子圖貢獻(xiàn)權(quán)值應(yīng)更大。以距離比為權(quán)值，將距離近的圖像在該點(diǎn)的值作為主元素進(jìn)行融合，稱為漸入漸出加權(quán)平均融合（后面簡(jiǎn)稱加權(quán)融合）［6］。

　　子圖垂直投影之后的圖像是扇形的，由于圖像是圍繞圓心旋轉(zhuǎn)的，權(quán)值便不應(yīng)用直線距離，而是該點(diǎn)到邊界的弧長(zhǎng)。對(duì)30°層進(jìn)行加權(quán)融合實(shí)驗(yàn)結(jié)果如圖6所示，可見該方法解決了子圖邊緣突變，效果良好。

　　對(duì)60°層進(jìn)行一次加權(quán)融合，僅顯示其重疊部分如圖7(a)所示，可見重疊部分還有重疊，再次加權(quán)融合，結(jié)果如圖7(b)所示，重疊消失。最后將一級(jí)融合后的像元和無重疊像元放入，形成60°層全景如圖7(c)所示。

　　先進(jìn)行同仰角層內(nèi)融合，再進(jìn)行層間融合。層間融合用半徑的變化量來計(jì)算權(quán)值，圖8（a）是30°、60°層間未融合圖像，圖8（b）是30°、60°層間加權(quán)融合，圖8（c）是融合后的全天空?qǐng)D像。

　　4等積柱面投影面的拼接

　　云量是按視面積來計(jì)算的，等積投影［5］為云量的計(jì)算提供了借鑒。前人證明軸射柱面投影下，任意形狀的球面與其在外切圓柱面上投影面積相等［7］。設(shè)A是球面上任意一點(diǎn)，A的柱面坐標(biāo)為(θA,ρA,zA)，A′是A的軸射柱面投影，其柱面坐標(biāo)為(θA′,ρA′,zA′)。投影公式為：

　　利用柱面坐標(biāo)與直角坐標(biāo)的轉(zhuǎn)換公式，以直角坐標(biāo)為媒介，將球坐標(biāo)轉(zhuǎn)換為柱面坐標(biāo)。以柱面底圓的弧長(zhǎng)為橫軸，z軸為縱軸，沿著0°方位角處的母線將柱面展開成投影平面。重合部位采用加權(quán)融合，先同仰角層內(nèi)融合，再層間融合，60°仰角仍需二次融合。

　　全天空?qǐng)D加權(quán)融合實(shí)驗(yàn)結(jié)果如圖9(a)所示，完整呈現(xiàn)了低仰角圖像信息，云量無失真地得以保留，但頂部（仰角90°附近）圖像拉伸產(chǎn)生了空洞。采用漸入漸出加權(quán)平均插值算法利用空洞周邊像素對(duì)空洞插值。實(shí)驗(yàn)結(jié)果如圖9(b)所示，空洞消除且圖像呈連續(xù)自然變化。

5結(jié)論

　　本文通過空間坐標(biāo)轉(zhuǎn)換和地圖投影方法，用垂直投影和等積柱面投影實(shí)現(xiàn)了全天空云圖拼接，采用漸入漸出加權(quán)平均融合算法解決了邊界融合。該方法成本低廉，能夠保留子圖的信息量，解決了低仰角圖像信息嚴(yán)重?fù)p失的弊端，為識(shí)別云狀、計(jì)算云量提供了參考。提高電機(jī)角速度可縮短子圖時(shí)差，根據(jù)經(jīng)緯度和時(shí)間計(jì)算位置，再插值可去除太陽，結(jié)合其他波段或雙目相機(jī)進(jìn)行聯(lián)合觀測(cè)可進(jìn)一步優(yōu)化。

　　致謝：感謝中國(guó)氣象局大氣探測(cè)中心馬舒慶研究員、陶法博士等對(duì)本課題的支持。

參考文獻(xiàn)

　?。?］ 高太長(zhǎng)，劉磊，趙世軍,等.全天空測(cè)云技術(shù)現(xiàn)狀及進(jìn)展［J］.應(yīng)用氣象學(xué)報(bào)，2010，21(1):101108.

　　［2］ 孫曉鋼，孫學(xué)金，牛珍聰，等．全天空云圖獲取的一種方式及算法實(shí)現(xiàn)［J］．氣象科學(xué)，2008，28(3):338341.

　?。?］ 唐爍,繆源.基于Harris角點(diǎn)的圖像匹配算法［J］.微型機(jī)與應(yīng)用,2013,32(2):4143.

　　［4］ 李曉輝,周蔭清,王祖林，等.基于曲面拼接的球面全景生成算法［J］.北京航空航天大學(xué)學(xué)報(bào),2007,33(6):668681.

　?。?］ 焦健,曾琪明.地圖學(xué)［M］.北京：北京大學(xué)出版社，2005.

　　［6］ 鄭永安,宋建社，曹繼平，等.基于小波變換的區(qū)域特征加權(quán)自適應(yīng)圖像融合算法［J］.計(jì)算機(jī)工程與應(yīng)用，2004，28(6):9495，127.

　?。?］ 楊道富,彭明.球面和圓柱面相等對(duì)應(yīng)關(guān)系的研究［J］.商丘師范學(xué)院學(xué)報(bào),2005,21(5):6163.