王亞林，劉 鵬，吳輝陽，黎 衡，趙 巍

（1.哈爾濱工業(yè)大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，哈爾濱 150001；2.中國人民解放軍63861部隊，白城 137001）

引 言

小行星是太陽系的重要組成部分，對地球環(huán)境的演化和人類的生存都有重大影響。一方面，小行星保存了較多太陽系早期的形態(tài)、特性，可以用來研究太陽系的形成和演變過程；另一方面，經(jīng)不懈努力所獲得的小天體知識將有助于人類進(jìn)一步對太陽系深度探索，如天體表面取樣用于評估小天體蘊含礦藏的可利用性等。

在過去的20年里，隨著航天技術(shù)的發(fā)展，小行星探測任務(wù)從地基、空基望遠(yuǎn)鏡觀測逐漸延伸到航天器的近距離觀測，比如2000年，NASA 的近地小行星交會（Near Earth Asteroid RendezvoHs，NEAR）探測器成功環(huán)繞近地小行星Eros（433），并最終降落在Eros 的表面；2005年日本“隼鳥號”（Hayabusa）成功在近地小行星Itokawa（25143）上進(jìn)行了一觸即走式著陸，采集了小行星表面物質(zhì)；2012年12月13日，我國的“嫦娥2 號”飛船近距離飛越小行星Toutatis（4179）并獲取了清晰的遙感圖像。

目前多個國家制定了明確的小行星探測規(guī)劃，如日本正在進(jìn)行的“隼鳥2號”任務(wù)（2014年12月3日發(fā)射），目標(biāo)為近地小行星1999JU3（162173）；美國正在開展的 OSIRIS-REx 任務(wù) （2016年 9月 8日發(fā)射），目標(biāo)為小行星1999RQ36（101955）；歐洲航天局（European Space Agency，ESA）MarcoPolo-R 任務(wù)，目標(biāo)為近地小行星2008EV5。

深空探測任務(wù)既需要充分掌握某潛在目標(biāo)的軌道、動力學(xué)、引力場等特性，也需要對目標(biāo)小天體的著陸取樣返回系列任務(wù)提供可靠的技術(shù)保障。小天體表面安全著陸是對探測器的著陸導(dǎo)航機(jī)構(gòu)和著陸導(dǎo)航策略的極大考驗。2005年11月25日日本的“隼鳥號”子飛船密涅瓦在登陸過程中迷失失聯(lián)；2014年11月22日ESA“羅塞塔號”的“菲萊”著陸器在登陸彗星67P/Churyumov-Gerasimenko 時反彈并進(jìn)入了一個光線昏暗的巨石堆/碎石堆區(qū)域。

對于尺寸很小（約1 km）的碎石堆小行星而言，其引力場低、逃逸速度低（約10-2~10-1m/s），對此類小天體開展著陸段導(dǎo)航和表面操作（如錨定、取樣、返回等）的實時性、安全性和精度有高要求。視覺輔助的地形相對導(dǎo)航可以提高探測器著陸的成功率。因此，需要構(gòu)建一個高精度的小行星地形模型用于著陸機(jī)構(gòu)安全性可靠性驗證和地形相對導(dǎo)航的算法驗證。

目前研究人員對行星、小天體表面地形的仿真驗證平臺的構(gòu)建方法主要有：Parkes 等人使用PANGU對星體表面進(jìn)行模擬，著重做較大行星表面仿真[1]；Sánchez等使用軟球離散單元法（Soft-Sphere Discrete Element Method，SSDEM）來模擬碎石堆小行星的旋轉(zhuǎn)重塑，重點分析組成小行星的碎石之間的作用力[2]；Van Wal 提出了一種高保真的小天體表面著陸仿真方法，重點在于著陸器與地形之間的碰撞分析[3]。這些方法多是建立較為粗糙的整體地形模型，沒有構(gòu)建出局部高精度地形。

本文根據(jù)小天體表面操作執(zhí)行機(jī)構(gòu)和著陸段導(dǎo)航方法驗證需求，選擇并仿真生成5 種典型局部地形。在不考慮小行星表面引力、熱慣量、自轉(zhuǎn)周期等物理特性的條件下，通過計算機(jī)仿真的手段評估仿真參數(shù)對地形生成結(jié)果的影響以及與真實地形間的相似程度，為機(jī)構(gòu)及方法驗證提供數(shù)據(jù)支持。

1 小行星地形研究

1.1 地形特性

研究表明，小行星是由重力和凝聚力結(jié)合在一起的碎石堆[4-5]，其表面覆蓋有風(fēng)化層。風(fēng)化層主要由厘米級至數(shù)米級大小不一的石塊構(gòu)成。Miyamoto等[6]通過研究小行星Itokawa 高分辨率的圖像發(fā)現(xiàn)這些石塊在干燥、真空和微重力的環(huán)境下，由于振動引發(fā)了整個小行星的顆?；^程；由于Itokawa的體積較小，這種顆粒化的過程成為了其主要的表面重構(gòu)過程，這也意味著該過程很有可能會發(fā)生在其它小行星上。小行星的體積越小，極低的重力允許相對較大的細(xì)顆粒得以逃脫，這些細(xì)顆粒可能是由微型隕石或者電離引起的，當(dāng)這些細(xì)顆粒被運輸漂散，地表就會腐蝕或“收縮”，由此產(chǎn)生了所謂的“巴西堅果效應(yīng)”[7]，即體積比較大的顆粒會上升到表層，而較小的顆粒會沉降到底部；并最終在小行星表面形成了，如圖1（a）所示的地形[8]。

圖1 小行星表面石塊圖片示例Fig.1 Image illustration of real boulders，cobbles and dust on asteroid surface

較小的碎石通常具有更高的流動性，因為它們有較低的摩擦，移動它們所需的振動的振幅也較小。碎石移動性的不同導(dǎo)致了它們的分離，進(jìn)而導(dǎo)致了灰塵池塘（光滑地形）的形成，如圖1（b）所示，這也是光滑地形和粗糙地形之間的邊界如此明顯的原因。

1.2 石塊特性

對碎石堆小行星表面石塊的進(jìn)一步研究表明：小行星表面的石塊尺寸概率密度函數(shù)PDFX(d)服從以下冪律分布

其中：α>0被稱作冪指數(shù)；dmin是石塊的最小尺寸。

顯然，對于不同地形區(qū)域其冪指數(shù)的值是不同的。通過對小行星Itokawa 表面石塊數(shù)據(jù)的研究Michikami 等得出結(jié)論[9]，Itokawa 表面的粗糙區(qū)域冪指數(shù)為2.8±0.1，光滑區(qū)域冪指數(shù)為3.2±0.1。Mazrouei指出[10]，在Itokawa表面直徑小于15 m的石塊的冪指數(shù)是2.1±0.1，進(jìn)而在最新的成果[11]中給出修正的冪指數(shù)為3.1～3.5，與其他學(xué)者的成果給出的2～4的范圍相當(dāng)。

另一方面，選擇石塊的最小尺寸dmin可以避免出現(xiàn)當(dāng)dmin趨近于0 時，PDF趨近于1 的畸形情況；同時，為了避免生成超大石塊，定義石塊的最大尺寸dmax。假設(shè)y是[0,1]之間的均勻分布，那么可以由以下公式得出服從冪律分布的d

小行星表面的石塊數(shù)量服從下述公式

其中：k0是單位面積上半徑大于d0的石塊數(shù)量；k是單位面積上半徑大于dmin的石塊數(shù)量。根據(jù)Mazrouei的研究結(jié)果[10]，小行星Itokawa 單位面積上半徑大于d0=2.5 m的石塊數(shù)量為k0=2.05×10-3m-2。

在對Itokawa 的石塊外形分析中[9]，指出石塊的長寬高比例大致為或簡記為a∶b∶c=1∶0.7∶0.5。

2 小行星地形模擬生成方法

2.1 地形模型

在給定了數(shù)目且服從冪律分布的石塊以后，還不能確定這些石塊在小行星表面的分布，原因是由公式（3）計算出的單位面積上大于某一尺寸的石塊數(shù)量是基于整個小行星的平均數(shù)據(jù)。但是大致的分布趨勢可以確定如下[12]：在重力勢較高的地方，均為厘米級大小的碎石鋪滿的“池塘”；而在重力勢較低的地方，大多是較大塊的石塊。為了進(jìn)一步區(qū)分不同引力場位置的局部地形，本文按照地形區(qū)域內(nèi)石塊數(shù)量及分布情況將地形分為以下5類：

1）灰塵池塘：只含有碎石；

2）石塊密布區(qū)：含有任意大小的石塊；

3）石塊區(qū)和灰塵池塘的過渡區(qū)域；

4）灰塵池塘中散布石塊；

5）石塊密布區(qū)中有小塊灰塵池塘。

圖2 5種典型地形的真實示例（圖片來源于Itokawa小行星，編號依次為st_2563537820、st_2530297837、st_2539423137、st_2563511720、st_2532629277）Fig.2 Image illustration of five typical real terrains from Itokawa asteroid with image label st_2563537820，st_2530297837，st_2539423137，st_2563511720，and st_2532629277

2.2 石塊模型

根據(jù)石塊的尺寸大小，將其分為3 種：碎石（dust）半徑0.02～0.05 m；鵝卵石（cobble）半徑0.05～0.5 m；巨石（boulder）半徑0.5 m以上。

由于石塊確切形狀既不是立方體也不是球體，或者其它立方體，因此石塊的模型必須有足夠的復(fù)雜性，同時必須有足夠的簡單性以供大范圍的模擬，選擇隨機(jī)形變的二十面體來表征石塊。石塊按照以下步驟生成：

1）生成一個直徑為d的正二十面體Stone=[Point,Face]；

2）對二十面體的12 個頂點進(jìn)行最大偏差為d/4的隨機(jī)伸縮Point=Point+StretchMatrix

3）對二十面體進(jìn)行a∶b∶c=1∶0.7∶0.5 的比例縮放Point=Point?ReshapeMatrix

4）在3 個方向上隨機(jī)選取旋轉(zhuǎn)角度，對石塊進(jìn)行旋Point=Point?RotationMatrix

式中α、β、γ是繞X軸、Y軸、Z軸旋轉(zhuǎn)的角度。

圖3依次展示了原始、隨機(jī)形變、比例縮放、隨機(jī)旋轉(zhuǎn)的石塊。

圖3 原始、隨機(jī)形變、比例縮放、隨機(jī)旋轉(zhuǎn)的石塊模型Fig.3 Illustration of regular icosahedron and its random deformation，non-uniform scaling，and random rotating in simulating boulders

2.3 地形生成方法

按照仿真需求，需要生成一個3 m×3 m 大小的仿真著陸區(qū)，每像素分辨率為區(qū)域邊長的千分之一。經(jīng)計算，一個1 024×1 024像素的圖片，每像素分辨率為3 mm 即可滿足仿真區(qū)域的需求。使用兩個1 024×1 024大小的圖分別存儲每個像素點的高度信息和石塊信息。

2.3.1 石塊放置方法

仿真地形按照以下步驟生成：

1）在1 024×1 024 的平面上，隨機(jī)選取位置放置石塊并以隨機(jī)深度下沉，更新高度和石塊信息（石塊編號），直至平面的鋪滿度達(dá)到某一閾值（例如95%）或者達(dá)到石塊數(shù)目上限；

2）在隨機(jī)選取位置的過程中，若該位置的石塊信息已被更改（初始為0），則再次隨機(jī)選取位置。此做法能保證石塊總是鋪在未曾放置石塊的位置，直至鋪滿度達(dá)到閾值；

3）石塊的詳細(xì)信息按照編號存儲在列表中。

地形I、II、IV 中石塊位置的隨機(jī)選取都服從均勻分布，為了實現(xiàn)地形III中石塊到灰塵池塘的過渡，大石塊的位置分布不再服從均勻分布，而是服從p(x,y)=a(x+y)-k的概率分布，其中a是歸一化因子，k=0.5，即：距離左上角越近，大石塊的分布概率越高；而在地形V 中設(shè)定了一個均值為d/3、方差為d/6 的圓形的區(qū)域，當(dāng)大石塊落入這一區(qū)域時，再次隨機(jī)生成一個位置。

石塊放置時，將石塊坐標(biāo)系的原點放置到選定的位置上，并將其隨機(jī)下沉，最后再將石塊的各點坐標(biāo)由石塊本體坐標(biāo)系o-xyz轉(zhuǎn)到為地形坐標(biāo)系O-XYZ下，如圖4所示。

圖4 地形坐標(biāo)系與石塊本體坐標(biāo)系Fig.4 Relation between terrain coordinate system and boulder’s body-fixed coordinate system

2.3.2 石塊高度更新方法

為了更新地面高程信息，需要得到石塊每個面上所有的XY坐標(biāo)值為整數(shù)值（代表一個像素點位置）的點(x,y,z)，稱為像素點。三維空間的任意平面可以由Kx+My+Lz+N=0 表示，對于一個石塊的某一個確定的面，只要確定x,y的值就可以根據(jù)平面公式計算出z值，則只需尋找所需的x,y值，尋找一個石塊一個面上所有的像素點通過以下步驟完成：

1）將石塊的面ABC沿A'B'軸方向投影到XOY平面上得到三角形A'B'C'，三角形A'B'C'邊上和內(nèi)部的所有XY坐標(biāo)為整值的點即為所求點，如圖5（a）；

2）計算三角形A'B'C'的外切矩形，將其中所有XY坐標(biāo)為整數(shù)值的點作為候選點；

3）要得到三角形A'B'C'所有邊上的點，因為采用像素的形式，XY坐標(biāo)只能是離散整數(shù)值，例如構(gòu)成邊A'B'的像素點不嚴(yán)格滿足其直線方程K'x+M'y+L'=0,故此處采用計算機(jī)圖形學(xué)中的直線掃描轉(zhuǎn)換方法之一：Bresenham 畫線算法（Bresenham's line algorithm）來確定構(gòu)成該條線段的像素點，如圖5（b）；

4）要得到三角形內(nèi)部的點，使用叉乘法判斷：平面內(nèi)一點P在三角形A'B'C'內(nèi)部當(dāng)且僅當(dāng)三組向量方向相同。對所有候選點進(jìn)行如上檢測，如圖5（c）；

5）根據(jù)所有已經(jīng)獲得的(x,y)，計算出其對應(yīng)在平面ABC的坐標(biāo)(x,y,z)，返回所有的(x,y,z)。

2.4 算法流程

總結(jié)2.1~2.3 中的小行星地形模擬生成方法，得到算法的流程圖，如圖6所示。

圖5 地形高度更新方法示意圖Fig.5 Main steps of terrain altitude updating

圖6 算法流程圖Fig.6 Flowchart of asteroids topography generating algorithm

3 實驗結(jié)果

在windows10 平臺下按照第3 部分的仿真方法，利用matlab實現(xiàn)了仿真算法，并對仿真關(guān)鍵參數(shù)：冪律指數(shù)α進(jìn)行了研究。由第2 部分小行星地形知識，學(xué)者得出的冪律指數(shù)范圍為2~4，但是結(jié)果來源于直徑5 m 以上[12]或直徑15 m 以下的石塊[11]的統(tǒng)計結(jié)果，本文的地形仿真中石塊尺寸在直徑1 m以下，故擴(kuò)大了冪律指數(shù)的取值區(qū)間，對冪律指數(shù)α按2~5的范圍進(jìn)行仿真驗證。其中圖7(a)α為 4.8，7(d)α為 4.2，7(b)、7(c)及7(e)的α均為2.8。

由圖7(a)(d)(b)可知隨著冪律指數(shù)α變小，仿真地形上大尺寸的石塊隨之增多，分別對應(yīng)著地形I、II、IV。圖7(b)(c)(e)中的冪律指數(shù)α均為2.8，大尺寸石塊數(shù)量大體相同，但由于其不同的放置策略得出了地形II、III、V。

由圖8可知，石塊尺寸分布直方圖大致相同，只是其斜率有所區(qū)別，由于坐標(biāo)系為對數(shù)坐標(biāo)，這表明生成的石塊尺寸服從冪律分布。

圖7 五種地形示例、仿真結(jié)果、地形仿真參數(shù)及石塊尺寸-頻率直方圖（對數(shù)坐標(biāo)）Fig.7 Five typical terrains，simulation results，terrain parameters and boulder-size-frequency histogram in log coordinates

4 結(jié) 論

本文通過對碎石堆小行星地表及石塊特性的研究，建立了地形仿真模型及方法。仿真結(jié)果表明該方法生成的仿真地形與實際探測得到的小行星表面局部地形相比一致性較好，可以用于著陸器的著陸機(jī)構(gòu)驗證和著陸器著陸導(dǎo)航及視景仿真。