0 引言
两个世纪前，人类发现了第一颗小行星[3]。时至今日，已经有超过十万颗小行星被人类所记载。小行星是我们了解宇宙的科学窗口，它们或它们的碎片自太阳系诞生以来已经存在了近46亿年[4]。作为太阳系历史的重要见证，其物质组成与内部结构为科学家深入剖析太阳系的诞生与演变过程提供了珍贵的数据。分析小行星的数据，可以洞察太阳系的形成时间与形成机制，同时也能阐明小行星构建的过程和体系的演变规则。小行星所含的化学构成与有机生物或许对科学家研究地球上生命的来源和演变过程有影响。科学家借助样本分析可以追寻地外生命的证据，揭示多元生命形态在宇宙中的可能性[5-6]。研究发现，小行星蕴含丰富的矿石资源和稀有金属资源，例如铱、铂、镍等含量是地球的数倍[7]，将其获取可以作为未来工业需求的大量原料。同时，小行星数量众多，有着非常高的科研价值与经济效益。随着航天技术的发展，开采小行星上的资源可以缓解地球资源面临枯竭的局面[8]。
小行星采矿是开发太空资源的第一步，也是走向太空、开发太阳系资源的基础[11]。航天大国虽然已经初步完成了小行星样本采集，拥有了近地小行星的绕飞、着陆的能力和经验[12-14]，具备了小行星资源开发的初步技术，但对小行星采矿仍停留在概念阶段。王巍和姚伟[9,11]从宏观层面对小行星采矿进行规划与评估，提出了近地小天体采矿体系和方案框架设想，太阳系资源开发体系架构和“天工开物”倡议设想，并初步建立了可开采矿物总量的评估方法和以资源净收益为目标的小天体采矿效益评估方法。Andreas 等[15]对小行星采矿进行了具体技术分析，为未来的技术发展和性能改进提供了基本建议。Zhang等[16]通过理论分析和数值仿真研究了使用柔性网航天器捕获小行星的过程，为小行星采矿提供了理论基础。本文首先构建了小行星的引力场和表面参数模型，其次提出了一种高效的表面开采方案，选择小行星101955 Bennu为目标小行星，分析其表面动力学特征；最后，构建了柔性绳网的动力学模型，采用集中质量非线性弹簧阻尼模型仿真模拟了柔性绳网捕获小行星的小行星采矿动力学过程，并对仿真结果进行分析。
1 小行星的引力场建模和表面参数化建模
多面体法是当前小行星引力场高精度建模的常用方法之一，虽然其计算效率不如球谐函数法等方法，但由于它的精度优势，适用于小行星任务的高精度离线仿真分析。传统多面体法在计算引力势和引力时，存在多面体表面上的奇异点，当质点运动到多面体侧面的边缘处时，由于多面体表面并不平滑，质点运动会发生速度违约现象。因此，必须将多面体相邻两个分块进行平滑化处理，以保证质点在穿越棱边时速度变化的连续性。为了解决多面体表面奇异问题，本节将多面体模型引力势和引力做了修正，如式（1）和式（2）所示。
式中，$ {\mathrm{G}} $ 为万有引力常数；$\sigma $ 为多面体密度；$FS$ 为多面体的所有侧面构成的集合；$ES$ 为多面体所有棱边构成的集合；变量${\boldsymbol{r}_e}$ 、${\boldsymbol{r}_f}$ 分别为从测试点指向棱$e$ 和面$f$ 任意一点的矢量。${P_e}\left( r \right)$ 和${Q_f}\left( r \right)$ 为分段函数，如式（3）和（4）所示。
球面谐波级数曲面模型是一种经典的参数化建模方法，在球面几何、天文学和计算机图形学等领域被广泛应用。球面谐波级数的基本思想是将球面上的函数表示为一组正交的基函数的线性组合，这些基函数称为球谐函数。考虑到球面谐波可以精确拟合形状复杂的表面，故本文探究此模型对小行星三维模型进行重构。
式中，$ \left( {\theta {\text{,}}\varphi } \right) $ 为小行星表面上任意位置的球坐标；$ R\left( {\theta {\text{,}}\varphi } \right) $ 为原点到该处的距离；$ c_l^m $ 为球谐复系数，在本文中由多面体坐标数据经最小二乘法[17]计算得到；$ Y_l^m\left( {\theta {\text{,}}\varphi } \right) $ 为球谐级数基函数，其计算方法如式（16）所示。
式中，$ P_l^m\left( \cdot \right) $ 表示$ l $ 阶Legendre多项式，如式（7）所示。
在实际应用中，若给定任意球坐标$ \left( {\theta {\text{,}}\varphi } \right) $ ，则可由式（6）给出其表面对应的径向距离。
2 表面开采
2.1 采矿方案设计
2.1.1 采矿阶段设计
本节针对小行星表面开采提出了一种步骤完备、矿物采集效率高、采集矿物种类多样的科学性的新方案。下面将具体开采流程（见图1 ）。
（1）了解目标小行星的矿物质成分和开采前景。确定小行星矿物成分的一个有效的方法是将它们与地球上的陨石进行比较。结合现代遥感技术，我们几乎可以确定一颗小行星含有的矿物质成分，主要方法[19]包括分光光度法、辐射测量法、高光谱影像、热建模等。确定小行星的矿物质成分后发射探测器，至近目标小行星环绕轨道，并向不同方向发射多个携带高清摄像机的球体巡视器探测小行星表面风化程度及矿产种类分布及地形分布。
（2）寻找采矿位置，选择采矿地点。根据巡视器的不同着陆情况，选择不容易发生相对滑动（即临界摩擦系数小），不存在大型岩石且碎石体积较小、分布均匀的位置作为采矿标记点，令探测器分别在这些标记点处作业。
（3）采集表面矿产资源。从标记点中选定一个目标点进行采矿作业，探测器通过内部磁源线圈通电使小行星磁，或通过磁化收集端将表面一些风化碎片（主要包括磁性金属碎片，如铁、镍等，含磁性矿物的岩石碎片，被磁化后的尘埃颗粒等）吸入存储装置中。将风化层用刮板刮下并收集风化碎片，此时刮板安装在探测器表面，不与目标小行星进行接触，方便下一步的操作。
（4）发射爆破装置。发射爆破装置撞击目标小行星表面，并在撞击小行星表面同时产生爆炸（类似隼鸟2号任务中撞击器对小行星162173 Ryugu的撞击）。同时，探测器自身提升轨道高度躲避至安全区域内，防止碎片等喷射物损坏装置，这样就在目标小行星表面形成了一个人造撞击坑。在爆破结束后，探测器重新回到人造撞击坑的上方，再次通过内部磁源线圈通电使小行星磁化，通过磁化收集端将风化碎片吸入存储装置中。
（5）锚定并固定收割航天器。将收割航天器锚定并固定在撞击坑内，并由探测器释放软层薄膜袋，薄膜边缘等间距挂有锚定器，锚定器在遇到小行星表面缝隙时会插入并锚固在其表面。防止挖掘过程中一些碎片颗粒的飘散，进而增大采集效率。
（6）收割航天器切割矿石并收集。对于大块矿石可以用探测器自身携带的切割装置进行切割处理，同时加大磁力从而平衡挖掘时钻井的反作用力。一部分航天器用于处理水和挥发物，另一部分航天器用于收集金属。切割装置应使用无磁化材料，避免采集过程中黏附过多颗粒，降低挖掘装置的工作效率。同时，切割过程中探测器继续通过内部磁源线圈通电使小行星磁化，通过磁化收集端将一些碎石吸入存储装置中，防止污染太空环境。该装置还能避免速度较大的喷射物对探测器的损害。切割完成后装载矿石，而对于小行星更深层处的资源如硫、铀等，可采用地下流体开采的方法，即通过在挖掘矿石留下的切割缝隙中注入蒸汽或溶剂融化并提取液态硫等资源。完成上述操作后用激光烧结采集后的小行星表面，使其形成光洁墙面，令探测器能够识别该采集坑，便于日后继续在此处开采。最后减小电流使磁力达到探测器可与小行星分离的阈值，施加速度增量使其能转移至环绕轨道，完成一次采集任务。随后探测器可就近抵达另一个合适的采矿目标点，重复上述（3）~（6）过程。
（7）探测器及收割航天器返回近地轨道。待开采完成全部预测采矿点后，探测器及收割航天器将带着收集到的矿产返回近地轨道。
（8）在近地轨道进行矿物处理并运回地球。在近地轨道建立空间生产站和垃圾站，将所收集到的矿产进行处理，使用3D打印等技术，从小行星获取材料，然后制造出各种东西，最后运输至地球。对于那些用不到的废料，可以放置在垃圾站。
2.1.2 采矿的关键技术要求和设备说明
2.1.2.1 采矿的关键技术要求
采矿的关键技术要求包括以下6个方面。
(1) 在技术方案选择方面，采矿机器设备采用太阳能甚至是核能供电，且在任务中使用电推进或离子推进技术以减少从地球往返小行星所需的燃料。
(2) 在小行星采矿过程中，所有探测器和采矿设备须紧紧固定在小行星上，以防因失重而飘走迷失于太空。
(3) 在采矿位置选择的过程中，需要计算小行星表面的有效势，分析物体在小行星表面的运动趋势。
(4) 在锚定并固定收割航天器时，大多数小行星处于姿态旋转状态。收割航天器与小行星之间属于刚性连接，有着非常强的动力学耦合特性，爆破过程的稳定控制问题非常复杂。
(5) 在小行星上的原位资源利用，即用小行星上的资源，如水、冰来生产燃料和其他必需品，减少对地球补给的依赖。因为将水从地球带到太空的成本很高。所以对于水的处理可以在太空中形成专门的加工设施，处理小行星风化层和挥发物，如水和碳氢化合物，生产有价值的产品和氧气，用于生命维持和推进。此外，甲烷和甲醇是潜在的火箭燃料，可以与液氧搭配用于强大的火箭推进器。
(6) 航天器的自主导航、制导和控制一直是一个关键问题。在小行星采矿任务中，需要航天器能够进行自主导航、制导和控制，这涉及先进的控制系统和人工智能算法。
2.1.2.2 采矿设备说明
采矿设备的说明包括以下7个方面。
(1) 收割航天器需要始终锚定并固定在小行星上，其中锚定装置由冰螺栓、冷气推力器、鱼叉装置组合完成。冰螺栓依靠着陆器的冲击力刺入小行星表面，冷气推力器同时喷气反推，保证冰螺栓刺入小行星表面，然后发射鱼叉装置，形成对小行星表面的多点刺入和完全钻入风化层，实现收割航天器与小行星的固定。
(2) 探测器释放的软层薄膜袋材料要满足在极低温下可以进行工作，并且可以进行延展，有一定的弹性，抗冲击强度高，可以抵挡从小行星上脱落的岩石造成的意外撞击。还要具有耐腐蚀性。
(3) 小行星采矿装置应该相对较小，而不是庞然大物[20]，从而减少发射成本。因此，紧凑型的动力技术是必不可少的，可展开的太阳帆技术是太阳系内部操作的一种可能，太阳热推进也是一种可能。
(4) 采矿装置的兼容性。由于小行星特性不同，要满足安全性和可靠性。采矿装置的兼容性对于在独特和恶劣环境下的操作至关重要。小行星采矿装置的设计在面对各种各样的小行星大小、形状和旋转速度时应该是兼容的。特别地，以便除特殊情况外不需要定制设计，航天器应适应合理的尺寸和长宽比范围。
(5) 在爆破和切割过程中，会产生高温，所以装置应具有排热系统。
(6) 航天器及其设备配备有适当的防护以防止辐射损害的装置。
(7) 空间生产站主要是利用在小行星上发现的铁、镍和钴等材料，通过3D打印在太空中建造结构。这些结构将为未来的太空探索提供运输和供应支持。
2.2 小行星表面动力学特性
寻找合适的采矿位置是小行星表面开采的重要环节，为找到合适的采矿点，需计算小行星表面的有效势和临界摩擦系数，以此来判断小行星采矿装置和开采矿物、石块等物体在小行星表面的运动趋势。本节在引力势$U$ 的基础上加入离心项，可得到有效势$V$ ，根据得出的小行星表面有效势的相对大小分布即可判断小行星表面物体自由运动的趋势[21]。
小行星Bennu上拥有丰富的铁资源，具有较高的开采价值，故本文选择Bennu作为目标小行星。小行星Bennu形状为陀螺形，直径为492 m，重量约7.32×108 kg，密度约1190 kg/m³。将其模型划分为具有1348个顶点和2692个面的多面体模型。通过选取每一个面的中心位置点计算该处的有效势和临界摩擦系数等物理量，以便近似代表整个面上所有的动力学性质。在计算中采取归一化计算方法，归一化长度单位为小行星Bennu的等效半径246 m，归一化时间单位为其自转周期4.288 h[22]。
故此，本文对小行星Bennu表面每个面上的中心位置处受力$ P $ 的切向分量$ {P}_{{\mathrm{t}}} $ 进行计算，并用箭头表示在图2 中。由图2 可知，Bennu表面物体具有向赤道周围区域运动的趋势，这与Zhang等在文献[22]中对Bennu表面粒子的大规模运动仿真结论相同。
对于物体在小行星Bennu表面的自发运动情况，通过计算发现Bennu表面物体受到合力$ P $ 的法向量Pn均指向表面内侧，因此，不会产生自发起飞现象。对于自发滑动情况，本文计算了小行星Bennu表面临界摩擦系数${\mu _{\mathrm{c}}}$ [17]，其分布如图3 所示。当物体与小行星表面实际摩擦系数小于临界摩擦系数时，物体会发生自发滑动。由图3 可知，在Bennu腰部，临界摩擦系数较小，物体不容易发生滑动。在局部地形凸起处（图3 中黄色区域），临界摩擦系数较大，物体容易发生滑动。
通过分析小行星Bennu合力切向分量分布和表面临界摩擦系数分布，可以研究物体在小行星表面的运动趋势及自发运动行为，更精准地评估小行星表面的地质特征和物质分布，从而有效确定最优的采矿位置。这包括判断表面物体的运动状态、监测岩石的分布情况和寻找合适的爆破地点，以便在实际采矿过程中选择最具潜力和经济价值的区域。结合第2.1节小行星表面开采流程可对小行星采矿任务提供以下具体指导意见。
（1）从动力学角度出发，爆破地点和矿物收集装置都应尽可能选择在低纬度地区，以避免中高纬度地区表面物体自由运动趋势造成的矿物收集失败。
（2）如果在Bennu表面局部凸起处进行表面开采，需要对采矿装置进行更好地锚固，并更加注意对矿物的吸附收集，以避免采矿装置和矿物的自发滑动现象。
3 掠飞捕获
3.1 柔性绳网动力学模型构建
3.1.1 柔性绳网航天器动力学模型
柔性绳网航天器（FNS）由探测器和绳网组成，在此研究中，本文暂不考虑探测器的形状，因此，FNS动力学模型的准确性取决于模拟网络变形的程度。本文采用了经典的离散化Kelvin-Voigt[24]方法对柔性绳网进行建模。该柔性绳网被离散化为一些列节点，绳段的质量由集中在该绳段两端的节点代替，在相邻的绳段之间加入非线性弹簧阻尼原件体现该绳段的弯曲特性，忽略扭转特性。考虑到柔性绳具有拉压不对称的特性，且只可承受拉力，节点i和节点j之间的绳段的张力可描述为
其中$ {l}_{0}^{\left(p\right)} $ 是绳段i的初始长度，$ {k}_{{\mathrm{p}}} $ 是节点间螺纹的轴向刚度，$ {c}_{{\mathrm{p}}} $ 是节点间螺纹的阻尼系数。它们可以被定义为
式中，$ E $ 为柔性绳的弹性模量；$ \zeta $ 为柔性绳的阻尼比；$ {A_i} $ 为第$ i $ 个绳段的截面积；$ {\rho _{\mathrm{s}}} $ 为柔性绳的密度；$ {l_i} $ 为第$ i $ 个绳段的当前长度。
3.1.2 碰撞检测
研究柔性绳网在小行星表面的碰撞动力学，首要任务是检测绳网与小行星之间的接触情况。FNS触地后，小行星与地面立即产生接触力，因此，在模拟过程中，每一个时间步骤都需要检测它们是否发生了接触碰撞。利用第1节的球谐级数小行星参数化模型，可以快速对二者的碰撞情况作出判定。定义$ {d}_{i} $ 为柔性绳网中节点$ i $ 与小行星表面点的径向距离之间的关系，如式（13）所示。
式中，$ \left( {{\theta _i},{\varphi _i}} \right) $ 为节点$ i $ 的球坐标角度；$ R\left( {{\theta _i},{\varphi _i}} \right) $ 表示球坐标为$ \left( {{\theta _i},{\varphi _i}} \right) $ 的小行星表面点的绝对距离；$ {{\boldsymbol{r}}_i} $ 表示节点$ i $ 的位置矢量。当$ {d_i} \geqslant 0 $ 时，二者发生碰撞，且$ {d_i} $ 为节点$ i $ 在小行星表面的嵌入量；当$ {d_i} < 0 $ 时，二者未发生碰撞。
为了更方便地描述柔性绳网的碰撞动力学特性，本研究定义了函数$ {\delta _i} $ 表示柔性绳网中节点$ i $ 与小行星表面的接触碰撞情况，当$ {\delta _i} = 1 $ 时，二者发生碰撞，当$ {\delta _i} = 0 $ 时，二者未发生碰撞，即
3.1.3 法向接触力
接触力由法向接触力和切向接触力（摩擦力）组成，如果发生碰撞，这两者都会极大地改变柔性绳网的动态行为。不同于经典的恢复系数模型，弹簧阻尼碰撞模型将碰撞视为一个动态接触过程，即碰撞过程中的物理参数（例如速度、加速度、碰撞力等）随着变形量和变形速度的变化而变化，从而能准确地描述碰撞物体之间的接触过程。相比于恢复系数模型，弹簧阻尼模型的另一个优势在于能够提供碰撞物体之间的法向碰撞力，即将弹簧力和阻尼力视为两者的法向碰撞力，其中弹簧力反应了实际碰撞中的变形回复作用，阻尼力反映了实际碰撞中的能量耗散情况。鉴于这一优势，本文采用弹簧阻尼碰撞模型研究柔性绳网与小行星的法向碰撞过程。
式中，$ {d}_{i} $ 为节点$ i $ 在小行星表面的嵌入量；$ {\dot{d}}_{i} $ 为节点$ i $ 在小行星表面的法向嵌入速度；$ {\delta }_{i} $ 描述了柔性绳网节点与小行星表面的接触情况；$ \alpha $ 为指数，$ \alpha \geqslant 1 $ ；$ {K}_{n} $ 为修正弹性指数，由Hertz碰撞理论推导得
式中，$ {\mathit{\text{κ}}}^{\mathrm{*}} $ 为等效曲率半径，$ {\mathit{\text{κ}}}^{\mathrm{*}}={\left(1/{\text{κ}}_{1}+1/{\text{κ}}_{2}\right)}^{-1/2} $ ，$ {\text{κ}}_{1} $ 为曲率半径，本文中取较小定值代替柔性网绳节点曲率半径，$ {\text{κ}}_{2} $ 为小行星表面在碰撞点处的曲率半径；$ {E}^{\mathrm{*}} $ 为等效弹性模量，$ {E}^{\mathrm{*}}={\left(\left(1-{v}_{1}^{2}\right)/{E}_{1}+\left(1-{v}_{2}^{2}\right)/{E}_{2}\right)}^{-1} $ ；$ {E}_{1} $ 为柔性绳网的弹性模量；$ {E}_{2} $ 为小行星表面物质的弹性模量；$ v_2 $ 为小行星$ v_1 $ 表面物质的泊松比；$ v_1 $ 为柔性绳网的泊松比。$ {C}_{n} $ 为修正阻尼系数，根据Hunt-Crossley阻尼理论可推导得
式中，$ e $ 为恢复系数；$ {v}_{0} $ 为柔性绳网节点刚与小行星表面接触时的速度大小。
3.1.4 切向接触力
当柔性绳网节点与小行星表面运动时，必然会受到摩擦力的影响。绳网节点可能处于黏滞状态，也可能处于滑动状态，而摩擦力随着运动状态的不同而不同。本文结合库伦摩擦定律和颗粒动力学理论进行考虑，为了计算各种运动下的摩擦力，需要建立两个动力学模型，当节点在小行星上滑动时，可以推导出切向摩擦力，而处于黏滞状态时，则需要引入弹簧阻尼器来描述切向接触模型。
基于柔性绳网处于黏滞或滑动的判断条件，不同状态下的切向摩擦力的计算表达式如式（18）所示。
式中，$ {k}_{i} $ 为黏滞状态下的切向拉伸弹簧刚度系数；$ {\mathit{s}}_{i} $ 为切向拉伸弹簧当前形变；$ c $ 为黏滞状态下的切向阻尼器的阻尼系数；$ {\mathit{v}}_{i}^{\tau } $ 为t时刻柔性绳网节点$ i $ 在小行星表面的单位切向速度；$ \mu $ 为滑动状态下的动摩擦系数；$ {t}_{0} $ 为初试碰撞时刻，$ t $ 为当前碰撞时刻。
3.1.5 柔性绳网的动力学方程
根据前文的分析，可以明确柔性绳网上各节点的受力情况。本文使用Kelvin-Voigt方法将柔性绳网离散为一些列质点，因此可以通过组合各个节点的动力学方程来获得柔性绳网的动力学微分方程组，如式（20）~式（22）所示。
式中，$ {m}_{i} $ 为绳网节点$ i $ 的质量；$ {\mathit{r}}_{i} $ 为绳网节点$ i $ 当前的空间位移；$ {\mathit{F}}_{{\mathrm{ex}}}^{i} $ 为绳网节点$ i $ 受到的外部作用力；$ {{F}}_{{\mathrm{in}}}^{i} $ 为绳网节点$ i $ 受到的绳网内部连接作用力；$ Np $ 为柔性绳网中节点总数；$ \omega $ 为小行星绕自身惯性主轴旋转的角速度；$ \nabla U\left({\boldsymbol{r}}_{i}\right) $ 为绳网节点$ i $ 受到的小行星引力；$ {N}_{i} $ 为绳网节点$ i $ 受到的支持力大小；$ {\boldsymbol{n}}_{i} $ 为绳网节点$ i $ 所处位置的单位法向量；$ {\mathit{f}}_{i} $ 为绳网节点$ i $ 受到的库伦摩擦力；$ {\delta }_{i} $ 为碰撞函数；$ Ne $ 为柔性绳网中与节点$ i $ 邻接的绳段总数；$ Nb $ 为柔性绳网中与节点$ i $ 邻接的弹簧总数；$ {\boldsymbol{T}}_{l}^{\left(p\right)} $ 为柔性绳网中与节点$ i $ 邻接绳段$ p $ 施加的轴向弹簧拉力；$ {\boldsymbol{D}}_{l}^{\left(p\right)} $ 为柔性绳网中与节点$ i $ 邻接绳段$ p $ 施加的阻尼力；$ {\boldsymbol{T}}_{b}^{\left(q\right)} $ 为柔性绳网中表征弯曲作用的弹簧$ q $ 施加的压力；$ {\boldsymbol{D}}_{b}^{\left(q\right)} $ 为柔性绳网中表征弯曲作用的阻尼$ q $ 施加的阻尼力。
通过计算柔性绳网在小行星引力场中的动力学方程，可得到柔性绳网在小行星引力场中的空间位移和运动速度等信息。前文推导了柔性绳网的微分方程组，但是柔性绳网离散之后的节点数量较多，未知数较多，计算量极大。针对这一问题，本文拟用基于四阶龙格-库塔法的一种求解动力学方程组的高效数值算法SnOdeEule算法进行数值求解[28]。
3.2 柔性绳网的运动仿真
3.2.1 模型参数确定
本节模拟了柔性绳网捕获小行星的情景。柔性绳网构型如图4 所示[28]，绳网尺寸为1064 m×90 m×1061m，节点总数为712，采用凯拉夫材质，在y轴方向上距离小行星50 m。
为更直观地表示绳网在小行星引力场中的运动情况，本文从柔性绳网上选取了部分节点，节点坐标如表1 所示，节点位置如图5 所示。
3.2.2 结果分析
4 结语
在地球资源日益枯竭的背景下，小行星采矿被认为是未来获取矿物资源的重要途径。本文设计了小行星采矿方案并研究了采矿过程中的关键动力学问题，通过采用多面体法建立小行星引力场模型，引入了参数曲面数学模型；提出了一种表面开采的高效采矿方案，并对其关键技术和设备说明作出了详细分析，选取Bennu为目标小行星，结合表面开采方案研究了Bennu小行星表面动力学特性；研究了柔性绳网的碰撞动力学，仿真了柔性绳网捕获小行星Bennu的动力学过程。