宇宙尘埃轨迹探测器信号提取方法及反演精度影响因素分析

冯桃君，焦子龙,2，姜利祥,2，郑慧奇，姜海富,2，刘宇明,2，王 鹏，刘学超，李 涛

（1.北京卫星环境工程研究所； 2.可靠性与环境工程技术重点实验室；3.航天东方红卫星有限公司：北京 100094； 4.中国航天科技集团有限公司，北京 100048）

0 引言

微米和亚微米级的宇宙尘埃粒子通常来自彗星、小行星、月球、星体间的碰撞碎片和行星际物质[1-2]。宇宙尘埃暴露在等离子体、太阳紫外辐射、高能粒子、电磁场等组成的复杂空间环境中，运动过程中会携带电荷，因此尘埃轨迹是行星和星际磁场、太阳和行星引力、等离子体阻力及辐射压力等因素共同作用的结果[3-4]。Xie 等[5]综述并研究宇宙尘埃的运动方式发现：小行星尘埃由于Poynting-Robertson 效应会以螺旋的方式缓慢接近太阳；β-流星体主要受到太阳辐射压的作用，以双曲线轨迹逃逸出太阳系；Ulysses 号航天器上的尘埃探测器发现木星尘埃粒子受到磁层电场的加速并以大于100 km/s 的速度逃逸出木星。并且指出，宇宙尘埃的轨迹和电荷信息有助于确定粒子起源，辨别不同起源的粒子，揭示带电尘埃粒子与行星际空间环境的相互作用。

对宇宙尘埃粒子进行探测的方法很多是利用尘埃探测器的原位测量[6-8]，高速粒子在探测器内部发生碰撞电离，可通过电离产生的总电荷量和电脉冲信号的上升时间分别推算出粒子质量和碰撞速度；并通过对粒子飞行时间的分析得出粒子的化学成分。该项技术已在Vega、Wind、Cassini、STEREO等航天器上用于探测尘埃的质量、成分和速度大小[9-11]，但速度方向只能通过统计数据分析得到，要求探测器进行全方位扫描，限制了碰撞探测器的应用。Grün 等[12]提出了尘埃天文学概念，涉及宇宙尘埃粒子的起源以及粒子化学成分、电荷、速度、轨迹的测量。尘埃轨迹探测器（dust trajectory sensor,DTS）是尘埃天文学的关键探测组件，当带电粒子穿过DTS 内部的金属丝阵列时，会产生诱导电荷，可通过分析金属丝输出的电信号提取粒子的电荷量和轨迹信息[1,13]。DTS 与不同的仪器结合可构成尘埃望远镜[14]、离子质量分析器[15]和静电质量探测器[1]等，适用于不同速度、不同尺寸的尘埃颗粒探测。

DTS 在几十年前就已被提出[5]，其几何构型已经过不断优化。Auer 等[1-2]利用COULOMB 软件对不同几何构型的DTS 进行数值仿真，分析金属丝的几何参数、间距及数量等对电荷分布、信号灵敏度的影响，获得了适用于速度在2～5 km/s 的微米级尘埃探测的DTS 结构参数；并得出DTS 的电荷噪声比（charge to noise ratio, QNR）≥6.25 才能从背景噪声中探测到尘埃信号。Xie 等[5]同样利用COULOMB 软件，发现探测器的QNR 越小电荷量探测误差分布范围越大，当QNR=10 时，尘埃电荷量探测误差为-1.5%～5%。为提高仪器对尘埃电荷量的探测灵敏度，Li 等[16]提出了分段式低电容金属丝阵列结构，使尘埃电荷量探测阈值降至0.2 fC。Voronov 等[17]介绍了一种基于网格面和金属丝阵列面交替的DTS 设计，不仅可以测量带电微流星体的电荷量、轨迹，还能计算微流星体的质量。

宇宙尘埃大都分布在亚微米到厘米量级，速度在每秒几千米到几十千米不等，DTS 的探测精度也因尘埃性质不同而显示差异。本文基于典型的DTS结构，针对尘埃轨迹和电荷量的测量建立DTS数值仿真模型，研究从DTS 金属丝信号提取尘埃粒子速度矢量和电荷量的方法，重点分析尘埃速度、带电量和入射方向对算法精度的影响，有助于探索适用于DTS 探测的尘埃范围。

1 DTS 结构及工作原理

DTS 一般由中间4 个金属丝阵列面和分别位于顶部和底部的2 个屏蔽栅网组成，每个金属丝阵列面包含7～16 条相互平行的金属丝，两相邻平面的金属丝方向正交。如图1[5]所示，两屏蔽栅网接地；每条金属丝都独自与一个电荷灵敏放大器（charge sensitive amplifier, CSA）连接，CSA 通道的输出信号被瞬态记录器以一定频率采集。当带电尘埃粒子穿过DTS 时，会产生诱导电荷。根据同一平面金属丝信号特性可确定粒子穿过当前平面的一维坐标（x1、y2、x3、y4），如图2[5]所示，进而可推导出尘埃粒子的速度矢量和电荷量[1,5]。

图1 DTS 结构示意[5]Fig.1 Structure schematic of DTS[5]

图2 带电尘埃粒子在DTS 中的穿越轨迹示意[5]Fig.2 Trajectory schematic of charged particles travelling in DTS[5]

本文研究所用DTS 模型含4 个平行于xy平面的金属丝阵列面，每个平面包含7 条金属丝，相邻平面间隔40 mm，金属丝直径0.4 mm，长140 mm，间隔20 mm； DTS 的整体尺寸为160 mm×160 mm×200 mm（长×宽×高），如图3 所示。平面编号从上到下依次增加，平面1、3 的金属丝平行于y轴，平面2、4 的金属丝平行于x轴。模型坐标原点位于结构的几何中心，图中以红点标记。

图3 DTS 数值模型Fig.3 Numerical model of DTS

探测器金属丝的诱导电荷量与尘埃粒子的入射位置、角度、探测器尺寸和构型有关。Auer 等[1]给出了当带电量为Q的球形尘埃颗粒穿过DTS时，金属丝i的诱导电荷qi的近似表达式

式中：ri为尘埃到金属丝i的垂直距离；j(1/rj)代表求和遍历了DTS 所有的金属丝。尘埃的位置(xp,yp,zp)随时间变化，因此尘埃到金属丝的距离和金属丝的诱导电荷均为时间的函数。尘埃电荷量可由所有金属丝诱导电荷相加近似得到。

Auer 等[1]在地面试验中采集到的随时间变化的DTS 输出信号如图4 所示，图中信号来自两金属丝阵列平面，每阵列包含16 条金属丝。输出信号的形状反映金属丝与尘埃粒子的距离，信号强度反映金属丝带电量，信号峰值表明尘埃粒子近距离地从某金属丝旁穿过该金属丝所在平面。

图4 Auer 在尘埃实验室记录到的DTS 信号[1]Fig.4 DTS signals recorded in the dust laboratory by Auer[1]

假设高速尘埃粒子在DTS 内的轨迹为一条直线，由尘埃粒子进、出DTS 的时间（tin、tout）和位置坐标（(xin, yin, zin)、(xout, yout, zout)）可以确定尘埃粒子的速度和方向。通过式(2)可由4 个一维坐标确定尘埃粒子进、出DTS 的坐标(xin, yin)、(xout, yout) [5]：

粒子速度的大小为[5]

粒子入射角为[5]

其中，θx、θy分别为尘埃轨迹在xz、yz平面的投影与z轴的夹角。

2 尘埃参数提取方法

2.1 计算流程

为提取尘埃粒子的速度矢量和电荷量，在已知探测器几何构型的前提下，表征尘埃轨迹和电荷量需要7 个独立参量：尘埃电荷量（Q）；尘埃进、出DTS 的时刻（tin、tout）；尘埃穿过4 个金属丝阵列平面的一维坐标（x1、y2、x3、y4），故设参数矢量为P=(Q,tin,tout,x1,y2,x3,y4)。文献[5]提出了一种基于探测数据与仿真数据拟合的参数提取方法，其反演流程如图5 所示。

图5 参数提取反演流程Fig.5 Inversion process for parameter extraction

参数提取的具体步骤为：

1）初始化参数矢量P0。尘埃初始电荷量Q可设置为DTS 所有金属丝电荷量之和；x1、y2、x3、y4根据式(1)估算得到。文献[5]中并未介绍tin和tout的取值方法，本文设DTS 的采样频率为r（单位为MS/s，即106sample/s），对应的信号采样个数为n，并将第1 个采样信号的时刻t1设置为0 作为参考，则tin∈[-1/r, 0]、tout∈[(n-1)/r,n/r]，单位为μs，tin和tout的初值在各自范围内均匀随机取值。

2）初始化参数的误差范围ΔP0。ΔQ为±10%，Δt为±1/r，Δx或Δy为±10 mm。

3）数据拟合。从P0±ΔP0的范围中随机选取一组P=(Q,tin,tout,x1,y2,x3,y4)，根据式(1)以及DTS采样频率仿真计算DTS 各金属丝的采样信号，然后计算仿真信号相对实测信号的误差，

式中：χ2为仿真信号误差；N为数据总个数；Dm为第m个实测数据；Sm为第m个仿真数据。

重复上述随机选取参数P即χ2计算103次，χ2值随P变化，χ2最小值对应的参数集为P1,min。

4）优化参数矢量P1。由于χ2也是单一参数的函数，依次改变P1,min中第l（l=1, 2, ···, 7）个参数P(l)的值，使χ2最小的参数值成为该参数的新估计值P1(l)。

5）缩小参数的误差范围ΔP1。ΔP1取2ΔP0/3和3 |P1-P0|中的大者。

6）重复步骤3）～5），直到χ2的变化小于0.1%。

2.2 数值仿真与结果分析

为验证2.1 节尘埃参数提取方法的有效性，针对典型算例进行计算分析。尘埃预设轨迹为入射点坐标(13, -22, 100)、出射点坐标(-7, 38, -100)，入射角θx=-5.7°、θy=16.7°，z方向速度分量vz=5 km/s，尘埃穿过各金属丝平面的一维坐标为x1=9 mm，y2=2 mm，x3=1 mm，y4=26 mm，涵盖了尘埃与金属丝之间的近距离、中等距离和远距离。另外，金属丝输出信号用诱导电荷与尘埃电荷之比表示，设尘埃电荷为1 C。设置QNR=10，采样频率为r=10 MS/s，尘埃穿过DTS 的时间为t=40 μs，因此每个金属丝有n=400 个采样信号。设尘埃进入DTS 的时间tin=0，则tout=40 μs，参数真值Ptrue=(1, 0, 40, 9, 2, 1,26)，本文用下标“opt”和“true”分别表示参数的最优估计值和真值。离尘埃近的金属丝将聚集大部分诱导电荷，因此在信号分析中，为节省计算时间只对每个平面中离尘埃最近的2 条金属丝的信号进行拟合，这样用于拟合的信号一共有3200 个。实际的观测信号包含了电荷灵敏放大器产生的随机噪声，为模拟真实信号，在式(1)计算出的电荷信号中添加正态分布的随机噪声，其均值为μ=0，标准差σ=1/QNR。根据上述仿真条件，随尘埃位置zp变化的无噪声金属丝模拟信号如图6(a)所示，图中虚线表示各平面对应的z轴坐标，顶端数字为平面序号。从图中可看出，当尘埃穿过平面时与金属丝距离很近（如y2=2 mm，x3=1 mm），该金属丝信号幅值骤升，随后下降；与此同时，尘埃另一侧相距较远的金属丝信号出现下凹。当尘埃穿过平面时与两侧金属丝距离相差不大（如x1=9 mm），两金属丝信号形状相似。图6(b)所示是叠加噪声后的金属丝信号，用于模拟DTS 的实测信号。

图6 尘埃粒子从(13, -22, 100)到(-7, 38, -100)穿越DTS的模拟信号Fig.6 Simulated signals generated by dust particles travelling from (13, -22, 100) to (-7, 38, -100) through DTS

按上述尘埃轨迹与电荷信号的分析方法，经过6 次迭代得到各参数的最优估计Popt=(0.99, -0.05,41.27, 9.04, 2.06, 0.92, 26.07)，χ2随迭代次数的变化如图7 所示。参数估计偏差ΔP=Popt-Ptrue=(-0.01,-0.05, 1.27, 0.04, 0.06, -0.08, 0.07)，尘埃电荷和轨迹的反演结果如表1 所示。其中，速度的反演结果误差最大，被低估了7.63%，其原因可能是在独立参数的估计中算法能获得精度较高的尘埃位置，但尘埃穿过DTS 的时间估值偏差较大（1.32 μs），导致反演的尘埃速度减小。

表1 尘埃电荷与轨迹的反演结果Table 1 Inversion results of the dust charge and trajectory

图7 χ2 随迭代次数的变化Fig.7 Variation of χ2 with iteration steps

参数初值P0与最优估计Popt对应的金属丝信号分别如图8(a)和(b)所示。

图8 不同参数估计对应的金属丝信号Fig.8 Metal wire signals corresponding to different parameter estimates

图中蓝色表示模拟的实际观测数据，红色表示根据P0和Popt反演出的金属丝信号。可以看到，经过6 次迭代后，Popt能复现实际的观测信号，说明时间估计的误差只影响尘埃速度大小的反演结果，对尘埃速度方向及DTS 金属丝信号的反演结果没有影响。

3 反演精度的影响因素分析

分别改变用于仿真的尘埃速度分量vz、入射角度θy、尘埃电荷量Q和探测器电荷噪声比QNR，分析它们对反演精度的影响。各参数取值如表2 所示，当改变其中1 个参量时，其他3 个参量按表2中粗体数字取值。根据采用的DTS 结构，尘埃穿过DTS 的最大入射角为38.6°，仿真设置θx=0°，θy在0°～30°范围内变化，表3 列出了与表2 中θy对应的预设尘埃轨迹（进、出DTS 的位置坐标）。

表2 各控制变量的取值Table 2 Values for control variables

表3 不同θy 对应的尘埃轨迹Table 3 Dust trajectories for different values of θy

按表2 依次改变各控制变量数值，针对每一组控制变量（vz、θy、Q、QNR），重复图5 的反演流程100 次，每次的模拟信号叠加不同的噪声，会反演得到100 组不同的尘埃电荷量（Q）及轨迹参数（vz、θy、θx）。尘埃入射角度的反演误差Δθ为反演最优值与真值之差，即

电荷量和速度大小的反演误差Δa以相对误差表征，即

其中符号a代表参数Q或v。

计算得到100 组反演结果的反演误差后，分析这100 组反演误差的均值及误差绝对值的均值随表2 中各控制变量的变化，结果如图9～图12 所示。其中，误差绝对值的均值（图中蓝线）反映反演结果相对 “真”值的整体偏离幅度，误差均值（图中粉线）反映反演结果相对“真”值的整体偏离方向（被高估或被低估）。

图9 反演误差均值及其绝对值均值随vz 的变化Fig.9 Variations of mean inversion error and mean absoluteinversion error with vz

3.1 控制变量为vz 的情况

反演误差随vz的变化如图9 所示。

3.2 控制变量为轨迹倾角θy 的情况

反演误差随θy的变化如图10 所示。

图10 反演误差均值及其绝对值均值随θy 的变化Fig.10 Variations of mean inversion error and mean absolute inversion error with θy

3.3 控制变量为尘埃电荷量Q 的情况

反演误差随Q的变化如图11 所示。

图11 反演误差均值及其绝对值均值随Q 的变化Fig.11 Variations of mean inversion error and mean absolute inversion error with Q

从图11 可以看出，尘埃电荷量的变化对各参数反演误差没有显著影响。保持在0.4%左右，在0.2%左右，和接近0.15°。在反演结果偏离方向上，尘埃电荷量和速度大小整体被低估0.1%左右，轨迹倾角相对于“真”值同样没有明显的偏移方向。

3.4 控制变量为QNR 的情况

反演误差随QNR 的变化如图12 所示。

图12 反演误差均值及其绝对值均值随QNR 的变化Fig.12 Variations of mean inversion error and mean absolute inversion error with QNR

通过上述分析可知，本文所提出的方法在其考虑的参数范围内，对尘埃电荷量的平均反演精度优于1.5%，速度大小的平均反演精度优于0.8%，轨迹倾角的平均反演精度优于0.6°。

4 结束语

尘埃穿过DTS 会诱导金属丝产生电信号，通过对观测信号和仿真信号进行拟合可提取尘埃的电荷量与轨迹信息。本文对DTS 探测尘埃速度矢量进行了数值仿真，结果表明：

1）尘埃电荷量和轨迹的反演误差随尘埃速度vz的增大而增大，随探测器QNR 的增大而减小，基本不受轨迹倾角变化和尘埃带电量的影响。

2）尘埃电荷量反演精度优于1.5%，速度大小反演精度优于0.8%，速度方向反演精度优于0.6°。

本研究通过DTS 探测获得了高精度的尘埃电荷量和轨迹信息，并对探测精度影响因素进行了仿真分析，可为DTS 研制与科学探测提供参考。需要说明的是：由于反演过程需要7 个独立参数，本文算法中通过随机计算103次来确定χ2的最小值很大可能是局部最小而不是全局最小；本文用于仿真的DTS 结构为较优结构，并未针对其他结构类型（如不同的栅网平面间距、金属丝间距等）进行研究；本文对尘埃颗粒形状进行了简化（以球形代替），且用于计算探测器信号的公式与实际情况间存在一定偏差。后续将重点开展金属丝诱导电荷信号的理论和试验研究，依据更符合实际情况的信号修正本文的计算方法，进一步开展反演精度分析。