初轨计算精度的bootstrap估计

第５４卷第１期　天　文　学　报　Ｖｏｌ＿５４　ＮＯ．１　２０１３年１月　ＡＣＴＡ　ＡＳＴＲＯＮＯＭＩＣＡ　ＳＩＮＩＣＡ　Ｊａｎ．，２０１３　初轨计算精度的ｂｏｏｔｓｔｒａｐ估计　王歆　，２十　（１中国科学院紫金山天文台南京２１０００８）　（２中国科学院空间目标与碎片观测重点实验室南京２１０００８）　摘要　从非参数统计角度，在仅有观测资料而没有其它信息的情况下，分别给出了初　轨计算的精度和置信区间的估计方法．该方法基于ｂｏｏｔｓｔｒａｐ方法，仅依赖于观测资料，　不需要精密定轨结果作为参考，也无需假设观测资料误差呈正态分布，并且对各种初轨计　算方法均适用．数值计算结果表明该方法应用简便，可作为初轨的评估和后续应用的重要　参考．　关键词航天器，天体力学，方法：统计　中图分类号：Ｐ　１３５；　文献标识码：　Ａ　１　引言　初轨计算，即二体意义下的短弧定轨，是天体力学中的经典问题，已有２００多年的　历史，至今仍是一项不可缺少的工作，特别是光学测角资料的初轨计算．初轨计算的主　要目的是通过采集到短弧数据计算出空间目标的轨道，从而给出空间目标的预报，以求　获取后续的观测数据．由于后续能否捕获目标取决于初轨计算的精度，人们一直在探索　提高初轨计算的精度的方法，但由于初轨计算方程的病态是本质的，各种方法只能改善　这种影响而不能彻底解决．　本文从另一个角度看待初轨计算，既然初轨计算的精度提高十分困难，那么不再追　求提高初轨的精度，而是在得到初轨后给出计算结果的不确定性的度量．特别是现在随　着单台设备的视场越来越大，观测网络中观测设备越来越多，使得空间目标的捕获不再　完全依赖于预报．只要能够给出空间目标方位范围的预报，对于当前设备能力捕获目标　并不困难．特别是近年我国有多套大天区覆盖设备投入使用，对于此类设备，空间目标的　观测已彻底无需预报．初轨计算也逐步拓展为服务于观测数据的关联和轨道识别等，此　时不确定性的度量尤为重要，区间估计比点估计在相关应用中更有优势和理论说服力．　文献［１－２】探讨了初轨计算的精度，引入了Ｃｒａｍ￣ｒ—Ｒａｏ下界（ＣＲ下界）作为初轨计　算的最佳精度．ＣＲ下界是最小方差无偏（ＭＶＵ）估计的方差的下界，ＣＲ下界并不直接　给出求得ＭＶＵ估计的方法，而是在得到一个可能的ＭＶＵ估计后，看其方差是否达到了　ＣＲ下界．由于初轨计算是非线性的，因此无法获得其精度的分析表达形式，在已知真值　２０１２—０６—２５收到原稿，２０１２—０７－２５收到修改稿　ｔ　ｗａｎｇｘｉｎ＠ｐｍｏ．ａｃ．ｃｎ　７４　天　文　学　报　５４卷　和误差分布情况下，可通过大量模拟计算方法求得初轨精度，文献『１－２］采用的也是这种　方式．如果已知初轨计算方法是ＭＶＵ估计，那么此时估计的精度接近ＣＲ下界，可以　通过计算ＣＲ下界来表示初轨计算的精度．　ＣＲ下界中涉及Ｆｉｓｈｅｒ信息量的计算，仍需　要已知观测误差的分布（一般假设为正态分布）以及无误差的轨道根数，实践中需要采用　精密定轨的结果．　可见，上述工作中初轨精度的获得都需要有额外信息，包括真值和误差分布（甚至要　求是正态分布）．而在实践中需要进行初轨计算时，往往是在没有其它信息情况下，本文　通过引入ｂｏｏｔｓｔｒａｐ方法，仅从采集到的观测资料本身出发，给出了一种初轨计算精度的　估计方法，并在此基础上给出了初轨的区间估计．　２　Ｂｏｏｔｓｔｒａｐ方法　令随机变量ｙ服从分布Ｆ，统计量　＝　（Ｆ）．由Ｍｏｎｔｅ　Ｃａｒｌｏ思想，可以通过模拟　的方法得到　的估计：　去　）一厂　（Ｆ）ｄＦ（ｙ）　遗憾的是Ｆ未知时无法进行模拟采样．　获取了一组观测量记为ｙｌ，　，…，　，　为观测数据的经验分布函数［３】．　（ｙ）＝　１∑　），　其中　（　）为指标函数，当　Ｙ时取１，其它情况取０．由插入原理（ｐｌｕｇ—ｉｎ　ｐｒｉｎｃｉｐｌｅ）［引，　在没有其它关于Ｆ的信息时，　可由　＝　（　）估计．这样只要从分布　中模拟数据，　即可用Ｍｏｎｔｅ　Ｃａｒｌｏ方法求得估计量，而岛在获取样本后便是已知的．　由　的定义，每个　都有相同的Ｉ／ｎ的概率，从　随机选取ｎ个点和可放回地　从原数据样本抽取样本量为凡是等同的．因此可形成模拟过程：　（１）有放回地从　，ｙ２　．，　中取出ｎ个数据的样本，构成ｂｏｏｔｓｔｒａｐ样本　，　．．．．，　（２）计算统计量　＝　（　，　…．，　）；　（３）重复Ｂ次步骤（１）和（２），获得磅，ｉ＝１，２，…，Ｂ．　上述过程为ｂｏｏｔｓｔｒａｐ方法的基本过程，　ｂｏｏｔｓｔｒａｐ方法是由Ｅｆｒｏｎ于１９７９年提出　的［４】＿本文用上标　表示估计量对应于重采样的ｂｏｏｔｓｔｒａｐ样本．　那么　的方差的ｂｏｏｔｓｔｒａｐ估计为　／　＼２　ＶＡＲ＝／（　（ｙ））。ｄＦ（Ｙ）一（／Ｔ（Ｙ）ｄＦ（Ｙ））　Ｐ　＝／　、２　／（　（ｙ））。ｄＦｎ（Ｙ）一（／　（ｙ）ｄ　（ｙ））　Ｂ　／　Ｂ　、２　　＝Ｂ∑（　）。一（　１∑醯）　１期　王　歆：初轨计算精度的ｂｏｏｔｓｔｒａｐ估计　７５　标准误差ｓＥ＝　　．３初轨计算精度的估计　初轨计算中参数一般为历元时刻ｔｏ的位置速度矢量或者轨道根数０，与观测值的关　系为　＝　Ｂ∑　一日　Ｙｉ＝ｈ（ｅ，ｔｉ）＋　ｉ＝１，２，…，Ⅳ，　（４）　／，●●其中Ⅳ为观测资料个数，ｅｔ为观测误差，仅需要假定ｅ　独立同分布　一一／　初轨计算即在获　取了一组观测量｛Ｙ１Ｂ∑　一　口　，ｙ２…．，ＹＮ｝后估计ｐ：　＝ｑ（Ｙ１，Ｙ２，¨．，ＹＮ）．　（５）　采用ｂｏｏｔｓ＼、／　ｔｒａｐ方法估计　的方差，需要对独立同分布的随机量进行重采样，显然初轨　计算中这个量为观测误差，观测误差可用残差表示：　Ｙｉ＝ｈ（ｅ，ｔｉ），　ｉ＝１，２，…，Ⅳ　（６）　乏ｔ＝Ｙｉ—Ｙｉ，　ｉ＝１，２，…，Ⅳ，　（７）　通过对乏ｉ的重采样，可以构建初轨计算的ｂｏｏｔｓｔｒａｐ方法，步骤如下：　（１）由观测数据｛Ｙ１，Ｙ２，．．．，ＹＮ｝计算初轨，得到　；　（２）根据９，求得｛Ｙ１，ｙ２，．．．，ｙⅣ）和｛乏１，舍２，．．．，　Ⅳ）；　（３）在［１，Ⅳ】区间内生成Ⅳ个均匀分布的随机整数／＇ｉ，ｉ＝１，２，…，Ⅳ，构成ｂｏｏｔｓｔｒａｐ　样本｛乏　舍　…，Ｅｒ　Ｎ）；　（４）生成伪观测量｛ｙ　，ｙ；，．．．，ｙ　）＝｛ｙ１＋仑…Ｙ２＋乏ｒ２ｊ…，ＹＮ十　）；　（５）由伪观测量计算初轨　＝ｑ（Ｙ７，Ｙｉ，．．．，ｙ★Ⅳ）；　（６）重复Ｂ次步骤（３）一（５），获得　，ｉ＝１，２，…，Ｂ．　根据（３）式可得到初轨计算的方差估计，即精度估计．　４初轨计算的区间估计　定义区间Ｃ＝　，　］和枢轴变量Ｒ＝　一　，根据区间估计定义，有：　Ｐ（ｃ　）　＝　ＪＦ）（　一　一　—　）　＝Ｐ（Ｏ—ｃｕ　一　—　）　＝Ｐ（Ｏ—　Ｒ　一ｃＩ）　＝Ｐ（Ｒ　一ｃＩ）一Ｐ（Ｒ　—　），　７６　天　文　学　报　５４卷　显然当　一　和　Ｃｉ分别为Ｒ的　／２和ｌ一　／２分位数ｒ。／２和ｒｌ一　／２时，有：　Ｐ（　６）　）＝１一　一　＝１－－Ｃｔ￣　（９）　即　Ｐ（ｏ∈Ｃ）＝１一ＯＺ，　此时Ｃ一［Ｃｔ，　］为０的１一　置信区间．　（１０）　由ｆ　Ｒ的分布未知，但可用　一　来估计［５】．仍采用ｂｏｏｔｓｔｒａｐ方法估计，可得：　＝　一　，　ｉ一１，２，…，Ｂ．　（１１）　可见　分位数与　的分位数是对应的，可得１一ＯＺ　ｂｏｏｔｓｔｒａｐ置信区间　＝［　，　］　＝　，　／２＝２　／２，　（１２）　（１３）　Ｇ　…０　ｒ＊￣／２　２０—０＊．／２．　５数值验证　由以上讨论可见，采用的初轨计算方法并不影响上述过程，本文数值验证采用改进　的Ｌａｐｌａｃｅ型基本方法，测量方程为　ｒ　Ｒ＝ｐＬ　ｆｌ４１　其中　为目标位置矢量，Ｒ为测站位置矢量，ｐ为目标的斜距，　方程两端叉乘　，并由ｒ＝ｆｒｏ＋－ｑ　０，可得到条件方程：　ｆｒｏ×Ｌ＋ｇ￣０×Ｌ＝Ｒ×Ｌ　为目标的方向矢量，　ｆ１５１　３个分量方程任选两个独立的方程．条件方程组由最小＿二乘法求解得到（ｒｏ，　ｏ），重新计　算，和ｇ后迭代，收敛后即获得初轨计算结果．　选取中国科学院空间目标与碎片光学观测网中昆明站的一圈实测资料进行了数值验　证，资料共采集２６３个点，弧长６　ｍｉｎ，测角精度优于５　．　表１给出了该圈资料的计算结果，其中０表示原始观测数据初轨计算结果；０为采　用多站多天精密定轨（微分改进）得到的结果，此处作为真值用作比对；　的精度用标准　误差ＳＥ　表示，ｂｏｏｔｓｔｒａｐ重采样１　０００次．计算中１　０００组ｂｏｏｔｓｔｒａｐ采样计算结果均　收敛，并且ＲＭＳ均在３　左右，结果说明ｂｏｏｔｓｔｒａｐ方法可以有效地估计初轨的精度，　计算结果符合初轨计算的规律，轨道面相关的两个根数ｉ和Ｑ精度较好，而轨道半长径　ａ的精度较差．　初轨计算中最难确定的量是轨道半长径。，往往用其来代表初轨计算的精度．表２针　对不同弧长情况给ｆｊ｛了ａ的计算结果．各组历元均归算到和表１相同时刻．表中“为原　１期　王　歆：初轨计算精度的ｂｏｏｔｓｔｒａｐ估计　始观’？贝４资料的初轨计算结果，ａ的精度和９５％置信区间由ｂｏｏｔｓｔｒａｐ方法给出．每组计　算ｂｏｏｔｓｔｒａｐ重采样１　０００次，所有５　０００组结果均收敛，并且ＲＭＳ都很小，如第１组　３５　Ｓ弧长的所有１　０００组计算结果虽然ａ差异较大，但对应相同弧段的原始资料的ＲＭＳ　仅为１．５４　２．１７　．众所周知，观测弧长越长，初轨的精度越好，表中的定量结果与此吻　合．精密定轨的结果都落在了置信区间内，说明区间估计对于初轨在后续应用或在此基　础上的推断是有效的．特别是对于空间目标搜索，在获得极短弧段后，初轨计算精度一　般不能满足后续直接跟踪捕获的要求，采用区间估计后可提供后续搜索范围．　表１初轨计算的精度　Ｔａｂｌｅ　１　Ｔｈｅ　ａｃｃｕｒａｃｙ　ｏｆ　ｐｒｅｌｉｍｉｎａｒｙ　ｏｒｂｉｔ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　表２　不同弧长的初轨计算精度　Ｔａｂｌｅ　２　Ｔｈｅ　ａｃｃｕｒａｃｙ　ｏｆ　ｐｒｅｌｉｍｉｎａｒｙ　ｏｒｂｉｔ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　ｗｉｔｈ　ｄｉｆｆｅｒｅｎｔ　ａｒｃ　ｌｅｎｇｔｈｓ　表３模拟了不同分布的观测误差，计算了相应的轨道半长径ａ的精度，给出了ｂｏｏｔ—　ｓｔｒａｐ估计ＳＥ　和采用已知分布模拟１　０００组求得的结果ＳＥ．模拟包括了正态分布和均　匀分布，分布均值为０，标准差　分别选取为２　、５”、１０　以及６０　．表中数据表明不　利用已知分布信息的ｂｏｏｔｓｔｒａｐ估计结果和利用已知分布计算的结果相当，　ａ的精度与　成正比，相同方差不同分布下的结果也不同．Ｂｏｏｔｓｔｒａｐ方法由于不依赖于观测误差的　特定分布，对各种误差情况也都有较好的响应．　表３　不同观测误差分布下的初轨计算精度　Ｔａｂｌｅ　３　Ｔｈｅ　ａｃｃｕｒａｃｙ　ｏｆ　ｐｒｅｌｉｍｉｎａｒｙ　ｏｒｂｉｔ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　ｗｉｔｈ　ｄｉｆｆｅｒｅｎｔ　ｄｉｓｔｒｉｂｕｔｉｏｎｓ　ｏｆ　ｎｏｉｓｅ　７８　天　文　学　报　５４卷　６结论与讨论　ｌ　２　３　４　５　综上所述，ｂｏｏｔｓｔｒａｐ方法可以“自助地”给出初轨计算精度，只需要观测数据，因　此在没有其它信息时可实现对初轨计算精度的估计．　Ｂｏｏｔｓｔｒａｐ方法理论比较复杂，但　应用却十分简便，对于各种初轨计算方法都适用，并且仅需要观测资料误差为独立同分　布，并不要求为正态分布．　Ｂｏｏｔｓｔｒａｐ方法应用中需要注意的是样本分布能够较好地描述总体分布，这需要有两　方面的保证．一方面需要有较多的观测资料，得益　探测手段的发展，目前空间目标观　测中一般能够接近１　Ｈｚ的采样率，很短的弧段也能够获取较多的观测资料；另一方面观　测数据中没有野值，当存在野值时，重采样的样本会受到污染，影响估计精度．因此需要　对剔除野值的样本进行重采样，而这并无实质的困难．　参考文献　贾沛璋，吴连大．天文学报，　１９９８，３９：３３７　Ｊｉａ　Ｐ　Ｚ，Ｗｕ　Ｌ　Ｄ．ＣｈＡ＆Ａ，１９９９，２３：３８４　Ｅｆｒｏｎ　Ｂ，Ｔｉｂｓｈｉｒａｎｉ　Ｒ．Ａｎ　Ｉｎｔｒｏｄｕｃｔｉｏｎ　ｔｏ　ｔｈｅ　Ｂｏｏｔｓｔｒａｐ．Ｎｅｗ　Ｙｏｒｋ：Ｃｈａｐｍａｎ＆Ｈａｌｌ，１９９３：３１—３８　Ｅｆｒｏｎ　Ｂ．ＡｎＳｔａ．１９７９　７：１　Ｓｉｎｇｈ　Ｋ　ＡｎＳｔａ，１９８１，９：１１８７　Ｔｈｅ　Ｂｏｏｔｓｔｒａｐ　Ｅｓｔｉｍａｔｉｏｎ　ｏｆ　ｔｈｅ　Ａｃｃｕｒａｃｙ　ｏｆ　Ｐｒｅｌｉｍｉｎａｒｙ　Ｏｒｂｉｔ　Ｄｅｔｅｒｍｉｎａｔｉｏｎ　ＷＡＮＧ　ＸｉｎＩ，２　（１　Ｐｕｒｐｌｅ　Ｍｏｕｎｔａｉｎ　Ｏｂｓｅｒｖａｔｏｒｙ，Ｃｈｉｎｅｓｅ　Ａｃａｄｅｍｙ　ｏｆ　Ｓｃｉｅｎｃｅｓ　Ｎａｎｊｉｎｇ　２１　０００８）　（２　Ｋｅｙ　Ｌａｂｏｒａｔｏｒｙ　ｏｙ　ａｃｅ　Ｏｂｊｅｃｔ　ａｎｄ　Ｄｅｂｒｉｓ　Ｏｂｓｅｒｖａｔｉｏｎ，Ｐｕｒｐｌｅ　Ｍｏｕｎｔａｉｎ　Ｏｂｓｅｒｖａｔｏｒｙ，Ｃｈｉｎｅｓｅ　Ａｃａｄｅｍｙ　ｏｆ　Ｓｃｉｅｎｃｅｓ，Ｎａｎｊｉｎｇ　２１　０００８）　ＡＢＳＴＲＡＣＴ　Ｆｒｏｎｌ　ｔｈｅ　ｐｏｉｎｔ　ｏｆ　ｖｉｅｗ　ｏｆ　ｔｈｅ　ｎｏｎ—ｐａｒａｍｅｔｒｉｃ　ｓｔａｔｉｓｔｉｃｓ．ａ　ｇｅｎｅｒａｌ　ｅｓｔｉｍａｔｉｏｎ　ｍｅｔｈｏｄ　ｏｆ　ａｃｃｕｒａｃｙ　ａｎｄ　ｃｏｎｆｉｄｅｎｃｅ　ｉｎｔｅｒｖａｌ　ｏｆ　ｔｈｅ　ｐｒｅｌｉｍｉｎａｒｙ　ｏｒｂｉｔ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　ｉＳ　ｐｒｏ—　ｐｏｓｅｄ　ｆｏｒ　ｔｈｅ　ｏｃｃａｓｉｏｎ　ｗｉｔｈｏｕｔ　ａｎｙ　ｏｔｈｅｒ　ｉｎｆｏｒｍａｔｉｏｎ．Ｔｈｅ　ｍｅｔｈｏｄ．ｂａｓｅｄ　ｏｎ　ｔｈｅ　ｂｏｏｔｓｔｒａｐ　ｍｅｔｈｏｄ．ｅｍｐｌｏｙｓ　ｏｎｌｙ　ｏｂｓｅｒｖａｔｉｏｎ　ｄａｔａ　ｗｉｔｈ　ｍｉｎｉｍａｌ　ａｓｓｕｍｐｔｉｏｎｓ　ａｂｏｕｔ　ｔｈｅ　ｄｉｓｔｒｉｂｕｔｉｏｎｓ　ｏｆ　ｏｂｓｅｒｖａｔｉｏｎａｌ　ｎｏｉｓｅ　ａｎｄ　ｄｏｅｓ　ｎｏｔ　ｒｅｑｕｉｒｅ　ｔｈｅ　ｒｅｓｕｌｔ　ｏｆ　ｐｒｅｃｉｓｅ　ｏｒｂｉｔ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　ａｓ　ａ　ｒｅｆｅｒｅｎｃｅ．Ｎｕｍｅｒｉｃａｌ　ｃｈｅｃｋ　ｓｈｏｗｓ　ｔｈａｔ　ｔｈｅ　ｍｅｔｈｏｄ　ｉＳ　ｖｅｒｙ　ｅａｓｙ　ｉｎ　ｉｍｐｌｅｍｅｎｔａｔｉｏｎ．ａｎｄ　ｔｈｅ　ｅｓｔｉｍａｔｉｏｎ　ｒｅｓｕｌｔ　ｉＳ　ｂｅｎｅｆｉｃｉａ１　ｆｏｒ　ｂｏｔｈ　ｔｈｅ　ｅｖａｌｕａｔｉｏｎ　ｏｆ　ａｃｃｕｒａｃｙ　ａｎｄ　ｆｏｌｌｏｗ—ｕｐ　ｉｎｆｅｒｅｎｃｅ．　Ｋｅｙ　ｗｏｒｄｓ　ｓｐａｃｅ　ｖｅｈｉｃｌｅｓ，ｃｅｌｅｓｔｉａｌ　ｍｅｃｈａｎｉｃｓ，ｕｌｅｔｈｏｄｓ：ｓｔａｔｉｓｔｉｃａｌ