莆田第十中学 欢迎您!
今天是:
当前位置: > 名师频道 > 骨干教师 > 教学研究
引潮位频谱展开算法研究
【字体:[大][中][小] 】【发布时间:2017-11-21】 【作者:/来源: 】 【阅读:次】【关闭窗口】

引潮位频谱展开算法研究

许剑伟 (福建莆田第十中学)

雷伟伟 (河南理工大学)

关键词:引潮位展开,固体潮,频谱展开,FFT,最小二乘法

  要:Wenzel制作了HW95引潮位展开表,适用于18502150。现行星历精度已有很大提升,适用年代更宽,有必要推演适用年代更宽、精度更高的谱表。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:平黄经多项式表达系数表。单位:°/儒略千年,力学时

K

常数项

t

t2

t3

t4

1

2

3

4

5

6

7

8

9

10

11

242.14980452999

218.31664562999

280.46645016002

 83.35324311998

234.95544499000

282.93734098001

252.25090551999

181.97980085000

355.43299958002

 34.35151874003

50.07744430000

127037328.88553056

   4812678.81195750

    360007.69748806

    40690.13635250

     19341.36261972

        17.19457667

   1494740.72172233

    585192.12953330

    191416.96370297

     30363.02774848

     12235.11068622

0.17696111

0.14663889

0.03032222

1.03217222

0.20756111

0.04568889

0.03034984

0.03101395

0.03105187

0.02232972

0.0519078

0.00183140

0.00185140

0.00002000

0.01249168

0.00213942

0.00001776

0.00001811

0.00001490

0.00001564

0.00003701

0.00002985

0.00008824

0.00015355

0.00006532

0.00052655

0.00016501

0.00003323

0.00006532

0.00006532

0.00006532

0.00005214

0.00009740

表2:各黄经项含义

K

参数

符号

本文标识

频率标识

(t项系数)

频率(°/小时)

1

2

3

4

5

6

7

8

9

10

11

平太阴时

月亮平黄经

太阳平黄经

月亮近点平黄经

月亮负升交点平黄经

太阳近点平黄经

水星平黄经

金星平黄经

火星平黄经

木星平黄经

土星平黄经

τ

s

h

p

N’

ps

Lmer

Lven

Lmar

Ljup

Lsat

τ

λ2

λ3

λ4

λ5

λ6

λ7

λ8

λ9

λ10

λ11

ω1

ω2

ω3

ω4

ω5

ω6

ω7

ω8

ω9

ω10

ω11

14.49205212018

0.54901651973

0.04106863991

0.00464181341

0.00220640687

0.00000196151

0.17051571090

0.06675703052

0.02183629520

0.00346372664

0.00139574614

 

迭代算法流程图描述

迭代算法描述(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)

HW95

LW600

项目

均方差

最大

最小

项目

均方差

最大

最小

Dgr

Dgx

Dgy

DV

1.516

1.499

1.008

4.739

9.420

 8.492

 4.852

29.725

-9.511

 -8.527

 -6.850

-27.868

Dgr

Dgx

Dgy

DV

0.859

0.711

0.608

2.421

  5.650

  3.586

  3.905

 12.023

-5.037

 -4.411

 -5.451

-15.598

 

 

 

 

 

 

作者:许剑伟,雷伟伟

地址:福建省莆田市秀屿区莆田第十中学,351146

电话:13850262218