DE406的切比雪夫多项式的压缩
把某一星体的某一坐标记作y = f(x),x表示时间量,在DE406中,借助切比雪夫多项式分段插值计算f(x)。其中水星,每段16天,所需系数有14个。为了方便计算,定义x的1单位对应水星16天,金星64天,地球32天……。下文以水星为例。
因此,计算 n到n+1之间的某一个x的值,要使用到DE406的第n段的14个水星系数。要计算n到n+m之间(m可取为某一整数)的某一个x值,则须准备m段系数,即14*m个系数。可以考虑把m段合并为1个新段,这个段所须系数至多为14*m,比如m取4,那么重新拟合的系数至多56个就可表达出原数据的所有高频细节。“合并算法”实际上就是在宽度为m的范围内利用傅里叶级数重新计算切比雪夫系数。事实上,DE406系数有很大沉余,只需35个系数就可以精确到1毫角秒,因此样可减小38%的系数。通过适当线性变换,每个系数可使用一个长整型表示,数据量比使用双精度减小一半。各阶系数随着阶数增加而减小,如下表所示,因此按实际所需的字节数保存数据还可减少40%左右。最后压缩率为0.62*0.5*0.6=0.19,本来水星数据量为44MB,压缩后约为8.5MB。对于木星、土星等长周期行星,m的值可取值为32或64甚至为128,压缩倍数可达几十倍到几百倍,压缩后的数据量往往只有100KB数量级。
系数号 最大系数 字节
0 668028793 4
1 1952107708 4
2 1454924545 4
3 749466213 4
4 293059421 4
5 134509767 4
6 55384297 4
7 27599110 4
8 12493479 4
9 6497368 3
10 3097608 3
11 1647385 3
12 814013 3
13 439547 3
14 222497 3
15 121500 3
16 62636 3
17 34473 3
18 18039 2
19 9985 2
20 5290 2
21 2944 2
22 1573 2
23 880 2
24 473 2
25 266 2
26 144 2
27 81 1
28 44 1
29 25 1
30 14 1
31 8 1
32 4 1
33 2 1
34 1 1
当增加系数(即增加阶数)到一定程度后,无法使用精度继续提高,原因有二:(1)已达到原系数的最高精度;(2)已达到截断误差的水平。一般是第一种原因。另外应注意,进行傅里叶级数计算切比雪夫系数时,建议使用n=2的辛普森公式进行数值积分求解傅里叶系数,数值积分所数的f(x)值采样密度最好为DE406每段系数个数的4倍或以上。比如,水星每16到14个系数,那么积分时,16天内采样56点为宜,如果取14点肯定不够,取28点马马乎乎。值得一提的是,不是按时间(即x)等间距采样,而是按角度量θ等间距采样,其中x=cosθ,如果按x等间距采样将造成数值积分上的困难。
切比雪夫的第k阶系数为:Ck=(2/π)∫f(cos(θ))*cos(kθ)dθ,积分范围是0到π,其中第0阶取C0的一半。系数具体积分算法如下,以第k阶系数为例:
h = π/n, H=2/n,(即0到π,对应1到-1,被分为n段,从头到尾共采样n+1点)
θi = h*i, i=0,1,2,……,n
Yi = f(cos(θi)* cos(kθi)
Ck = (Y0+4*Y1+Y2 + Y2+4*Y3+Y4 +Y4+4*Y5+Y6 + ……)/3*H,辛普森数值积分。
如果DE406中的11个坐标量压缩时,m分别取值为4,4,4,8,32,32,32,64,64,4,2,平均精度控制在0.5毫角秒以内,那么压缩后的文件大小由原来190MB变为33MB。以上压缩算法还可很大的改进空间,进一步分析还可以再压缩1半左右,但算法变得比较复杂。适当增加m的取值,可压缩到29MB。如果按比特位压缩(而不是像上面所述的所有坐标按字节压缩),平均精度降到0.5至1毫角秒,那么最后的数据量为23MB以内。如果把坐标变换到黄道坐系,结合比特压缩,可以大大减小z轴坐标数据。太阳的太阳系质心坐标可以考虑使用解析法计算,可减少1.5M的数据量,最后得到的的数据量在20M以内。
建议压缩到33M,即只做3步压缩(1)重建切比雪夫系数减小系数量(2)浮点转整型(3)整型转字节。
阿剑 于家里 设计于2009年2月19日——21日