鐘 成，李 卉，黃先鋒

（1.中國地質大學 教育部長江三峽庫區(qū)地質災害研究中心，湖北 武漢 430074；2.中國地質大學 行星科學研究所，湖北 武漢 430074；

3.武漢大學 測繪遙感信息工程國家重點實驗室，湖北 武漢 430079）

傳統(tǒng)正射糾正未考慮建筑物、樹木等地面目標的投影變形，在城市大比例尺正射影像上，建筑物傾斜遮擋的情況非常嚴重，影響了影像信息的有效解譯和利用.20世紀90年代以來，隨著世界范圍內城市三維仿真項目的興起和城市大比例尺地圖的應用需求日益強烈，利用城市三維模型制作“真正射影像”的研究和生產逐漸展開.

真正射影像制作的核心問題是高效穩(wěn)健的遮蔽檢測過程.根據攝影測量的基本原理，非正射投影造成的高層建筑傾斜會導致對低矮目標的遮蔽，若直接利用數(shù)字微分糾正方法為DSM（數(shù)字表面模型）格網尋找對應影像信息，低矮目標將獲得高層建筑的影像信息，其結果是造成高層建筑在正射影像上多次出現(xiàn)，學術界稱之為“鬼影”現(xiàn)象.這樣的正射影像無法談及幾何精度或光譜精度等問題.解決該問題的思路是，首先檢測出高層建筑對低矮目標的遮蔽情況，在正射糾正過程中這些遮蔽地區(qū)不分配影像信息，最后利用影像序列的重疊信息或利用仿真方法修復遮蔽地區(qū).

已有遮蔽檢測方法大多基于計算機圖形學的ZBuffer消隱技術，比較目標距光源點的距離判斷其遮蔽情況［1－4］.這類方法在處理大量復雜的矢量建筑物模型時效率很低［5］；而在處理柵格數(shù)據時，因難以捕獲相機畸變和投影成像造成的柵格變形，將導致錯誤的判斷競爭關系，造成大量偽遮蔽和偽可見噪聲［6－7］.此外，該方法無法處理被高窄墻體遮蔽目標的偽可見問題，即所謂“M-Portion”問題［3，8－9］，改進的方法需要將墻面納入遮蔽檢測計算過程，大大增加了計算工作量和復雜度.產生這些問題的根本原因在于復雜多樣的現(xiàn)實世界難以用簡單的三維體塊模型表示，基于圖形學的Z-Buffer技術也未顧及攝影測量的特殊情況［10］.

1 快速遮蔽檢測方法

從計算效率的角度出發(fā)，優(yōu)先考慮采用柵格DSM參與計算；從計算穩(wěn)健性的角度出發(fā)，需克服各類偽可見和偽遮蔽噪聲.為此可采用徑向投影高度角的變化情況進行遮蔽區(qū)域分析，其原理是，同一方位射線上逐漸遠離天底點時目標的投影高度角逐漸增大；在進入被遮蔽地區(qū)時投影高度角突然變小，直到離開被遮蔽地區(qū)時又回到之前的數(shù)值.據此，可以容易發(fā)現(xiàn)被遮蔽地區(qū).該方法對遮蔽地區(qū)和遮蔽源的掃描都是連續(xù)的，對遮蔽地區(qū)的判斷僅需高層建筑的邊緣柵格，可避免因柵格變形造成的錯誤競爭，也可避免因缺乏墻面信息造成的“M-Portion”問題［11］.

為全面檢測遮蔽情況，從天底點到DSM邊緣所有柵格的射線都需要分析太陽高度角的變化情況.此時，鄰近天底點的柵格將被多次訪問，導致效率低下.可行的方法包括徑向掃描和螺旋掃描兩種方法（圖1），前者首先以不同半徑將DSM劃分為多個圓環(huán)，在每個圓環(huán)內根據完整掃描的要求制定合適的掃描和方向間隔；后者從天底點出發(fā)以螺旋方式掃描DSM柵格，最終完成對DSM的完整掃描.本文力圖通過分析兩種方法的計算步驟和算法復雜度，推導最優(yōu)的配置方案，最終提出具有最高效率的遮蔽檢測方法.

兩種方法基本原理一致，區(qū)別在于對DSM掃描的方式，以及由此導致的不同的掃描面積、方位射線數(shù)和進行可見性判斷的面積.因此，通過分析計算過程的復雜度，找到最優(yōu)的分配方法，可找到復雜度最低的遮蔽檢測方法.首先，基于算法原理利用嚴密的幾何關系可以推出兩種方法不同的掃描面積、射線數(shù)據和可見性判斷面積；然后，利用計算數(shù)學理論推導出各個步驟包含的計算復雜度；最后利用微積分的極值理論求出具有最小復雜度的方法.

1.1 算法過程

首先，對于給定DSM 柵格（x，y，z），投影中心坐標為（Xc，Yc，Zc），計算并記錄每個DSM 柵格的投影高度角a.

圖1 徑向掃描和螺旋掃描方法Fig.1 Illustration of radial sweep and spiral sweep method

然后掃描DSM，分析沿掃描射線的投影高度角變化情況.對于徑向掃描，計算當前范圍內由邊緣柵格決定的方位射線，沿著射線方向逐個比較a變化，同時記錄柵格的可見性，若發(fā)現(xiàn)a突然變小，則記錄當前的amax值，投影高度角未超出amax的柵格都被認為是被遮蔽的地區(qū)，其可見性值設為零.

對于螺旋掃描方式，需要計算每個柵格與天底點確定的方位射線，然后比較其與射線上相鄰且靠近天底點的柵格的a變化，同樣記錄當前的amax值，投影高度角未超出amax的柵格被認為是被遮蔽的地區(qū)，其可見性值設為零.

由上述分析可知，兩種方法效率的影響因子包括不同的掃描面積、方位射線數(shù)和可見性判斷面積，需要研究這些因子的優(yōu)化方法，以提高算法效率.

1.2 計算復雜度分析

假定DSM為圓形且未進行分割，設定柵格尺寸為1，徑向掃描的每條掃描射線的面積為

式中：r為DSM半徑；d為DSM柵格的直徑.而掃描射線的數(shù)量可以由下式獲得：

可計算得到徑向掃描的掃描面積為

從式（4）可以看到，徑向掃描面積兩倍于DSM實際面積，存在較大程度的重復掃描.考慮到若一個柵格的可見性已經判定，那么在重復掃描的過程中，不需要再重復判斷，因此需要進行可見性判定的面積為πr2.

若DSM進行了n次分割，根據圖1，其掃描面積、射線數(shù)量和可見性判斷面積計算公式演變?yōu)?/p>

式中：Sr為徑向掃描面積；Lr為徑向掃描射線數(shù)；Vr為徑向掃描需要進行可見性判斷的面積.

對于螺旋掃描，容易判斷其掃描面積和可見性面積等于DSM的面積.考慮到每個柵格都需要計算方位射線，可以得到

式中：Ss為螺旋掃描面積；Ls為螺旋掃描射線數(shù)；Vs為螺旋掃描需要進行可見性判斷的面積.

通過計算數(shù)學分析，可知計算每一條射線參數(shù)需要5個基本操作（加減乘除），確定柵格的可見性需要4步操作.在徑向掃描中需要3步操作來掃描一個柵格，而螺旋掃描中僅需要1步可實現(xiàn)柵格定位.考慮到這些基本操作的復雜度是常量，兩種方法的計算復雜度可以估計為

式中：Cr代表徑向掃描方法的復雜度；Cs代表螺旋掃描方法的復雜度.

式（7）表明徑向掃描方法的復雜度主要依賴于DSM劃分方式，而螺旋掃描的復雜度依賴于DSM的大小.在給定DSM的情況下，螺旋掃描的復雜度是一個常量，而徑向掃描的復雜度是可以通過修改DSM劃分方法而改善的.

1.3 求解最小復雜度

兩種方法計算復雜度的差異可以通過下式計算：

若采用等差數(shù)列方法劃分DSM，即ri－ri－1＝m（m＞0）可得

在給定DSM的情況下，r為常量，此時求最小復雜度演變?yōu)榍蠼鈫巫兞亢瘮?shù)的最小值問題.考慮到n和ΔC在［1，r］是連續(xù)的，且ΔC是可微的，對ΔC求極值可得到

則可獲得最小ΔC的分割數(shù)為

帶入式（9），計算得到最小的ΔC為

當r大于10時，ΔC的值小于0，也就是說通過選擇最佳分割數(shù)，徑向掃描的復雜度此時小于螺旋掃描.這證明螺旋掃描效率高于徑向掃描的一般性結論不成立.

若采用等比數(shù)列方法分割 DSM，即ri／ri－1＝m（m＞0），式（8）演變?yōu)?/p>

由于求解冪指數(shù)的極值比較復雜，反而會影響計算復雜度，本文不再討論.同樣的，也可以采取其他方法劃分DSM，進而比較算法復雜度.

上述推導已經證明通過求解最優(yōu)分割數(shù)，徑向掃描方法的復雜度可以得到大幅度改善并獲得比螺旋掃描更高的計算效率.因此，基于最優(yōu)分割數(shù)的徑向掃描方法可認為是一種高效穩(wěn)健的遮蔽檢測方法.

2 實驗與討論

考慮到LiDAR（激光掃描數(shù)據）可以快速有效地獲取地表三維模型，已經日益成為城市三維建模的重要數(shù)據源，實驗采用LiDAR點云構建高分辨率DSM；作為基礎實驗數(shù)據，引入Z-Buffer方法對比算法的穩(wěn)健性和效率.

實驗數(shù)據如圖2所示，包含加拿大多倫多地區(qū)的LiDAR點云數(shù)據（圖2a）和航空影像數(shù)據（圖2b）.航空影像的分辨率為0.5m，LiDAR點云的分辨率為1.0m.航空影像記錄了城市目標的灰度和紋理信息，LiDAR數(shù)據記錄城市目標的三維坐標.在圖2b中左下部分有重大偏移和遮蔽的地區(qū)，從圖2c展示的局部點云可知，這些地區(qū)點云中不存在變形和遮蔽.因此，利用LiDAR構建的DSM實施航空影像的真正射糾正是可行的.

圖2 實驗數(shù)據Fig.2 Experimental data

圖3展示了利用典型方法實施遮蔽檢測的實驗結果.圖3a展示了Z-buffer方法產生的可見性矩陣，圖3b展示其局部放大圖，圖3c展示了對其進行Z-buffer遮蔽檢測的結果，圖3d展示了徑向掃描方法產生的可見性矩陣.從圖3a～3c中，可以發(fā)現(xiàn)大量的離散噪聲，以及由噪聲組成的黑色曲線.這些是因為Z-buffer方法無法適應DSM起伏和分辨率差異引起的.如果采用濾波方法消除噪聲，將降低影像的細節(jié)層次.而在圖3d～3f中則沒有這些離散噪聲和曲線，表明基于徑向投影高度角變化的檢測方法可以有效地適應城市地區(qū)復雜的DSM起伏和分辨率差異.圖3也表明，由于LiDAR數(shù)據分辨率低于遙感影像，由其構建的DSM生產的真正射影像中缺乏精確的斷裂線，建筑物輪廓出現(xiàn)鋸齒.

圖3 遮蔽檢測實驗Fig.3 Detecting occlusions with typical methods

高窄的結構容易產生M－portion問題，在數(shù)據源中較難發(fā)現(xiàn)此類結構，為了驗證算法對于 M－portion問題的處理能力，圖4展示了兩種方法對于高窄結構的檢測結果.圖4a展示了實驗中人為添加了模擬的高窄結構.圖4b是利用Z-buffer方法檢測的結果.可見建筑頂面的遮蔽范圍被檢測出來，而部分被墻面遮蔽地區(qū)未被檢測到.這是因為考慮到效率問題，一般未將墻面納入DSM參與計算和分析，若墻面遮蔽地區(qū)大于建筑本身面積，則會導致偽可見的情況，即所謂“M-portion”問題.高且窄的建筑物最容易導致這類問題的出現(xiàn).圖4c是利用徑向掃描方法檢測的結果，可見建筑物頂面和墻面遮蔽的地區(qū)都被檢測出來了.根據該方法的原理，僅需掃描建筑物邊緣的柵格即可完成對頂面和墻面遮蔽的檢測，無需額外改進和其他步驟.而絕大部分基于Z－buffer原理的方法及改進算法往往很難高效地發(fā)現(xiàn)并解決M－portion問題.從這一點來說，該方法極大地提高了遮蔽檢測的穩(wěn)健性.

圖4 對高窄結構的檢測Fig.4 Detecting occlusions with simulated DSM

本文對典型方法的計算效率做了對比實驗，包括Z-Buffer方法，徑向掃描方法、螺旋掃描方法和最優(yōu)分割數(shù)徑向掃描方法.實驗設備是具有雙核2.8G CPU和3.0GDDR內存的微機.圖5展示了效率實驗的結果.

圖5 效率實驗Fig.5 Performance of typical methods

從圖5可見，隨著給定DSM尺寸增大，各類方法消耗的時間增加.其中，對于任何尺寸，Z-Buffer總是消耗最長時間，徑向掃描方法次之，螺旋掃描方法又次之，最優(yōu)分割徑向掃描方法消耗的時間最短.這說明最優(yōu)分割徑向掃描方法在這些方法中具有最高效率.考慮到其他基于Z-Buffer原理的方法或改進算法往往比基于柵格數(shù)據的Z-Buffer消耗更長的時間［11］，可認為最優(yōu)分割徑向掃描方法是具有最高效率的遮蔽檢測方法.此外，值得注意的是，最優(yōu)分割徑向掃描方法的效率曲線表現(xiàn)為線性函數(shù)，而ZBuffer方法的效率曲線表現(xiàn)為高次函數(shù).隨著給定DSM尺寸逐漸增大，兩種方法消耗的時間差距會越來越大.

3 結語

近十幾年來，國內外學者努力嘗試尋找一種高效、穩(wěn)健、精確的遮蔽檢測方法.已有的基于ZBuffer技術的遮蔽檢測方法并不適應攝影測量的特殊要求，已經被證明存在諸多難以克服的問題，而一些改進算法往往過于復雜，計算不穩(wěn)定且效率低下.實驗表明，基于徑向投影高度角變化的方法可以有效解決因地形起伏、分辨率差異和高窄墻體引起的各類偽可見和偽遮擋問題，實現(xiàn)穩(wěn)健、精確的遮蔽檢測.對比實驗表明，同樣數(shù)據量下該方法消耗時間最少，且數(shù)據量增大時的時間增量也最少.本文的理論推導和實驗表明，最優(yōu)分割的徑向掃描方法可以實現(xiàn)穩(wěn)健、精確、快速、高效的遮蔽檢測，有利于大規(guī)模高質量真正射影像制作.

［1］ Amhar F，Jansa J，Ries C.The generation of true orthophoto using a3D building model in conjunction with a conventional DTM ［J］.International Journal of Photogrammetry and Remote Sensing，1998，32（4）：16.

［2］ Skarlatos D.Orthophoto production in urban areas［J］.Photogrammetric Record，1999，16（94）：643.

［3］ Rau J Y，Chen N，Chen L.True orthophoto generation of built-up areas using multi-view photos［J］.Photogrammetric Engineering and Remote Sensing，2002，68（6）：581.

［4］ Bang K，Habib A F，Kim K，et al.Comprehensive analysis of alternative methodologies for true ortho-photo generation from high resolution satellite and aerial imagery［C］／／The 5th International Symposium on Mobile Mapping Technology MMT’07.Padua：International Society for Photogrammetry and Remote Sensing，2007：1-5.

［5］ 鐘成，黃先鋒，李德仁，等.真正射影像生成的多邊形反演成像遮蔽檢測算法［J］.測繪學報，2010，39（1）：59.

ZHONG Cheng，HUANG Xianfeng，LI Deren，et al.Polygon based inversion imaging for occlusion detection in true orthophoto generation［J］.Acta Geodaetica et Cartographica Sinica，2010，39（1）：59.

［6］ Braun J.Aspects on true orthophoto production［C］／／Proceedings of 49th Photogrammetric Week.Heidelberg：［s.n.］，2003：205-214.

［7］ Sheng Y W.Minimising algorithm-induced artefaccts in orthophoto generation：a direct method implemented in the vector domain［J］.The Photogrammetric Record，2007，22（118）：151.

［8］ Rau J Y，Chen N Y，Chen L C.Hidden compensation and shadow enhancement for true orthophoto generation［C／CD］∥Proceedings the of 21st Asian Conference on Remote Sensing.Taipei：［s.n.］，2000.

［9］ Sheng Y W，Gong P，Biging G.True orthoimage production for forested areas from large scale arial photograph［J］.Photogrammetric Engineering and Remote Sensing，2003，69（3）：259.

［10］ Zhong C，Li H，Li Z，et al.A vector-based backward projection method for robust detection of occlusions when generating true ortho photos［J］.GIScience and Remote Sensing，2010，47（3）：412.

［11］ Habib A F，Kim E，Kim C.New methodologies for true orthophoto generation［J］.Photogrammetric Engineering and Remote Sensing，2007，73（1）：699.