謝 鵬，汪超亮，李子揚

（1.中國科學院 定量遙感信息技術重點實驗室，北京 100094；2.中國科學院 光電研究院，北京 100094；3.中國科學院大學，北京 100049）

基于陰影圖的航空遙感地面覆蓋快速分析方法

謝 鵬1,2,3，汪超亮1,2*，李子揚1,2

（1.中國科學院 定量遙感信息技術重點實驗室，北京 100094；2.中國科學院 光電研究院，北京 100094；3.中國科學院大學，北京 100049）

航空遙感系統(tǒng)在地形復雜地區(qū)側視成像作業(yè)時，有可能會因地形遮擋而出現成像盲區(qū)，在任務規(guī)劃階段進行載荷成像盲區(qū)檢測是解決該問題的有效方法。鑒于現有基于地形可視域分析的盲區(qū)檢測方法存在計算效率低、通用性差等缺點，提出了一種基于陰影圖算法的航空遙感地面覆蓋快速分析方法。實驗結果表明，該方法不僅具有效率高、通用性強的特點，而且可定量分析載荷成像地面覆蓋率，為航空遙感作業(yè)任務規(guī)劃的盲區(qū)檢測與覆蓋評估提供了很好的支持。

陰影圖算法；載荷成像地面覆蓋分析；航跡評估；三維可視化仿真

航空遙感作為一種重要的高分辨率遙感數據獲取手段，在資源探查、災情評估、軍事偵察等方面得到廣泛的應用。但航空遙感系統(tǒng)在地形復雜地區(qū)側視成像作業(yè)時，可能會因地形遮擋而出現成像盲區(qū)，因此在航空遙感作業(yè)任務規(guī)劃階段，需要先進行載荷成像的盲區(qū)檢測，再采取針對性的措施保證數據獲取的完整性，以提高航空遙感作業(yè)任務的成功率。

地形視域分析方法是目前航空遙感任務規(guī)劃階段載荷成像盲區(qū)檢測的重要手段，其為復雜地形環(huán)境下航空遙感任務地面覆蓋分析提供了很好的技術支持。為提高地形視域分析的效率，國內外學者開展了大量的研究工作?，F有的視域分析方法大多以LOS（line of sight）算法[1]為基礎，該算法利用點與點的幾何關系判斷地形可視性，通過比較視線上相應位置的高度和對應地形位置的高度判斷是否可視。高程內插方法是影響LOS效率的關鍵所在，常見內插算法有R3算法[2-4]、R2算法[4]、Xdraw[4]算法等。DYNTACS算法[5]采用非均等劃分視線段的雙線性插值算法進行視域分析，王智杰等對DYNTACS算法進行了改進[6]，根據視線的傾斜方向，只計算與橫向格網或與縱向格網的交點，提高了算法效率。應申等基于DEM數據的特征，提出了一種雙增量地形可視域[7]分析方法。此外，還有一些研究采用并行計算方法進行可視計算，如基于參考面的可視域算法[8]，通過點與點形成的空間關系判斷點與所有目標點是否可視，很大程度上提高了視域的計算效率。上述可視分析算法均需進行地形插值計算或格網的相交計算，計算效率不高，當地形范圍和分辨率提高時，會產生過量計算的問題，難以實現地形可視域的快速分析。

現有的地形視域分析方法大多基于CPU實現，在分析地形復雜、覆蓋范圍較大的區(qū)域或高分辨率地形數據時，計算冗余、效率低下，無法直接應用到航空遙感成像地面覆蓋快速分析。本文為滿足航空遙感成像盲區(qū)快速檢測的迫切需求，通過研究三維渲染管線和陰影圖技術[9]，提出一種基于陰影圖的地面覆蓋快速分析方法，將大量圖形計算任務通過GPU分析執(zhí)行，充分發(fā)揮GPU三維渲染能力，實現載荷在復雜三維場景的覆蓋快速分析，大大提高了視域分析的效率和三維可視化的顯示效果。

1 航空遙感地面覆蓋快速分析

1.1 載荷仿真成像深度圖構建

航空遙感載荷可分為被動型和主動型。被動型載荷接收成像視場內的地物反射或發(fā)射的電磁波信號；主動型載荷自身發(fā)射電磁波，然后接收由地物反射的電磁波信號。不論是哪種類型的遙感載荷，在側視成像時都有可能受地形遮擋的影響，導致電磁信號被遮擋而無法成像，形成載荷的成像盲區(qū)。在計算機圖形學中，通過光源模型進行場景陰影渲染時也有類似的過程，光線由光源主動發(fā)射，通過與場景地形計算獲取到陰影區(qū)，最終在屏幕上進行顯示。兩種類型載荷成像時的地面覆蓋計算均可通過光源模型的地形陰影渲染過程進行仿真計算，這里采用的光源模型為聚光燈模型。聚光燈模型作為一種位置性光源，按照一定的位置和方向進行場景的光照陰影處理，在進行三維地形場景渲染過程中，通過載荷成像時的位置和姿態(tài)計算得到聚光燈模型的位置和方向參數。三維渲染管線通過聚光燈模型按照攝像機成像模型[9]仿真載荷成像地面覆蓋情況，并將成像視場內的三維場景信息投射到二維的成像平面，其對應的瞬時視場內地形覆蓋仿真計算借鑒陰影圖[10]算法。

利用三維渲染的方法進行載荷的盲區(qū)檢測，首先需要使用載荷的成像參數構建三維渲染場景的成像深度圖，深度圖構建過程涉及不同的坐標系變換，其對應的變換公式為：

式中，Pv'為變換后的圖像紋理坐標；Mr為歸一化矩陣；Mproj為投影矩陣；Mper為視圖矩陣；Pv為變換前的世界坐標。Mr矩陣將復合變換結果中所有的分量取值范圍限定為[0,1]，使之適合紋理坐標索引的合理范圍。文中以面陣載荷為例，通過分析圖1所示的載荷成像模型，構建出其成像所涉及的變換矩陣。

假定面陣光學載荷系統(tǒng)的焦距為f，底片的高度為H，可根據載荷成像的幾何模型（圖1）計算瞬時視場角，其計算方法為：

面陣載荷成像視錐體投影矩陣Mproj描述了載荷成像視錐的大小，可用視場角fovy和成像視錐體寬高比aspect表示， 即

載荷的成像位置和姿態(tài)對載荷成像地面覆蓋的影響可通過視圖矩陣近似表達。假定載荷成像位置為eye，載荷成像中心位置為center，以及垂直于成像相機向上的向量up，3個參數唯一確定了成像相機的位置和姿態(tài)。在構建視圖矩陣之前通過這3個矢量構建相機的成像坐標系，即

然后在相機成像坐標系的基礎上構建出載荷成像的視圖矩陣，即

通過上述投影變換矩陣和視圖變換矩陣，可近似模擬面陣光學載荷的成像模型，同時結合式（1）、（2），渲染得到載荷在當前成像狀態(tài)下的深度圖（圖2）。

圖1 載荷成像三維幾何模型

圖2 模擬成像深度圖

1.2 載荷成像可視區(qū)域三維渲染

構建出載荷成像瞬時視場內三維場景的深度圖后，再根據深度信息判斷場景區(qū)域是否可見。具體方法為：首先通過將當前視點下的像素轉換到光源坐標系下，計算每個像素與光源的距離；然后將距離值和陰影圖中對應的深度值進行比較，確定當前像素是否處于陰影中；最后根據比較的結果對載荷成像盲區(qū)和載荷成像區(qū)域分別進行著色處理，從而實現對載荷成像盲區(qū)的檢測。

三維渲染管線對載荷成像區(qū)域進行三維場景渲染著色時，需區(qū)分成像可視區(qū)域、成像盲區(qū)和成像背景區(qū)域，文中采用3種不同的顏色進行著色處理，如綠色表示成像可視區(qū)域、紅色表示成像盲區(qū)，藍色表示成像背景區(qū)域。通過著色渲染后，不但能區(qū)分地形區(qū)域可見與否，而且能區(qū)分載荷的成像范圍和背景范圍，便于后續(xù)的定量分析。載荷成像盲區(qū)檢測渲染所對應的偽代碼如下：

算法名稱：載荷成像盲區(qū)檢測渲染算法

輸入：三維場景的頂點，深度紋理，紋理坐標gl_ TexCoord；

Fragment Shader:

{

Vec4 shadowcolor，visiblecolor，outscopecolor； //定義3種不同著色顏色；

//深度比較，當前片元是否可見，結果保存在compareResult中；

float compareResult = shadow2DProj( osgShadow_ shadowTexture, gl_TexCoord[1] )；

//歸一化坐標；

Vec3 outXYZ=vec3(gl_TexCoord[1].x/gl_TexCoord[1]. wgl_TexCoord[1].y/gl_TexCoord[1].w，gl_TexCoord[1.z]/ gl_TexCoord[1].w)；

//判斷位置進行著色；

If((outXYZ.x＞1)||(outXYZ.x＜0)||(outXYZ. y＞1)||(outXYZ.y＜0) ) //表明不在成像范圍；

{

gl_FragColor=outscopecolor；

}

Else if((outXYZ.x＞0)&&(outXYZ.x＜1)&(outXYZ. y＞0)&(outXYZ.y＜1)) //表明在成像范圍；

{

//判斷是否可視，并著色；

gl_FragColor=mix(shadowColor,visibleColor,compar eResult)；

}

}

1.3 載荷成像覆蓋參考面構建

雖然通過三維渲染管線的渲染區(qū)分了載荷的地面覆蓋區(qū)域內的成像范圍和成像盲區(qū)，但由于渲染管線作為一個狀態(tài)機，保持了高度封裝的特性，未提供訪問每一塊渲染區(qū)域的接口，無法直接獲得渲染區(qū)域的坐標信息。因此本文在場景中增加了一個額外的“虛擬”相機，用于記錄渲染管線的渲染結果，為下一步的定量分析提供二維柵格數據。此相機采用平行投影的投影方式，垂直于地面向下投影，如圖3所示。采用平行投影的原因有兩方面：在平行投影方式下，成像視錐是一個平行的長方體，可最大范圍覆蓋成像區(qū)域，不會產生成像盲區(qū)；在平行投影方式下，視口圖像不會發(fā)生扭曲變形。

圖3 “虛擬相機”的平行投影方式

1.4 載荷成像覆蓋率計算

覆蓋率是指載荷在瞬時視場內可觀測到的地面面積占載荷成像投影于參考面所涵蓋范圍面積的百分比。覆蓋率計算的可以有兩種定義方式：①將載荷成像區(qū)域的地形作為參考面，載荷當前位置下投影于地形參考面的成像外邊界范圍作為覆蓋率計算的總面積，如圖4a所示黃色折線圍起來的綠色區(qū)域和紅色區(qū)域；②采用成像區(qū)域平均海拔高度平面作為參考面，將載荷當前位置下投影于海拔平均高度參考面上的成像范圍作為覆蓋率計算的總面積，如圖4b所示黃色折線圈起來的綠色區(qū)域。在實際應用中，根據需求選擇不同的覆蓋率計算方法。

載荷成像覆蓋率計算是通過分析§1．2三維渲染后的二維柵格圖像獲取的。假定柵格圖像的寬度與高度分別為width和height，共有width×height個像元。設二維柵格圖像中紅色像素數目為a（表示載荷成像盲區(qū)的不可見像素數目），綠色像素數目為b（表示載荷成像可視區(qū)域的可見像素數目），藍色的像素數目為c（表示不在成像范圍的像素數目），a+b+c=width×height；將載荷在成像區(qū)域平均海拔高度平面上成像的可見像素數目記為d。根據覆蓋率的兩種定義，其對應的計算公式分別為：

式中，Φa表示第一種參考面定義下的載荷成像地面覆蓋率；Φb表示第二種參考面定義下的載荷成像地面覆蓋率。

圖4 參考面的三維渲染結果

2 實驗與結果分析

2.1 實驗環(huán)境

實驗計算機使用Windows 7操作系統(tǒng)，硬件配置為CPU i3 3．4GHz、顯卡NVIDIA GeForce GT620 M、內存4 GB，開發(fā)環(huán)境OpenSceneGraph 3．0．1，開發(fā)工具Visual Studio 2010。實驗采用內蒙古包頭地區(qū)的DEM數據：經度范圍108．775 6～109．203 8，緯度范圍40．618 4～40．766 0。

2.2 實驗結果

實驗采用面陣CCD載荷參數：焦距34．51 mm、底片的元件長7．1 um。成像位置為（42．708 2N，109．100 7E，3 000 m），其視場角計算公式為：

面陣CCD載荷45°側視角成像的效果如圖5所示。

圖5 側視成像效果圖

通過陰影圖算法的渲染分析，得到載荷在當前成像位置下的仿真計算結果，如圖6所示。其中，a圖像為視域分析的結果，紅色區(qū)域為載荷成像盲區(qū)，綠色區(qū)域為載荷成像可視區(qū)域，b圖像為左側對應方框的局部放大圖。

圖6 視域分析結果

通過本文方法對載荷成像地形區(qū)域的三維渲染，得到載荷地面覆蓋二維柵格圖像，如圖7所示。其中，a圖為實際地形參考面的成像渲染結果，b圖為采用平均海拔高度為1 270 m的參考面的渲染結果。二維柵格圖像相對應的像素統(tǒng)計信息如表1所示。

圖7 載荷成像渲染結果

表1 像素數目統(tǒng)計表

按照第一種參考面的定義，計算得到當前載荷瞬時成像覆蓋率為：

按照第二種參考面的定義，計算得到當前載荷瞬時成像覆蓋率為：

2.3 算法效率對比實驗

為了驗證本文方法的計算效率，將本文提出的方法分別與傳統(tǒng)的R3視域分析算法、商業(yè)軟件ArcGIS10．0的可視性分析模塊的可視域分析效率進行了比較。實驗隨機抽取不同視點，采用大小為325×224、650×507、1 446×520的高程數據進行了比對，3種算法的處理時間如表2所示。

表 2 不同分辨率的DEM數據下算法的處理時間/s

由表2可知，隨著地形數據精細程度的增加，R3算法和ArcGIS分析算法的計算時間會越來越長，本文算法在計算效率上有明顯優(yōu)勢，且對DEM分辨率不敏感，處理時間比較穩(wěn)定。

由于在實際應用中需考慮盲區(qū)分析的實時性，因此本次實驗同時也分析了深度圖分辨率對算法效率的影響，分別分析了512×512、1 024×1 024、2 048×2 048、4 096×4 096四種分辨率下的深度圖算法效率。

從圖8可知，隨著深度圖分辨率的增加，算法的效率呈明顯降低的趨勢，當深度圖分辨率達到4 096×4 096或更高時，很難進行實時的覆蓋分析。所以在進行載荷盲區(qū)實時檢測時，需根據需求和計算機配置條件選擇合適分辨率的深度圖，在效率與精度方面進行折衷考慮。由于實時進行場景三維渲染分析的幀率至少為24 Hz，文中實驗計算機配置條件下選擇1 024×1 024的深度圖分辨率較為合適。

圖8 深度圖分辨率對計算效率影響

3 結 語

本文借鑒了陰影圖算法，將其成功應用于航空遙感成像覆蓋快速分析，實現了成像過程的盲區(qū)分析檢測，在保持一定精度的基礎上，實時顯示當前載荷位置下的成像盲區(qū)以及快速計算成像覆蓋率，提高了航空載荷盲區(qū)檢測分析的效率。在未來的工作中，可從以下3個方面進行深入研究：①在深度圖分辨率較低時，渲染結果會出現鋸齒形邊緣現象，可進行渲染顯示的抗鋸齒優(yōu)化；②結合載荷的飛行航跡，估算載荷在整條航跡上的成像覆蓋情況；③該方法可擴展應用到建筑規(guī)劃、景區(qū)觀測點選址、手機基站布設等領域，可根據具體應用領域的特點進行實時分析算法優(yōu)化。

[1] Chen C, Zhang L, Ma J, et al． Adaptive Multi-resolution Labeling in Virtual Landscapes [J]． International Journal of Geographical Information Science, 2010, 24(6)∶ 949-964

[2] Shapira A． Visibility and Terrain Labeling [J]． Master's Thesis, Rensselaer Polytechnic Institute, 1994, 8(1)∶ 15-17

[3] Franklin W R, Ray C K, Mehta S． Geometric Algorithms for Sitting of Air Defense Missile Batteries Research Project for Battelle, Columbus Division[J]． Delivery Order, 1994, 27(5)∶41-44

[4] 繆賁術,屈軍亮． 基于Xdraw 算法的SAR 航線規(guī)劃[J]． 艦船電子工程,2011,31(10)∶ 52-55

[5] De Floriani L, Marzano P, Puppo E． Line-of-sight Communication on Terrain Models [J]． International Journal of Geographical Information Systems, 1994, 8(4)∶ 329-342

[6] 王智杰,印欣,朱良家． 地形通視性算法改進[J]． 情報指揮控制系統(tǒng)與仿真技術, 2004 (2)∶ 54-57

[7] 應申,李霖,梅洋,等．增量法地形可視計算與分析[J]．測繪學報, 2007,36(2)∶192-197

[8] Zhi Y, Wu L, Sui Z, et al． An Improved Algorithm for Computing Viewshed Based on Reference Planes [C]．Geoinformatics, 2011 19th International Conference on IEEE, 2011

[9] 蘇國中,鄭順義,張劍清,等． OpenGL 模擬攝影測量方法研究[J]．中國圖像圖形學報,2006,11(4)∶ 540-544

[10] Williams L． Casting Curved Shadows on Curved Surfaces[C]．ACM Siggraph Computer Graphics, ACM, 1978

6 結 語

在本檢定實例中，電子羅盤顯示的定向角度值為109°30'16"，示值誤差小于最小允許誤差的15'，遠遠大于本校準方法的擴展不確定度，對于GPS電子羅盤，是一種可靠的校準方法。

參 考 文 獻

[1] 田湘,范勝林,劉建業(yè)．基于GPS的短基線姿態(tài)系統(tǒng)研究[J]．傳感器與微系統(tǒng),2007,26(10)∶29-31

[2] 王建斌．GPS方位角及其在導線測量中的應用研究[J]．北京測繪,2005(3)∶39-42

[3] 熊建明,王宇飛．GPS短邊方位測量若干問題的探討[J]．四川測繪,2002,25(3)∶117-118

[4] 熊建明．GPS短邊方位測量的精度分析[J]．測繪通報,2000(11)∶18-20

[5] 王解先,劉慧芹,唐立軍．不同站心地平坐標系下的坐標歸算[J]．工程勘察,2005(5)∶58-59

[6] 孔祥元,郭際明,劉宗泉．大地測量學基礎[M]．武漢∶武漢大學出版社,2001

[7] 曹金國,廖望東,李建偉,等．GPS短邊方位角測量精度定量分析研究[J]．全球定位系統(tǒng),2013,38(1)∶18-19

第一作者簡介：彭友志，工程師，主要從事GPS和全站儀等測繪儀器的計量和相關科研工作。

P237

B

1672-4623（2016）01-0071-05

10.3969/j.issn.1672-4623.2016.01.021

謝鵬，碩士，主要從事遙感信息處理的研究。

2014-10-21。

項目來源：國家高技術研究發(fā)展計劃資助項目（2013AA122105）。（*為通訊作者）