|
引潮位频谱展开算法研究 许剑伟 (福建莆田第十中学) 雷伟伟 (河南理工大学) 关键词：引潮位展开，固体潮，频谱展开，FFT，最小二乘法 摘 要：Wenzel制作了HW95引潮位展开表，适用于1850到2150。现行星历精度已有很大提升，适用年代更宽，有必要推演适用年代更宽、精度更高的谱表。Wenzel原文对展开原理只做了简要叙述，难以重现算法，为此本文重新设计了一款基于“FFT+最小二乘法”的频谱展开方案，获得成功。当采用600年星历，只须展开17000项就可获得0.86pm/s2的精度(均方差)，最大误差为5.5pm/s2 一、引力位球谐函数的频谱概念 1、基本公式表达 引力位球谐函数表示为：
式中r是测站与地心距离，是地心纬度，a是地球半径，GMb是天体引力常数，δ是天体赤经，R是天体的地心距，H是天体的时角，式中是完全规格化缔合勒让德函数。 是大地系数，与地理坐标相关，与时间无关，不必做频谱展开。 公式的后面一串记为：，它是时间相关的量，下文称“时域公式”。所谓“引潮位频谱展开”就是把“时域公式”展开为如下形式：
或者
二、谱线提取的思想方法 利用“时域公式”计算出时域序列，FFT变换后即得频谱。由于谱线之间的互相干扰，没能一次性将所有谱线精确提取出来。如果找到最大谱线的频率，然后利用最小二乘法拟合该谱线的准确的振幅及相位，并在时域序列中扣除该谱线，得时域残余值，残余值序列重新按上述方法提取谱线，直到残余值足够小为止。这样就可以得到所须的谱线并列表。 进行最小二乘法操作时，要选取适当的正弦拟合函数。不同的拟合函数，效果也有所不同。本文采用
这是四参数拟合法，在拟合范围内可以收敛到很高的精度，拟合范围之外，精度较差。注意，t项系数参与拟合是有必要的，否则难以收敛。HW95采用单t系数拟合，本文采用C1与S1两个t系数拟合。 三、时域计算 进行频谱展开前，要事先计算时域序列。而且计算前，还需要天体坐标。我们采用DE406计算天体坐标，然后进行坐标变换，转换到当日赤道坐标系，再转换到时角坐标系。方法是： 用DE406算出天体地心直角坐标 => 转到J2000黄道直角坐标 => 直角坐标转球面坐标 => 转到当日坐标系(岁差改正) => 黄经章动改正 => 转到赤道坐标(含交角章动改正) => 转到时角坐标(dT=0)。 dT的影响与地理经度L的影响是类似的，生成谱表之后，计算时引潮位时只须在角度项中补回dT引起的角偏移即可。 四、频谱展开算法描述 下表是计算过程所须的一些参数，沿用HW95使用的参数，Simon et al.(1994) 对这些参数没有特别的精度要求，谱表的最终精度是由迭代的截断阈值以及时域基准序列精度决定的。 表1：平黄经多项式表达系数表。单位：°/儒略千年，力学时
表2：各黄经项含义
迭代算法流程图描述
迭代算法描述(1个星体的某阶次迭代算法) 1、在时域中计算出某阶次潮波的序列值，记为<Y>。 2、序列加汉宁窗，得<y>序列（只对月亮加窗，行星不加窗），并把<y>导入FFT 3、利用FFT把时域序列转换为频域序列。 4、在频域序列中找到最大振幅谱线所对应的频率ω并做频率插值。 利用Voglewede方法做频率插值，得到最大谱线的准确频率ω 针对月亮，如果是20至39次迭代，频率直接使用0到19次迭代用的频率。 5、做角度项最佳整系数拟合，重新确定最大谱线的频率ω 一般把频谱的角度项表达为各项黄经的组合。
第m次(指阶次中的m)的频谱展开，可以认为频谱是以mτ的频率(即mω1)为中心展开的。所以整数组合中习惯性规定必须含有mτ，相当于N1=m 即： 组合须满足频率关系： 拟合误差尽量做到小于，是FFT的频谱分辨力。 用最小二乘法确定系数时，ω采用整数拟合后的ω 如果是m=0阶次，当遇到第1行迭代，ω应置0并取消加窗操作。 6、确定谱线的振幅 把谱线表达为： <y>与h(t)进行最小二乘拟合得到该谱线的振幅参数 针对月亮2阶，如果是前20次迭代，采用拟合，把t项的系数置为0 。月亮2阶前20次迭代得到的谱线振幅很大，谱线之间干扰比较严重，对t项系数的值影响较大，造成拟合区之外的精度下降过多，所以前20项只拟合常数项C0和S0，等20条大谱线吸收之后，再重算这些谱线的t项系数。 7、在<Y>序列中减去谱线h(t)，得到新的<Y>序列 8、按指定格式输出谱线到内存缓存区。 输出参数为：N1、N2……N11、参数，同时将输出到屏幕 输出时，检查该频率谱线是否已在缓存表中出现过，如果已存在，则累加合并即可。 9、如果h(t)的振幅小于阈值b则结束迭代 记当前迭代谱线振幅，当时结束迭代 迭代过程中，a跳动较大，如果意外跳到接近于0就会造成提前结束迭代。所以迭代过程中对a进行加权累加，即，当时结束迭代。 10、重复以上2到9步骤 11、把缓存数据输出到谱表文件 输出振幅小于阈值b的项，振幅，只输出的项。 12、谱表验算 用缓存谱表计算引潮位与时域计算结果多点比对，算出截断误差(方差的开方)。 五、关于谱展开的一些细节 关于加窗 月亮引潮位最大，对计算的相对精度要求最高，须做加窗处理。 FFT过程与最小二乘过程，二者对谱线的幅度提取过程有相似之处的。从FFT的角度看，矩形窗采样非同步的信号，频谱泄漏严重，影响了振幅的准确提取。所以，有必要对月亮进行加窗。加窗之后，收敛速度明显加快。 不过，经过几百次加窗迭代计算之后，信噪比很低，继续加窗迭代反而会降低收敛速度。因为，加窗之后，频率分辨反而下降，当信噪比很低时，频率插值精度过低，插值得到的频率并不对应最大振幅项，最小二乘法的吸收谱线的效率下降，收敛速度下降。因此只须前300次迭代加窗。 频率是未知的，一般使用汉宁窗即可。 直流或超低项不能加窗，在m=0的阶次中要特别注意。 行星计算不必加窗。行星计算的谱线只需几条，而且本身已收敛很快，只需前几次迭代加窗，加窗不适当会适得其反，所以没有加窗的必要。 整数拟合： 采用穷举法、贪心算法。但要注意搜索的优先顺序及搜索范围，以得简洁组合。 采样间距与采样周期： 不同“次”的采样间距不同，m=0对应低频，采样间距控制在20小时左右即可。而1次的应使用5小时或更小间距。2次与3次的，使用2.5小时间距，更高次的使用1.25小时。 截断误差： 当迭代得到的振幅小于阈值b后停止迭代，这就产生了截断误差。设当前已输出了N项，那么截断误差可以估计为经验公式 当项数较少时，收敛速度快，实际误差会比公式计算结果小一倍。 总截断误差为各阶次截断误差的平方和的开方。 六、精度验算和比对： 截断阈值采用2*10-8，时域采样1700年到2300年。 以基于DE406星表计算的时域引潮位为基准测算误差。我们分别对HW95展表表以及我们的展开表进行验算，验算结果如下： 起日日数 = 2433295.000(日) 终止日数 = 2469795.000(日) 时间间距 = 87.600(小时)，测试点数 = 10000
引潮位误差DV (单位μm2/s2)，重力误差Dgr Dgx Dgy(单位pm/s2)
作者：许剑伟，雷伟伟 地址：福建省莆田市秀屿区莆田第十中学，351146 电话：13850262218
|
||||||||||||||||||||||||||||||||||||||||||||||||
莆田第十中学 欢迎您！
今天是:
秀屿公安备案35030502000133号