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日