本发明涉及激光测距中的探测效果分析与仿真技术，特别涉及一种回波数量低、回波展宽大的空间碎片激光测距效果仿真分析方法和装置。
背景技术
激光测距技术是一种高精度距离测量方式，通过测量激光脉冲在目标与探测器之间的飞行时间获得目标距离，精度优于雷达一至两个数量级，在空间目标测量领域中具有广泛应用。对激光测距效果进行仿真分析，对获得先验测距结果及评估实测结果具有重要意义。第一，可以获得回波的散点分布仿真图，便于直观分析，选择后续回波识别处理的方法；第二，可以根据预估的回波和噪声强度获得探测概率和虚警概率，从而设计探测器的响应阈值；第三，可依据实测结果与仿真结果的差异，排查测距过程可能出现的问题。
空间碎片无角反射器，存在反射回波数量低、回波展宽大的特点。空间碎片激光测距如果采用传统的基于划分区间计算探测概率的仿真分析方法，会因为过于依赖区间划分长度，出现仿真结果与实测结果不一致的虚假周期性现象，使仿真结果失去参考价值。因此需要研究建立一种针对空间碎片激光测距的有效的仿真分析方法，提高仿真结果与实测结果的一致性。
技术实现要素：
有鉴于此，本发明提供了一种空间碎片激光测距效果仿真分析方法和装置，能够在无需对光子进行区间划分的前提下，完成对空间碎片激光测距效果的仿真分析，提高仿真效果与实测效果的一致性。
为了解决上述技术问题，本发明是这样实现的。
一种空间碎片激光测距效果仿真分析方法，该方法包括：
针对光子时刻序列，按照光子先后顺序，判断光子是否被探测；探测方式为：在0-1的均匀分布范围内，生成随机数δ，判断随机数δ与单个光子探测效率η＝1-e-1的大小关系；如果δ≤η，那么该光子被成功探测，如果δ>η，那么该光子未被探测，则继续判断后一个光子，直到所有光子判断完成；
基于光子探测结果，生成仿真的探测器输出信号。
优选地，在按照顺序判断光子是否被探测时，若当前光子判定为被探测，则不再探测探测器死时间内到达的光子，直接探测死时间之后的光子；若当前光子不被探测，则继续判断后一个光子，直到所有光子判断完成。
优选地，所述基于光子探测结果，生成仿真的输出信号为：在成功探测的光子上叠加探测器时间抖动的影响，生成仿真的输出信号。
优选地，所述叠加探测器时间抖动的影响为：认为探测器时间抖动符合高斯分布，随机生成符合时间抖动分布的随机数δ，那么仿真的探测器输出信号序列为lx＝rxi+δ，其中rxi为到达探测器的光子序列中的第i个光子的光子本身距离，x表示该光子是探测器成功探测到的第x个光子。
优选地，该方法进一步包括：将仿真的探测器输出信号进行显示，获得空间碎片激光测距的散点分布仿真图。
优选地，所述光子时刻序列由目标反射光子和噪声光子组成，并对这些光子按照时刻进行排序。
优选地，该方法进一步包括：按照以下公式计算探测过程的探测概率pd和虚警概率pfalse；
pd＝e-w·(1-e-s-n)
pfalse＝1-e-w+e-w-s-n-e-w-s-n-q
其中，w为回波分布范围之前的噪声光子数，s为经空间碎片反射后的回波光子数，n为回波分布范围内的噪声光子数，q为回波分布范围之后的噪声光子数。
优选地，该方法进一步包括：改变经空间碎片反射后的回波光子数s和噪声强度vnoi，得到探测概率和虚警概率随二者的变化曲线图；其中，噪声强度vnoi的改变会影响q、n和w的值。
本发明还提供了一种空间碎片激光测距效果仿真分析装置，该装置包括探测模块和仿真模块；
所述探测模块，用于针对光子时刻序列，按照光子先后顺序，判断光子是否被探测；探测方式为：在0-1的均匀分布范围内，生成随机数δ，判断随机数δ与单个光子探测效率η＝1-e-1的大小关系；如果δ≤η，那么该光子被成功探测，直到所有光子判断完成；
所述仿真模块，用于基于光子探测结果，生成仿真的探测器输出信号。
优选地，该装置进一步包括概率计算模块，用于按照公式(8)和(9)计算探测过程的探测概率pd和虚警概率pfalse；
pd＝e-w·(1-e-s-n)(8)
pfalse＝1-e-w+e-w-s-n-e-w-s-n-q(9)
其中，w为回波分布范围之前的噪声光子数，s为经空间碎片反射后的回波光子数，n为回波分布范围内的噪声光子数，q为回波分布范围之后的噪声光子数。
有益效果：
(1)本发明在考察每个光子是否被探测时，通过引入随机数来对比该光子探测概率大小以判断光子被探测与否，无需强制探测器进行不合物理实际的“光子累积”，无需对光子进行复杂的区间划分，只需对光子简单排序并进行探测概率计算，但却能反映更真实的探测结果。
(2)本发明在判断某个光子被成功探测后，略过对探测器死时间内到达的光子进行判断的过程，符合探测器探测过程中需要一定死时间来恢复到待测状态的物理特性，并且相比于对每个光子判断后再剔除死时间内被探测的光子的方式，具有更高的仿真效率，可简化操作、节约时间。
(3)本发明考虑了探测器时间抖动随机数带来的光子输出响应差异，能够体现输出过程中探测器响应时间带来的统计性差异，更符合探测过程物理本质。
(4)本发明在对光子进行排序时，将回波分布情况假定为高斯分布，而将噪声分布情况假定为均匀分布，极大地简化了仿真步骤，却能够在较大程度上反映光子分布的实际情况。
(5)本发明可通过对统一的光子序列中每个光子探测概率的计算，得到测距过程的探测概率和虚警概率，从而为实验过程中探测器响应阈值的设计、数据处理算法的选择提供参考，也能评价整个激光测距过程的难易度。
总之，本发明综合考虑了探测器本身的即时响应的物理特性，避免了人为设立区间进行光子的强制累积和分段，更符合物理过程实际，并且能够真实客观体现探测概率/虚警概率随回波和噪声的变化，得到的测距结果散点分布仿真图，有助于后续设计探测器响应阈值，方便选择合适的数据处理方法，以及排查测距过程可能出现的问题，是一种灵敏度高、可信度强的空间碎片激光测距效果仿真分析方法。
附图说明
图1为探测概率与虚警概率随回波和噪声强度的变化关系仿真图。
图2为统一的光子序列{ri}示意图。
图3为空间碎片no.24809探测光子的散点分布图；其中(a)为仿真结果；(b)为实测结果。
图4为本发明操作流程示意图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本发明提供了一种空间碎片激光测距效果仿真分析方法，该方法的前期分析和实际使用过程如下：
步骤一、将空间碎片材料反射率、距离、大气透射率等参数代入空间碎片回波光子数方程，获得回波光子数，并估计回波光子分布；依据天光背景亮度、探测器暗噪声等参数估计噪声数量及分布。
其中，认为回波分布服从高斯分布，噪声分布服从均匀分布。空间碎片等非合作目标仅依靠其表面漫反射特性反射激光，其平均回波光子数方程为
其中，λ为激光波长，h为planck常数，c为光速，et为激光发射能量，t为大气透过率，kt，kr分别为发射系统和接收系统光学效率，d为望远镜接收口径，θt为入射激光发散角，r为空间碎片距离，s为空间碎片等效散射截面积，r为空间碎片等效反射系数。例如，某空间碎片平均回波光子数为0.1的统计意义是，每发射10个脉冲平均才有一个回波光子到达探测器。
步骤二、认为探测器类似具有传递函数功能的“黑箱子”，对到达光敏面的光子进行“即时判断”，以考察其是否被探测。若该光子被探测，那么在探测器死时间内不再探测其他光子，即跨过死时间的时长，直接探测死时间之后的光子；若该光子不被探测，那么探测器继续判断下一光子是否被探测，整个过程不等待光子累积。依次按照光子序列顺序计算第x个光子被探测到的概率px为：
px＝e-(x-1)·(1-e-1)(2)
步骤三、在单个光子被探测概率的基础上，计算得到激光测距的探测概率pd和虚警概率pfalse，并得到其与回波和噪声强度的变化关系图，如附图1，可以发现空间碎片的探测概率随回波光子数增加以及噪声光子数降低而增大，虚警概率则反之。
特别需要注意的，在输出信号中，处于回波展宽范围内的回波信号和噪声相互不可分辨，这些光子获得的响应概率之和即为回波探测概率。
其中，探测概率pd和虚警概率pfalse的推导方式为：
第1步：计算回波展宽区间内每个光子的探测概率为：
其中经空间碎片反射后的回波光子数为s，假设噪声速率为vnoi，回波展宽时长为techo，门控时长为tgate，回波在门控中的位置为γ，取值为0～1。那么回波分布范围内的噪声数为n＝vnoi·techo，回波分布范围内所有不可分辨的回波和噪声的光子个数之和为(s+n)，回波分布范围之前的噪声光子数为w＝vnoi·(tgate-techo)·γ，p1内，p2内，…，p(s+n)分别为回波分布范围内第1个，第2个，…，第(s+n)个光子的探测概率。
第2步：将第1步中得到的概率相加，得到探测概率pd为：
第3步：计算回波分布范围之外的所有噪声各自的探测概率为：
其中q＝vnoi·(tgate-techo)·(1-γ)为回波分布范围之后的噪声光子数。
第4步：将第3步中得到的概率相加，得到虚警概率pfalse为：
第5步：改变回波光子数s和噪声强度vnoi,便可得到探测概率和虚警概率随二者的变化曲线图。
步骤四、综合考虑所需判断的激光测距过程中的发射激光脉冲总数，依据每个脉冲的目标反射光子数s和噪声光子数n乘以脉冲数量，得到待探测光子集合以供后续散点分布图仿真；将回波光子和噪声光子按照时间序列进行排序，得到统一的光子时刻(距离)序列{ri}。在统一的光子序列中，按照光子先后顺序，判断光子是否被探测。若被探测，则探测器不再探测其死时间内到达的光子，若不被探测，则继续判断后一个光子，直到所有光子判断完成。而后叠加探测器时间抖动的影响，生成仿真的输出信号。在实际中，也可以每一个光子都进行判断，最后再去除死时间内的光子数据，这样的计算效率不如前者。
其中，光子是否被探测的具体判断方式为：
第1步：在0-1的均匀分布范围内，生成随机数δ。
第2步：判断随机数δ与单个光子探测效率η＝1-e-1的大小关系。如果δ≤η，那么该光子被成功探测，在探测器死时间内，不再探测后一个光子。如果δ>η，那么该光子未被探测，则继续判断后一个光子。
本步骤使用的公式η＝1-e-1由式(2)变形而来，比较二式，去掉了式(2)的e-(x-1)一项。这是因为，在实际仿真过程中，但凡需要对第i个光子进行判断，说明其前面的光子未被成功探测，所以式(2)中表征前序光子未被探测概率的e-(x-1)一项的概率已确定为1。
第3步：重复第1步和第2步，直到所有光子遍历完成。
之后，生成散点分布仿真图的具体方式为：
认为探测器时间抖动符合高斯分布，随机生成符合时间抖动分布的随机数δ，那么探测器输出信号序列为lx＝rxi+δ，其中rxi为到达探测器的光子序列中的第i个光子的光子本身距离，x表示该光子是探测器成功探测到的第x个光子。可以发现探测器输出信号完全由探测到的具体光子距离以及探测器时间抖动决定。将输出信号进行显示，可获得空间碎片激光测距的散点分布仿真图。
在实际应用中，如果希望获得空间碎片激光测距的散点分布仿真图，则可以采用步骤四；如果希望获得探测过程的探测概率pd和虚警概率pfalse，则直接利用公式(8)和(9)，进而可以获得探测概率和虚警概率随回波光子数s和噪声强度vnoi的变化曲线图。前文的步骤一～步骤三是推导过程。
为了更加清晰在实际仿真分析时执行哪些步骤，下面结合图4给出如下执行步骤：
步骤一、评估空间碎片激光测距中回波和噪声分布情况，首先估计回波展宽和强度，得到回波光子的高斯分布情况；其次估计噪声强度，结合天光背景噪声、探测器暗噪声等，得到噪声光子的均匀分布情况；然后结合回波和噪声各自分布，以及测量过程中经历的激光脉冲个数，得到统一的光子分布序列{ri}。
步骤二、根据统一的光子分布序列，判断每个光子的探测概率与生成随机数的大小关系，叠加符合探测器时间抖动分布的时间抖动随机数，得到探测器对所探测光子的实际响应距离，得到探测光子的散点分布仿真图。根据散点分布仿真图，可以评估空间碎片激光测距效果，回波越强，噪声越弱，表明响应阈值可以设置更高，空间碎片越容易被识别。
步骤三、进一步的，根据本发明推导的探测概率和虚警概率的物理表达式，将步骤一评估获得的回波和噪声强度代入，得到探测概率与虚警概率的相应值，可供测距效果评估及探测器响应阈值设立。
下面结合数据举一个实例。
本发明的仿真对象是空间碎片激光测距过程。以云南天文台实验系统对于空间碎片no.24809进行激光测距实验为例，其激光测距过程的有关参数为：激光波长λ为532nm，planck常数h为6.63×10-34j·s，光速c为3×108m/s，激光发射能量et为0.8mj，大气透过率t为0.5，发射系统和接收系统光学效率kt，kr分别为0.8，望远镜接收口径d为1.2m，入射激光发散角θt为12”，空间碎片距离r为1000km，空间碎片散射截面积s为12m2，空间碎片等效反射系数ρ为0.4。将上述数据代入下述空间碎片回波光子数方程：
得到该碎片反射光子数s为0.1个，意味着每发射10个脉冲，平均有1个反射回波到达探测器接收面。因为碎片本身为漫反射，其回波展宽根据经验估计为6m，对应回波展宽时长techo为0.04μs。
根据观测经验，天光背景噪声速率约为2mhz，探测器暗噪声为300khz，因此整个测距过程噪声速率vnoi约2.3mhz。设置门控时长tgate为1.3μs，回波在门控中的位置γ为0.5。那么回波分布范围内的噪声数为n＝vnoi·techo为0.09个，回波分布范围之前的噪声光子数为w＝vnoi·(tgate-techo)·γ为1.45个，回波分布范围之后的噪声光子数为q＝vnoi·(tgate-techo)·(1-γ)为1.45个。
将s、n、w、q代入本发明推导的探测概率和虚警概率计算式：
pd＝e-w·(1-e-s-n)(8)
pfalse＝1-e-w+e-w-s-n-e-w-s-n-q(9)
计算可得空间碎片的探测概率为3.86％，虚警概率为91.39％，与附图1中所示探测概率与虚警概率随回波和噪声强度的变化规律相一致。
进而，选取5000个激光脉冲以得到统计的直观反映，那么统一的光子序列中，共有15000个噪声光子，500个目标反射光子，共计15500个光子。将上述按照在门控内的距离先后进行排序，如附图2所示。
按照光子先后顺序，从第1个光子开始判断，第1个光子被探测的概率η＝1-e-1为0.6321，生成0-1范围内满足均匀分布的随机数δ为0.1343，那么第1个光子未被探测，继续判断第2个光子。第2个光子被探测的概率也为0.6321，重新生成0-1范围内满足均匀分布的随机数δ为0.4274，那么第2个光子未被探测，继续判断第3个光子。第3个光子被探测的概率同样为0.6321，冲洗生成0-1范围内满足均匀分布的随机数δ为0.8825，因此第3个光子被成功探测。
第3个光子的距离r3为0.037m，需在此基础上叠加探测器时间抖动随机数作为输出信号。探测器的时间抖动为100ps，时间延迟为0.67μs，换算成距离，距离抖动为0.015m，距离延迟为0.1m，对应的高斯分布均值为0.1，均方差为0.0035。生成符合时间抖动高斯分布的随机数δ为0.107，因此与第3个光子对应的输出信号l1＝r13+δ＝0.037m+0.107m＝0.144m。
探测器死时间为20ns，对应距离为3m。因为第4个-第251个光子的距离与第3个光子的距离差小于3m，因此需要继续从第252个光子开始判断是否被探测。
第252个光子被探测的概率也为0.6321，重新生成0-1范围内满足均匀分布的随机数δ为0.3295，那么第252个光子未被探测，继续判断第253个光子。重复上述光子判断过程，直到第15500个光子判断完成，并对其中成功探测的光子叠加探测器时间抖动作为输出信号。
将所有生成的输出信号{lx}进行绘图，可得到空间碎片激光测距散点分布仿真图，如附图3(a)所示，其中纵坐标为输出信号叠加碎片预估轨道距离后的距离残差值。在背景噪声之中，能够观察到一条光子密度较高的曲线，该曲线就是回波曲线，可供后续数据处理和回波识别。与空间碎片实测结果附图3(b)对比可以发现，本发明提供的仿真分析方法能够很好地仿真空间碎片激光测距效果。
基于上述方法，相应的，本发明还提供了一种空间碎片激光测距效果仿真分析装置，包括探测模块、仿真模块和概率计算模块。
探测模块，用于针对光子时刻序列，按照光子先后顺序，判断光子是否被探测；探测方式为：在0-1的均匀分布范围内，生成随机数δ，判断随机数δ与单个光子探测效率η＝1-e-1的大小关系；如果δ≤η，那么该光子被成功探测，如果δ>η，那么该光子未被探测，则继续判断后一个光子，直到所有光子判断完成。
这里，在按照顺序判断光子是否被探测时，若当前光子判定为被探测，则不再探测探测器死时间内到达的光子，直接探测死时间之后的光子；若当前光子不被探测，则继续判断后一个光子，直到所有光子判断完成。
仿真模块，用于基于光子探测结果，生成仿真的探测器输出信号。具体来说，是在成功探测的光子上叠加探测器时间抖动的影响，生成仿真的输出信号。将仿真的探测器输出信号进行显示，就获得了空间碎片激光测距的散点分布仿真图。
概率计算模块，用于按照公式(8)和(9)计算探测过程的探测概率pd和虚警概率pfalse；
pd＝e-w·(1-e-s-n)(8)
pfalse＝1-e-w+e-w-s-n-e-w-s-n-q(9)
其中，w为回波分布范围之前的噪声光子数，s为经空间碎片反射后的回波光子数，n为回波分布范围内的噪声光子数，q为回波分布范围之后的噪声光子数。进一步地，通过改变经空间碎片反射后的回波光子数s和噪声强度vnoi，可以得到探测概率和虚警概率随二者的变化曲线图。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。例如，本发明限于回波为高斯分布、噪声为平均分布的情况，同样适用于其他的回波、噪声分布形式。