碎石堆構(gòu)造小行星表面地形分析與仿真驗(yàn)證

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

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

引 言

小行星是太陽系的重要組成部分，對地球環(huán)境的演化和人類的生存都有重大影響。一方面，小行星保存了較多太陽系早期的形態(tài)、特性，可以用來研究太陽系的形成和演變過程；另一方面，經(jīng)不懈努力所獲得的小天體知識將有助于人類進(jìn)一步對太陽系深度探索，如天體表面取樣用于評估小天體蘊(yù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)航策略的極大考驗(yàn)。2005年11月25日日本的“隼鳥號”子飛船密涅瓦在登陸過程中迷失失聯(lián)；2014年11月22日ESA“羅塞塔號”的“菲萊”著陸器在登陸彗星67P/Churyumov-Gerasimenko 時反彈并進(jìn)入了一個光線昏暗的巨石堆/碎石堆區(qū)域。

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

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

本文根據(jù)小天體表面操作執(zhí)行機(jī)構(gòu)和著陸段導(dǎo)航方法驗(yàn)證需求，選擇并仿真生成5 種典型局部地形。在不考慮小行星表面引力、熱慣量、自轉(zhuǎn)周期等物理特性的條件下，通過計(jì)算機(jī)仿真的手段評估仿真參數(shù)對地形生成結(jié)果的影響以及與真實(shí)地形間的相似程度，為機(jī)構(gòu)及方法驗(yàn)證提供數(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ì)顆?？赡苁怯晌⑿碗E石或者電離引起的，當(dāng)這些細(xì)顆粒被運(yùn)輸漂散，地表就會腐蝕或“收縮”，由此產(chǎn)生了所謂的“巴西堅(jiān)果效應(yīng)”[7]，即體積比較大的顆粒會上升到表層，而較小的顆粒會沉降到底部；并最終在小行星表面形成了，如圖1（a）所示的地形[8]。

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

較小的碎石通常具有更高的流動性，因?yàn)樗鼈冇休^低的摩擦，移動它們所需的振動的振幅也較小。碎石移動性的不同導(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）計(jì)算出的單位面積上大于某一尺寸的石塊數(shù)量是基于整個小行星的平均數(shù)據(jù)。但是大致的分布趨勢可以確定如下[12]：在重力勢較高的地方，均為厘米級大小的碎石鋪滿的“池塘”；而在重力勢較低的地方，大多是較大塊的石塊。為了進(jìn)一步區(qū)分不同引力場位置的局部地形，本文按照地形區(qū)域內(nèi)石塊數(shù)量及分布情況將地形分為以下5類：

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

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

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

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

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

圖2 5種典型地形的真實(shí)示例（圖片來源于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 個頂點(diǎn)進(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)計(jì)算，一個1 024×1 024像素的圖片，每像素分辨率為3 mm 即可滿足仿真區(qū)域的需求。使用兩個1 024×1 024大小的圖分別存儲每個像素點(diǎn)的高度信息和石塊信息。

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ī)選取都服從均勻分布，為了實(shí)現(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)系的原點(diǎn)放置到選定的位置上，并將其隨機(jī)下沉，最后再將石塊的各點(diǎn)坐標(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ù)值（代表一個像素點(diǎn)位置）的點(diǎn)(x,y,z)，稱為像素點(diǎn)。三維空間的任意平面可以由Kx+My+Lz+N=0 表示，對于一個石塊的某一個確定的面，只要確定x,y的值就可以根據(jù)平面公式計(jì)算出z值，則只需尋找所需的x,y值，尋找一個石塊一個面上所有的像素點(diǎn)通過以下步驟完成：

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

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

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

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

5）根據(jù)所有已經(jīng)獲得的(x,y)，計(jì)算出其對應(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 實(shí)驗(yàn)結(jié)果

在windows10 平臺下按照第3 部分的仿真方法，利用matlab實(shí)現(xiàn)了仿真算法，并對仿真關(guān)鍵參數(shù)：冪律指數(shù)α進(jìn)行了研究。由第2 部分小行星地形知識，學(xué)者得出的冪律指數(shù)范圍為2~4，但是結(jié)果來源于直徑5 m 以上[12]或直徑15 m 以下的石塊[11]的統(tǒng)計(jì)結(jié)果，本文的地形仿真中石塊尺寸在直徑1 m以下，故擴(kuò)大了冪律指數(shù)的取值區(qū)間，對冪律指數(shù)α按2~5的范圍進(jìn)行仿真驗(yà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é)果表明該方法生成的仿真地形與實(shí)際探測得到的小行星表面局部地形相比一致性較好，可以用于著陸器的著陸機(jī)構(gòu)驗(yàn)證和著陸器著陸導(dǎo)航及視景仿真。