徐辛超，徐愛功，劉少創(chuàng)，徐宗秋

（1.遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院，遼寧阜新 123000；2.中國科學(xué)院遙感應(yīng)用研究所，北京 100101）

一、引 言

我國的探月工程是《國家中長期科學(xué)和技術(shù)發(fā)展規(guī)劃綱要（2006—2020年）》中的16個重大專項之一，也是中華民族歷史上重大的科技工程之一。2013年發(fā)射的嫦娥三號衛(wèi)星將攜帶月面著陸器和月面巡視探測器實現(xiàn)月球軟著陸和巡視探測任務(wù)［1］。由于月面巡視探測器在行進過程中可能會因地形判斷錯誤陷入撞擊坑或發(fā)生側(cè)翻事故，為確保月面巡視探測器安全開展科學(xué)探測任務(wù)，首先需要以著陸區(qū)附近的三維地形為依據(jù)生成移動指令，完成巡視器的路徑規(guī)劃［2］。在著陸區(qū)中，作為月面最典型的地貌特征，撞擊坑地形的出現(xiàn)是不可避免的，且該地形也容易導(dǎo)致側(cè)翻事故的發(fā)生，因此，撞擊坑地形的恢復(fù)具有非常重要的意義。

以攝影測量為基礎(chǔ)的測圖技術(shù)是測繪領(lǐng)域中獲取地形的主要手段，即利用不同位置拍攝的影像，通過內(nèi)外方位元素計算得到匹配點對的空間三維坐標，并最終由這些匹配點對實現(xiàn)地形恢復(fù)［3-4］。但在紋理缺乏、反射率單一的月表，只能得到較少量的匹配點，通過匹配點也只能得到精度較低的地形，不能滿足巡視器路徑規(guī)劃的需求。SFS（shape from shading）是計算機視覺領(lǐng)域中的一種三維恢復(fù)技術(shù)，即通過影像中的各個像素灰度和物體表面反射模型為基礎(chǔ)，建立反射方程，進而求解得到物體的三維形狀［5］。SFS以影像像素為單位進行三維恢復(fù)，但由于SFS算法在恢復(fù)時要求物體表面的反射特性一致，因此其應(yīng)用領(lǐng)域大大受到限制，并且現(xiàn)有關(guān)于SFS方面的算法較少，恢復(fù)精度有限。因此，研究一種適合月面撞擊坑地形恢復(fù)的SFS算法可以為探月工程提供一定的技術(shù)保障。

二、透視投影下的SFS技術(shù)

SFS又被稱為明暗恢復(fù)形狀，主要用于計算機視覺領(lǐng)域。在地面影像中，由于地表物體反射率往往不一致，因此不能運用SFS技術(shù)進行恢復(fù)。但月球表面較為均勻的月壤成分為SFS技術(shù)的應(yīng)用提供了基礎(chǔ)。在傳統(tǒng)SFS技術(shù)中，一般假設(shè)影像為正射投影，但在著陸器下降過程中獲得的降落影像，由于相機焦距和角度不變，當著陸器高度較低時，影像的投影關(guān)系應(yīng)視為透視投影，因此本文需要研究一種透視投影下的SFS技術(shù)。

1.傳統(tǒng)SFS技術(shù)

運用SFS進行三維恢復(fù)時，要求物體表面的光照反射率已知，并且表面各點對應(yīng)的光照反射率幾乎一致，光源為已知位置的無窮遠處的點光源；此外，影像拍攝過程可近似認為是正交投影［6］。在滿足上述要求的基礎(chǔ)上，按照一定的光照反射模型即可進行三維恢復(fù)［7］。

本文月面地形恢復(fù)過程中，采用探月一期處理影像時采用的反射模型——Lommel-Seeliger反射模型建立反射方程并求解［8］。圖1為Lommel-Seeliger反射模型示意圖，其中，e為觀測者所在方向與物體表面法向量的夾角；i為光照入射角；L為光源方向；E為觀察方向；N為物體表面法向量。

圖1 Lommel-Seeliger反射模型

假設(shè)光源方向向量為N（ps，qs，-1），相機的方向向量為N（pe，qe，-1），物體表面點的法向量為N（nx，ny，nz）。為減少反射方程中的未知數(shù)，通常采用梯度形式的法向量N（p，q，-1）代替N（nx，ny，nz）。根據(jù)反射模型建立反射方程

在反射方程中，p、q是未知的，但對于影像中的每個像素，只能建立一個方程。因此，SFS實際是個病態(tài)的問題，需要將其正則化后，通過方程的求解得到最終三維恢復(fù)結(jié)果。

2.透視投影下的SFS

當相機位置與月表距離較小時，影像的投影方式已經(jīng)不能簡單地近似為正射投影［9］。因此，需要研究透視投影情況下的SFS。

如圖2所示，在像平面內(nèi)的兩條直線c1、c2分別平行于v軸和u軸，其在實際曲面S中對應(yīng)的兩條曲線分別為C1、C2。像平面內(nèi)P點對應(yīng)曲面S中的Q點。

圖2 透視投影模型

將曲面方程寫為參數(shù)形式后，可得透視投影下像平面內(nèi)的曲線方程

假設(shè)f為相機焦距，由透視投影與正射投影關(guān)系可得對應(yīng)正射投影下的曲線方程

對應(yīng)曲線C（s）的切線為

由于過點P的直線c1、c2平行于坐標軸且分別交坐標軸于v0和u0處，因此c1、c2可記為c1（u）、c2（v），實際曲面中對應(yīng)的曲線方程記為C1（u）、C2（v），求取切線方程

將二者叉乘即可得到實際曲面上點Q的法向量N，將其改寫為梯度形式，即

由梯度定義可得p=zu/z，q=zv/z，以 Lommel-Seeliger反射模型建立透視投影下的反射方程，并整理得

得到透視投影下的反射方程后，經(jīng)過一定的約束求解過程，即可進行三維恢復(fù)。

三、透視投影下的撞擊坑三維恢復(fù)

由于SFS建立的反射方程實際上是病態(tài)的方程，因此在建立透視投影下的反射方程后，需要增加一定的約束條件才能實現(xiàn)正確的三維恢復(fù)。

1.邊緣線約束

邊緣信息是影像中一個最基本的重要特征，在影像分析中起著重要的作用。邊緣是指影像中局部特征突變的部分，在一定程度上反映出地形的特點。邊緣線可以反映不同地形的變化，如撞擊坑與周圍地形的交界、斜坡面與月溪的交界等。

為方便問題描述，可以將撞擊坑地形看做由垂直于其邊緣線的寬度為ds的微小面塊組成。如圖3所示，建立OXYZ坐標系統(tǒng)，假設(shè)某邊緣點處為原點O，沿著邊緣線方向為Y軸方向，Z軸指向天頂方向，X軸和Y軸、Z軸構(gòu)成右手空間直角坐標系。

圖3 OXYZ坐標系統(tǒng)

在與邊緣線相鄰的微小區(qū)域內(nèi)，地形可近似看做是微小斜面，且該斜面過Y軸。假設(shè)其與XY平面夾角為α，其中0＜α＜90°，則該斜面方程在此坐標系下可以描述為

可得該區(qū)域內(nèi)任意點（x，y，z）處的梯度形式法向量為（sin α/cos α，0，-0/cos α）。該向量在XY平面內(nèi)的投影為（sin θ/cos θ，0），由此可得該區(qū)域內(nèi)任意地形點的梯度比例k在XY平面內(nèi)的投影為

即該微小斜面內(nèi)的各點法向量在XY平面內(nèi)的投影垂直于Y軸。

在OXYZ坐標下，其對應(yīng)邊緣線與Y軸重合，因此邊緣線上的任意點法向量在XY平面內(nèi)的投影也垂直于Y軸，即與邊緣線相鄰的微小斜面內(nèi)各點的法向量在XY平面內(nèi)的投影與邊緣線點法向量在XY平面內(nèi)的投影平行。因此，該微小斜面內(nèi)的梯度比例即可由其附近邊緣線的梯度比例求得。

2.連續(xù)性約束

一般情況下撞擊坑內(nèi)部地形表面都是連續(xù)的。為求取整個撞擊坑內(nèi)部其他點的梯度比例，本文算法在此運用了Horn提出的連續(xù)性約束條件。假設(shè)Ω為影像（x，y）的覆蓋范圍，連續(xù)地形表面模型可以表示為

式中，px、py和qx、qy分別為p和q關(guān)于x、y方向的偏導(dǎo)數(shù)。

式（10）也可稱為光滑約束條件，該條件下地形表面相鄰各點的表面法向量是相似的。因此，相鄰各點的梯度比例也是相似的。撞擊坑內(nèi)部除邊緣線附近點外的其他點的梯度比例可按照連續(xù)性約束條件演化得到。

至此，根據(jù)本文提出的邊緣線約束和連續(xù)性約束條件，即可運用撞擊坑邊緣線點法向量的梯度比例演化求得撞擊坑內(nèi)部任意點處的梯度比例k。將本文算法提出的梯度比例因子代入式（7），得到加入本文約束比例k后的反射方程為

顯然通過梯度比例與地形表面連續(xù)性約束后，本文算法在透視投影下的反射方程是正則化的，通過一定方法進行求解后，即可實現(xiàn)三維恢復(fù)。

四、試驗分析

為了對本文提出算法的可行性進行論證，采用模擬影像進行了算法正確性測試，然后對月面影像進行了恢復(fù)測試，并與經(jīng)典SFS方法進行了比較分析。

1.算法實現(xiàn)

由第三節(jié)得到的撞擊坑的邊緣線點處的梯度比例可求得整個撞擊坑內(nèi)部點的梯度比例，從而實現(xiàn)反射方程的約束，使得SFS問題能夠正確求解。撞擊坑地形恢復(fù)流程如圖4所示。

圖4 透視投影下的撞擊坑地形恢復(fù)流程

首先通過灰度卷積模板，得到撞擊坑初步的邊緣位置；然后對每個邊緣像素周圍的邊緣點進行擬合，得到該邊緣點處的梯度比例；此后按照曲面連續(xù)性約束條件，得到整個撞擊坑的梯度比例；最后將得到的梯度比例代入透視投影下的反射方程，求解得到撞擊坑的三維地形。

2.影像三維恢復(fù)

本文首先采用模擬影像對算法的可行性進行驗證。圖5為本文提出算法對模擬影像的恢復(fù)結(jié)果及恢復(fù)誤差。模擬影像對應(yīng)光照方向向量為（0.05，0.1，1）。

圖5 模擬影像三維恢復(fù)結(jié)果及恢復(fù)誤差

將圖5中本文算法恢復(fù)結(jié)果與模擬影像真實值進行對比可以得出，本文算法的恢復(fù)結(jié)果與模擬影像結(jié)果非常接近，證明了本文算法的正確性和可行性，但仍存在一定程度的恢復(fù)誤差。本文算法的平均誤差為0.065，最大誤差為0.157（由于SFS恢復(fù)結(jié)果為相對關(guān)系，在此只采用相對數(shù)值，無具體單位）。撞擊坑內(nèi)部點恢復(fù)誤差較小，而誤差較大的部分集中在撞擊坑邊緣附近。本文算法是基于地形特征邊緣線的梯度約束展開，因此恢復(fù)誤差的產(chǎn)生必然與特征邊緣線有關(guān)。由于算法的局限性，根據(jù)灰度算子提取得到的特征邊緣線與影像的實際邊緣不可避免地存在位置偏差；此外，在梯度比例因子求取過程中還需要將特征邊緣線進行擬合，在此過程中會產(chǎn)生擬合誤差。因此，在真實邊緣點與所提取得到的邊緣點附近必然會存在較多的恢復(fù)誤差。

本文采用真實月面影像對本文提出的算法及SFS算法中對真實野外影像恢復(fù)效果最好的Tsai算法進行了恢復(fù)效果測試，并得出了測試結(jié)果。真實月面影像中光照方向未知，因此在進行地形恢復(fù)時需要首先運用現(xiàn)有光照估計方法對光照方向進行估計。本文采用文獻［11］中的方法進行估計，得到圖6中月面撞擊坑影像的光照方向分別為（1.0，1.2，1.5）和（-0.1，-1.2，1.5）。Tsai算法與本文算法的最終恢復(fù)效果如圖6所示。

圖6 月面影像三維恢復(fù)

由圖6可以看出，這兩種算法對于陰影區(qū)域的恢復(fù)效果都比較差。由于SFS恢復(fù)過程中，需要首先根據(jù)影像中的像素灰度分別建立反射方程，而陰影區(qū)域的像素灰度實際上不符合月面反射模型特征，在此基礎(chǔ)上建立的反射方程也是錯誤的，從而導(dǎo)致最終的三維恢復(fù)結(jié)果存在較大誤差；此外，在影像中存在某些點的反射率突變，直接造成恢復(fù)結(jié)果的突變。因此，為保證恢復(fù)結(jié)果的精度，運用SFS技術(shù)進行三維恢復(fù)時，要求物體表面的反射率近似一致。雖然Tsai算法在現(xiàn)有SFS算法中對于實際地形的恢復(fù)效果最接近實際地形，但從測繪角度而言，該方法對于實際地形的恢復(fù)效果仍不能滿足生產(chǎn)需要，而本文算法由于存在邊緣梯度約束與撞擊坑內(nèi)部表面連續(xù)性約束，因此得到的恢復(fù)結(jié)果與實際地形更加接近。但由于透視投影下本文采用的反射模型與實際反射情況的偏差、光照方向估計偏差，以及SFS約束條件的求取等因素導(dǎo)致本文算法仍存在一定的恢復(fù)誤差。

五、結(jié)束語

本文研究了月面影像的特點，以及通過匹配進行月面地形恢復(fù)的局限性，分析了通過計算機視覺中SFS技術(shù)進行月面地形恢復(fù)的可行性。在對透視投影情況下SFS算法進行深入研究的基礎(chǔ)上，建立了透視投影下的反射方程，以Lommel-Seeliger反射模型為基礎(chǔ)，通過月面影像中存在的特征信息，實現(xiàn)了對反射方程的正則化約束，最終實現(xiàn)了對月面地形的三維恢復(fù)。采用Matlab模擬影像和真實月面影像對本文算法的正確性進行了測試。試驗證明，本文算法在模擬影像情況下的恢復(fù)結(jié)果與真實值非常接近，驗證了算法的正確性。此外，相對于其他SFS算法，本文算法對于實際月面影像有更好的恢復(fù)效果。SFS技術(shù)用于地形恢復(fù)不僅可以豐富測繪手段，也可為我國深空探測領(lǐng)域中測繪方面的問題提供一定的技術(shù)儲備。

［1］ 歐陽自遠，李春來，鄒永廖，等.中國月球探測二期工程科學(xué)目標的研究［C］∥中國宇航學(xué)會深空探測技術(shù)專業(yè)委員會第一屆學(xué)術(shù)會議論文集.北京:［s.n.］，2005:1-8.

［2］ 葉培建，張火高，饒煒.積極應(yīng)對深空探測的技術(shù)挑戰(zhàn)［C］∥中國宇航學(xué)會深空探測技術(shù)專業(yè)委員會第二屆學(xué)術(shù)會議論文集.北京:［s.n.］，2005:1-8.

［3］ 孫敏，胡爭.基于向量匹配的稀疏深度圖生成算法［J］.測繪學(xué)報，2011，40（1）:90-95.

［4］ 李壯，雷志輝，于起峰.基于梯度徑向夾角直方圖的異源圖像匹配［J］.測繪學(xué)報，2011，40（3）:318-325.

［5］ PRADHAN R，GHOSE M K，JEYARAM A.Extraction of Depth Elevation Model（DEM）from High Resolution Satellite Imagery Using Shape from Shading Approach［J］.International Journal of Computer Applications，2010，7（12）:40-46.

［6］ BURT P J，ADELSON E H.High-frequency Shape and Albedo from Shading Using Natural Image Statistics［C］∥IEEE Conference on Computer Vision and Pattern Recognition.Los Angeles:［s.n.］，2011:2521-2528.

［7］ KHUU S K，KHAMBIYE S.The Influence of Shape from Shading Information on the Perception of Global Motion［J］.Vision Research，2012，55（15）:1-10.

［8］ 歐陽自遠，胡浩，王世杰，等.月球科學(xué)概論［M］.北京:中國宇航出版社，2005:38-40.

［9］ TANKUS A，SOCHEN N，YESHURU Y.Shape from Shading under Perspective Projection［J］.International Journal of Computer Vision，2005，63（1）:21-43.

［10］ TSAI P S，SHAH M.Shape from Shading Using Linear Approximation［J］.Image and Vision Computing，1994，12（8）:487-498.

［11］ 楊杰，鄧志鵬，郭英凱，等.兩種新的光照方向估計方法［J］.上海交通大學(xué)學(xué)報，2002，36（6）:894-896.