引潮位频谱展开算法研究 许剑伟 (福建莆田第十中学) 雷伟伟 (河南理工大学) 关键词:引潮位展开,固体潮,频谱展开,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
|