陳青青1,2，湯志亞1，楊玲1,2
?。?.成都信息工程大學，四川 成都 610225；2.中國氣象局 大氣探測重點開放實驗室，四川 成都 610225）
摘要：采用地基鏡頭旋轉掃描拍攝不同俯仰角的天空圖像，通過全景圖拼接獲得全天空圖像，避免了魚眼透視畸變，保留了低仰角信息。目前全景圖拼接技術主要局限于單層柱面，通過空間坐標轉換及以漸入漸出加權平均融合算法和插值算法，借鑒地圖學投影原理，將28張圖像在垂直投影面和等積柱面投影面進行無縫拼接得到360°×180°的全天空全景圖，為識別云狀、等積計算云量提供了連續(xù)完整的全天空圖像。
0引言
圖1原始圖像隨著大氣探測技術的發(fā)展，云的觀測逐漸自動化，國內外研發(fā)了大量不同波段不同視角的測云儀器［1］，全景拍攝方式簡單快捷但低仰角畸變大；另外還有全景圖拼接方式［2］。本文采用監(jiān)控攝像機，無盲區(qū)掃描拍攝可見光波段天空圖像，得到36張具有重疊部分的子圖，通過還原物體在天球中的相對位置，并借鑒地圖投影轉換為平面圖像，對重疊部分進行融合，可得到高清連續(xù)的全天空云圖。
1觀測平臺
云臺型號為HY900B，相機為索尼FCBEX1010P，采集器為HB18204T，圖像傳感器CCD為EXview HAD 1/4英寸，成像面積為3.2 mm×2.4 mm，固定焦距為3.4 mm，視角為57.8°，圖像尺寸為704×576，長邊視角為44.70°，寬邊視角為36.58°，圖像上下左右都存在重疊，無盲區(qū)。實驗用機配置Intel Core i3 CPU，2 GB內存，Windows XP系統(tǒng)。原始圖像如圖1所示，拍攝俯仰角順序如表1所示。
令拍攝次序t=1~36，其中1~9表示相機光軸仰角為0°的圖像，10~18表示仰角30°圖像……由此可見第t張圖像方位角θt仰角φt與t存在如式(1)所示的關系，其中［·］表示向下取整，%表示求余。
相機標定以標準網格角點作為匹配點進行多項式擬合［3］，將擬合關系式應用到整幅圖像并對坐標進行雙線性插值，截去無效邊緣，如圖2所示。
圖2相機標定
2圖像在球面上的配準
圖像拼接包括配準和融合兩個環(huán)節(jié)。本文圖像有確定的空間坐標，能實現(xiàn)圖像在空間位置上的鑲嵌，從而找到配準點。
以相機位置為原點，拍第一幅圖像時相機光軸方向為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)
聯(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ù)字圖像坐標是離散值，直接由（j，i，t）計算（x，y，z），會導致有些空間坐標找不到與之對應的圖像坐標，投影后圖像出現(xiàn)空洞。采用反算插值的方法可避免出現(xiàn)這種情況，即根據投影面上每個點來反查其在原圖中的位置和像素值［4］。
當相機光軸θt=0°，φt≠0°時，聯(lián)立式（6）、式（8）反解j、i得：
聯(lián)立式（7）、式（9）即得到θt=0°，φt≠0°時球面各點（x，y，z）所對應的第t張圖的圖像坐標（j，i，t）。
同理，當θt=0°、φt=0°時，聯(lián)立公式（5）、式（8），反解出
聯(lián)立式（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所示，銜接無誤位置準確，但邊界亮度躍變，這是相機自動曝光引起的。數(shù)字圖像以矩陣形式存儲，反算下標會取整而產生誤差，應盡量減少取整次數(shù)，在最后一級投影才對下標取整和圖像融合。
加權平均融合是指在匹配點重合的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)所示，完整呈現(xiàn)了低仰角圖像信息，云量無失真地得以保留，但頂部（仰角90°附近）圖像拉伸產生了空洞。采用漸入漸出加權平均插值算法利用空洞周邊像素對空洞插值。實驗結果如圖9(b)所示，空洞消除且圖像呈連續(xù)自然變化。
5結論
本文通過空間坐標轉換和地圖投影方法，用垂直投影和等積柱面投影實現(xiàn)了全天空云圖拼接，采用漸入漸出加權平均融合算法解決了邊界融合。該方法成本低廉，能夠保留子圖的信息量，解決了低仰角圖像信息嚴重損失的弊端，為識別云狀、計算云量提供了參考。提高電機角速度可縮短子圖時差，根據經緯度和時間計算位置，再插值可去除太陽，結合其他波段或雙目相機進行聯(lián)合觀測可進一步優(yōu)化。
致謝：感謝中國氣象局大氣探測中心馬舒慶研究員、陶法博士等對本課題的支持。
參考文獻
?。?］ 高太長，劉磊，趙世軍,等.全天空測云技術現(xiàn)狀及進展［J］.應用氣象學報，2010，21(1):101108.
?。?］ 孫曉鋼，孫學金，牛珍聰，等．全天空云圖獲取的一種方式及算法實現(xiàn)［J］．氣象科學，2008，28(3):338341.
?。?］ 唐爍,繆源.基于Harris角點的圖像匹配算法［J］.微型機與應用,2013,32(2):4143.
?。?］ 李曉輝,周蔭清,王祖林，等.基于曲面拼接的球面全景生成算法［J］.北京航空航天大學學報,2007,33(6):668681.
?。?］ 焦健,曾琪明.地圖學［M］.北京：北京大學出版社，2005.
［6］ 鄭永安,宋建社，曹繼平，等.基于小波變換的區(qū)域特征加權自適應圖像融合算法［J］.計算機工程與應用，2004，28(6):9495，127.
?。?］ 楊道富,彭明.球面和圓柱面相等對應關系的研究［J］.商丘師范學院學報,2005,21(5):6163.