设计基于VLBI观测数据对日长变化和月地距

BACHELOR'SDEGREETHESISOFWUHANUNIVERSITYOntheCorrelationbetweenChangesinLODandEarth-MoonDistancewithVLBI ：GeodesyandGeomatics ：GeomaticsEngineering ：ChangLiuDirectedby：Prof.ErhuJune学术示例对本所涉及的研究工作做出贡献的其他个人和集体，均已在文中以本人签名 日期 日长变化的研究对气候预测、防治、空间导航、深空探测、大地测量等域的研究有着重要的作用。日长变化与系天体特别是月球与地球的相互作用密切相关，为提高日长变化预报精度，分析这些相关性显得尤为重要。VBI以其极高的相对精度和角分辨率在日长变化等领域发挥了无可替代的作用。基于此，本了一下研究工作：OCCAM6.209-10VLBI数据得到日长变化信息，并简JPL星历计算月地距时间序列；选择合适的插值方法对非等距的日长：VIBI；日长变化；JPL；月地距；互相TheysisoftheLODisofprofoundsignificanceonclimateforecasting,disasterprevention,deepspaceexploration,geodesyandotherfields.Withextremelyhighprecisionandangularresolution,VLBIystheirreceableroleontheysisofthechangesofthelengthofday.HereisacloserelationshipbetweenΔLODandtheinctionofsolarsystemobjectsespeciallytheMoonandtheearth.So,ysisofthesecorrelationsisparticularlyimportantinenhancingtheaccuracyofLODprediction.ΔLODdataisgotthroughthehandingofthe09-10VLBIdatabyOCCAM6.2,andhowtochoosemodelsandsoftwareischieflyyzed;thetimeseriesofdistancebetweentheMoonandEarthbyJetPropulsionLaboratory(JPL)ephemerisarecalculated.Selectedappropriationinterpolationmethodonnon-equidistanttimeseriesofΔLODinterpolation，andprovedtheeffectivenessoftheselectedinterpolationmethod.Trendremovalofthetwotimeseries,andgotperiodsthroughspectralysisandsinecurvefitting.UsedwaveletysistoextracttheperiodiccomponentsofΔLODanddonecross-correlationysisofthem.TheresultsshowthatmonthperiodiccomponentofΔLODaresignificantlynegativecorrelatedwiththedistanceoftheMoonandEarth,semi_monthperiodiccomponentofΔLODaresignificantlypositivecorrelatedwiththeabsoluterateofchangesofdistancebetweentheMoonand:VIBI；ChangesinLOD；JPL；Earth-MoonDistance；Cross-1章绪研究背景和研究现 研究目的和研究意 本文主要内 2章VLBI测量及数据处理原VLBI测量基本几何原 VLBI涉及的坐标系统转 坐标转换流 转换矩阵说 VLBI观测量估计日长变化原 天文和大地参数解 日长变化计 OCCAM对VLBI观测数据的处 OCCAM简 OCCAM解算参数和模型选 3章月地距的计JPL星历简 JPL星历使 基本文件获取及简 JPL星历的使用方 JPL星历计算结 4章数据分析方法、原理及结原始处理结果分 插值方法的选择与原始数据插 插值方法的选 原始数据插 趋势项去 周期项分 频谱分 曲线拟 相关性分 日长变化周期分量提 互相关分 5章结论与结 展 参考文 致 研究背景和研究现

1VLBI是20世纪60年代末期发展起来一种用于射电天文学的无线电技术。全球有40多个地面VLBI站，VLBI的发展有SVLBI（空间甚长基线测量）和e-VLBI（实时甚长基线测量）。国际上正在进行的有关VLBI的项目主要有：地测量观测系统）、SGP（空间大地测量计划）；已经不再进行的项目有CDP(地壳动力学项目)、DOSE（NASA的固体地球研究）；其他的项目还有宇宙航空（ISA）组织的VSO（VLBI空间观测项目）项目、由NSA出资JPI实行的SLBI项目（VO）。VBI以得到更高的数据处理结果。目前VBI数据处理的方法和软件主要有OCCM、Calc/SLK、MOEST、GORIA、、GOAT、SKD、FielSystem、S、SteelBreze、VieV、c5++、AIADNA等。在IVS296个中OCCACalc/OVK软件[1]。日长变化周期的时间尺度从日、年、十年到千年甚至亿年。日长变化方面的研究目前主要集中于日长变化周期的研究及其产生机理的研究。Sphenson和Morrison（1980）系统的计算了日长变化的周期[2]。产生机理的研究这方面的研究主要有冰期和末次冰消期循环、地幔对流和大陆漂流对日长变化千年和亿年尺度周期的影响；地球核幔电磁耦合，地磁场对日长变化的10量对日长的准两年变化，季节性变化，月周期和半月周期变化的激发；还有少量的海洋、光压、潮汐摩擦、水储量、月球等对日长变化影响的研究，(Lmbeck,K.1980)系统总结了日长变化的机理[3]。在日长变化对气候、环境等的影响的研究方面也有一些进展，主要是日长变化对大气环流、、厄尔尼诺、拉尼娜、南方涛动指数的研究，（DEMICHEIS,TZZI,和MELNI,2005）研究了日长变化与温度的相关性[4]。在日长变化的预报这一方面的研究相对较少,日长变化预报方法目前来说有两类,一类是线性模型，其代表是基于自回归的预报模型，一类是非线性模型，其代表是神经网络预报模型(,2007)[5。找出影响日长变化的因素，对日长变化预报来说尤为重要。研究目的和研究意合成为等效直径为最长基线长度的综合孔径望远镜。由于VLBI采用高稳定度原子钟作为独立本振系统，其基线长度理论上可以扩充到所需的任意距离，因而VLBI技术能获得极高的角分辨率。VLBI具有极高的相对精度和分辨率其广泛应用于天体物理，地球自转速率变化，主要是由日长变化来表征的。日长变长则地球自转速率减小，日长变长则地球自转速率变长。目前VBI对日长的测量精度总体可以优于0.05ms。目前许多的资料已经初步证明日长变化与大气角动量、、厄尔尼诺、拉尼娜等存在着紧密的联系。日长变化的研究对气候预测、防治、空间导航、深空探测、大地测量等领域有着重要的作用。如果我们能掌握日长变化的周期信息，影响日长变化的因素及这些因数的影响程度，我们就能够建立完善的日长变化预测的模型，对日长变化进行短期甚至长期的精确预报，这方面的研究具有具有重要的科学意义和实际应用价值。在目前预测模型方面,不论何种模型大都是先去除潮汐项的影响,然后在模型中考虑大气角动量等的影响对日长变化进行预测。笔者认为因为潮汐和大气角动量都是与月球、等星体对地球的影响和地球物质结构密切相关的，大气、海洋固体地球之间相互摩擦，相互产生影响，因此大气角动量与潮汐效应间也是相互影响相互作用的，他们之间可能相互依赖或者反映的是复杂系统的同一特征。在预测模型中加入相关性较强的因素进行预测会导致信息冗余且增加计算的的复杂程度，模型中两因素相关性较大时，虽然有时看起来提高了预测精度，但是这样的预测模型在原始数据精度提高时也可能会使预测结果更差[6。因此本文以月球对日长变化的影响为切入点分析地月距与日长变化的相关性，希望可以为以后日长变化的预测提供新的思路和参考。本文主要内第二章介绍VLBI测量及数据处理原理，日长变化基本原理并分析OCCAM软件 JPL包括精度分析，插值方法选择，基于最小二乘的趋势项去除，利用频谱分析和正弦曲线拟合进行周期提取，小波分析提取日长变化周期分量，日长变化与月地距的相关性分析等内容。2VLBI测量及数据处理原VLBI测量基本几何原甚长基线测量（VLBI）主要是指，相距为D具有独立本振的射电望远镜A和B同时对某些射电源进行观测，并将接受到的信号通过相关处理器进行处理，得到不同射电源信号到达A、B的时间延迟和时间延迟率，从而精确的测定射电源的位ABEOP参数的一整套利用VLBI技术测定射电源位置和基线向量的基本原理如图2.1所示（、、、彭碧波等2010）[7]。2.1VLBI观测原理（图错误某一时刻t射电望远镜A和B同时对某射电源进行观测。用于甚长基线测量的射电源与地球相距都在万光年以上，而射电望远镜的距离①1万多公λ/D(λ为射电波波长，D为基线长度)为到达两测站的波前面是垂直于信号方向的平面。那么在图2.1中当射电信号到B时也同时到达了a，该信号经过时间gA。我们称g为时间延迟，g①本文涉及的VLBI观测都是地基VLBI，测站都在地球上，空基VLBIgSaA

cosb1b1b

SaA2.1aAbABS射电源方向上的单位向量，θ为b和S我们通过对多个源的时间延迟g和时间延迟率g的多次测量，就可以通过国最小二乘或者卡尔曼滤波等方法得到b和源的位置。r0(x0,y0z0)rr0rtidersolrocnratm

式中rtide为潮汐对测站坐标的总影响，其由固体潮改正rsol、海潮负载改正rocn大气负载潮改正ratm、极潮改正

设射电望远镜AB在协议地球坐标系中的的实际位置矢量分别为rAxAyAzA和rBxByBzBb地球rB

2.2节中介绍到协议地球坐标系中坐标可以通过式r（PNSW）r地球协议天球坐标系中。式中上式中P表示岁差矩阵，N表示矩阵，S表示地球自转W表示极移矩阵。协议天球坐标系中的基线向量为：b

纬表示(αb,δb)射电源在协议天球坐标系的赤经和赤纬分别为αsδs。②cosbcosbcosscossbcossincossin c

b

s

问 VBI改正包括固体潮、海潮、极潮和大气负载潮等；另一类是针对信号环境的改正模型，主要包括，大气延迟改正（包括电离层、对流层等）、地球自转改正、系统②式（2.1）s为射电源指向测站的方向的单位向量，而式（2.6）中用到的为地心指向射电源的方向，因此设备引起的延迟改正（包括天线结构改正、电缆延迟改正、改正等）、射电源VLBI涉及的坐标系统坐标转换流阵球坐标阵球坐标阵球坐标阵球坐标阵阵球坐标2.2转换矩阵说2.3O-XYZITRS坐标系，O-XYZ为瞬时真地球坐标系，（xp,yp）2.3W(t)Rx(yp(t))Ry(xp

xp(t)和yp(tcosxp sinxp(t)Ry(xp

(t)) cosxp

Rx(yp

(t))

cosyp

sinyp

(t)

sinyp

cosyp

瞬时地球坐标系与瞬时天球坐标系Z轴都与地球瞬时自转轴重合。不同的是瞬时地球坐标系的X轴指向平子午线，而瞬时天球坐标系的X轴指向真春分点，两者间相差一个真春分点的恒星时角GAST，GAST可以理解为瞬时极2.4所示。2.4 0R(t)RZ(GAST) 0

图2.5中O-XYZ坐标系为瞬时天球坐标系，O-X'Y'Z'坐标系为瞬时平天球坐标系。Xγ0；X'γ；图中ε为黄赤交角；Δε为交角；ΔΨ为。IAU序列又678个日月项和687个行星项组成，IERS网上获得。图2.5改N(t)Rx()Rz()Rx(

R() cos(

cos(

R()sin(

sin()

R()

OCCAM6.2IAU1976IAU2000模型。本文的数据处理IAU2000模型。(McCarthyDDPetitG,2004）IAU2000岁差模型中，3A0.0000327''t4 .10.0000601''t4 0.0000470''t

2.6中，γ0J2000.0时的平春分点；Q0Q0'J2000.0时的平赤道；E0E0'为J2000.0时的平黄道；γ1t时刻的真春分点；Q1Q1't时刻的真赤道；E1E1'为J2000.0时刻的真黄道。坐标转换的岁差基本量：AAO0;AW(t)Rz(A)Ry(A)Rz(zA

cosW(t)sin

2.6

VLBI观测量估计日长变化原天文和大地参数解VBI建立对应的观测误差方程，利用最小二乘或卡尔曼滤波原理解算误差方程，得到所VBI的基本观测方程为：OtCX,t

式中Ot是在t时刻所获得的延迟和延迟率观测量；X是由与观测量有关的参数组成的向量CX,t为参数X对观测量Ot影响的数学模型，也成观测量的理论值。vt是观测量的噪声残差向量。设参数X由先验值和改正值x两部分组成，则对观测方程线性化有Ot

,tCX,

ytAtx

其中ytOtCXtt时刻观测值与理论值之差。At为偏导系数矩阵，它与

参数解算主要包括三个方面的工作：（1）建立计算延迟和延迟率的理论模型C(X,t)C(X0,t)；（2）At的计算；（3）参数改正值x的解算，可以采用最小二乘法或者卡尔曼滤波法。本文采用卡尔曼滤波日长变化计通常用日长LOD（lengthofday）的观测值相对于标准日长（86400国际时秒）的变化LOD来表示地球自转速率的变化[10]。通过对VLBI观测量进行解算便可得到相应的UT1。LOD的计算方法如(2.22) ，1990）[11]，当LOD大于零时，日长变短，地球自转速度加快，反之，日长变长地球自转速度减

OCCAMVLBI观测数据的OCCAM简VLBI数据处理运用最广的软件是Calc/SOLVK软件和OCCAM软件，由于OCCAM软件是不断更新和对外公开发布的③故本文的数据处理采用OCCAM软件。 是一种VLBI 数据处理分析软件,目前最新版本为OCCAM6.3,OCCAM_6.3与OCCAM6.2数据处理界面及方法并无大的变化，仅在结果文件的格式有较大变化,下面仅对OCCAM6.2作简要介绍。OCCAM6.2支持滤波、最小二乘（LSM）和最小二乘配置（LSCM）这三种平差方法。LSM和LSCM提供包括射电源历表、VLBI测站地理位置的时间序列、EOP和补偿在内参数的全球解，滤波接口主要用于EOP管理服务[12]。OCCAM6.2软件的模块结构2.7计算岁 数据准备GAUSS-MARKOV多基线KalmanlOCCAM6.2NGS2.3节已经介绍到测站坐标的参考框架为ITRF2000，射电源坐标的参考框架为ICRF。OCCAM6.2数据模型可采用IERS96或IERS2000标准，模型可以选择IAU1980或MBH2000，几何模型可IERS协议（1992）建立的模型或者IAU-1997模型。计算岁 数据准备GAUSS-MARKOV多基线Kalmanl获取观测-获取OCCAM图获取观测-获取OCCAM③OCCAM软件地址 .au/geodesy-EPHEMDAT，*EPHEMDAT，*文开数据准备：读入原始数据，生成OCCAM准数NGS,EOPXXXX,*CAT文输EPHEMDAT文EPHEMDAT文 架转换后的坐标数各类改正支持文件基线，测站坐标输测站改正微分：潮汐改正、对流层改正等系列改OCCAM标准数据库台输输BASTIM文OCCAM标准数据库台输滤波法进行多基线的数分 输SINEXSESSDAT，几何模型选择：计算每个原始数据文件中测值的理论延迟和延最小二乘配置全球解第一步:求解射电源置输最小二乘配置全球最小二乘配置全球解第二步:解算周日数，残差验后优化射电源位输STATISTSOURCESDAT等结果记录文OCCAM标准数据库台输INPUT8TXT文结结信息提2.8OCCAM6.2OCCAM解算参数和模型选用于VLBI观测波段主要有S波段（13cm）、P波段（92cm）、P波段（49cm）L波段（21cm）、L波段（18cm）、C波段（6cm）、X波段（3.6cm）、K（1.3cm）S/X双频波段主要用于大地测量学，L、CK等波段主要用于天体物理学，本文用到的是S/X波段的数据。我们从IVS获得的数据中通常有SA、SE、SH、SK、SU、SN、XE、XA、XH、XK、XF、XU、XNXSXS波段受到的影响要小得多，因此本文采取X波段的数据进行后续的分析。XUXK类型的数据只有23个测站的观测结果,不能满足本文的需求,XAXE这两种类型的结果精1XAXE这两种类型的观测数据进行后OCCAM_6.3OCCAM_6.2抽样计算结果表明，OCCAM_6.2计算得到1OCCAM_6.32.1所示。此外OCCAM_6.3系统处于调试阶段较不稳定故本文所用到的VLBI数据选择用OCCAM_6.22.1OCCAM软件处理结果方差因子对数据类处理的方差因处理的方差因两种计算海洋潮汐对极移和LOD影响的先验模型，RAY（GOT99.2b）EANES（CSR4.0）模型。可以选择的模型有IAU1980和MBH2000（等同于IAU2000A模型），为了从以上模型组合中选择更好的模型组合，本文选择了41.RAY+IAU19802.RAY+MBH20003.EANES+IAU1980EANES+MBH2000102.2所2.2CHI-SQδ（mas）LOD的标准差，CHI-SQ列的平均值12.2OCCAMRAY+IAUEANES+IAUCHI-CHI-CHI-CHI-0.1259从表2.2中可以看到从方差因子的角度考虑EANES+MBH2000组合最优，从标准差的角度考虑除了RAY+IAU1980组合较差，其他三种组合没有明显差距。基于EANES+MBH2000组合计算结果质量更高，本文的数据处理采用此组合模型。此外本文参考（，2003，2006）的研究结果[13]，[14]，选择提取周期小于35天短期变化和IERS推荐的日插值方法以得到更准确的U1插值。数据预处理完毕之后，本文采取多基线滤波法来进行多基线分析获取测站位置、EP以及补偿参数的最小二乘解。采取多基线滤波法时，已知参数包含固体潮参数、海洋荷载潮汐改正数、大气改正数和岁差参数，对流层以及钟补偿做随机处理。非随机参数包含VBIEP(地球定向参数和大气天顶延迟等。JPL星历简

3地距的计JPL（喷气推进）一直都在从事高精度的星历创建工作，JPL星历是主要用于空间导航和天文学的一种天文星历。JPLDE86、DE96、DE102、DE200、DE202、DE403、DE405、DE406、DE410、DE414、DE418、DE421等版3.1JPL3.1JPL星创建时基本情覆盖时备.5(-1410.04.16)J1950.5(3002.5(1599.12.09)J2000.5(1899.12.04)J2000.5(1599.04.29)ICRF.5(2199JED305424.50(1599.12.09)ICRFICRF.5(1900.02.06)ICRF.5(1899.12.04)ICRF.5(1899.12.04)\.5(1899.07.29)\.5(1899.07.29)ICRF.5(-3000.12.07)ICRF.5(3000JED378480.5(1799.12.16)ICRF2.0DE405IERS推荐的星历也是目前应用最为广泛星历。DE406DE405的结多项式的精度有所降低；DE403主要适用于行星和月球激光测距的数据处理；④①JED，JED（ET）47131112（JD）DE410主要用于火星探测车导航；DE413主要是为了更新冥王星轨道和研究冥王星最大的Charon而建立；DE414主要为了服务于火星探测而创建；DE421主要适用于行星和月球激光测距的数据处理；DE422MESSENGERDE406的补充，也承担一部分火星探测的星历服务；DE423MESSENGER探测器的探测任务。本文月地距的计算采用的DE405星历。JPL星历以每32天一组切插值多项式的形式将行星在系质心系中三维位置矢量提供给用户，每个星历记录包含一定时间跨度的切插值多项式系数。JPLACIID405ACII格式的数据可从ftp//d.jpl.naa./pub/ep/nets/ascii/e405/获，格式的数据每20年的数据。文件名为“ascSY.”，其中“YYYY表示该文件星历数据的起始年，"XX"JPL星历的REDME.txt[15说明文档。JPL星历使JPLFortran版、CJAVA版星历应用程序，CJAVA程序大部分版的星历应用程序，因此本文仅对Fortran版的星历应用程序做简要介绍[16]。基本文件获取及简在JPL的ftp服务器(ftp://ssd.jpl.nasa./）上获取一下地址的文件/pub/eph/nets/test-testeph.f——星历和插值程序。将这个这个程序的运行结果与JPL提供的检selcon.f——获取星历常数程序asc2eph.f——将ASCII格式的星历文件转化为二进制格式的星历文件的程序ascSYYYY.XXX——星历文件。这些以20asc2eph.f转换成二进制文件使用。转化后的二进制文件可以也可binmerge.f合并成一个单独的文件使用，也可以利用binshort.f从星历文件中提取一小部分单独使用。header.XXX——DEXXX的头文件信息，asc2eph．fASCII文件转化为二进制testpo.XXX——JPL的测试计算结果，它被testeph调用以检测星历安装是否正确，testpo.XXX文件。asc2eph.f、binmerge.fbinshort.f是创建和处理二进制星历文件的程序文件，此外用户从二进制星历文件中获得相应的信息时testeph．f中的子程序PLEPH、SELCON、CONST、DPLEPHSTATE是必须的。PLEPH主要作用是给定某一儒略日时刻，得到目标行对某行星中心的位置和速度；SELCON主要作用是获得读者所需的星历表某一常数；CONST主要作用是获取所有的有关星历表的常数；DPLEPHPLEPH相似，但输入的儒略日时刻更加精确。JPL星历的使用方199920113月的数据，因此需要用到的星历文月地距的计算方法[17]如下第一步：数据准备。首先要合并“header.405”、“ascp1980.405”和“ascp2000.405”文件，可以在WINDOWS的命令行界面输入“copyheader.405+ascp2000.405pre.405”产生ASCII格式的合并文件“pre.405”。然后将合并文件“pre.405”转换为二进制格式的“JPLEPH”文件。用asc2eph.f文件之前要修改a文件的一些内容，将NRECL定义并赋值注意要赋值成PARAMETER形式，本文NRECL赋值为4。可以在WINDOWS令行界面输入“DFasc2eph.f asc2eph<pre.405”产生第二步：修改testep.f并检查安装是否正确。在检测前先将FSIZE1和FSIZE2子函数的NAMFI、NRECL参数赋值。将FSIZE3子程序中各参数值做以下赋值“NEC=4，NFILE=12，NMFIL="JPLEP"，IZE2036”FSIZE3子程序。之后在WINDOWS命令行界面中输入“DFtesteph.f testeph<testpo.405”此时界面中会输出一系列的常用参数，以及和标准数据的比较结果，其中对应“difference”第三步：星历的使用。首先将selcon.f的全部代码和第二步中修改好的testeph.fFSIZER3、PLEPH、STATE、INTERP、SPLIT、STATE、CONST子程序到一个自己新建的fortran工程的源文件中。然后调用PLEP子程序得到目标行对X的M参数赋值为TRUE以得到km基本为单位的结果，否则得到的结果将以U为基本单位，如果需要U单位的结果可以不更改。JPL星历计算结为了与OCCCAM6.2处理得到的日长变化数据进行相关分析，计算2009年1月4.4（a）所示。4据分析方法、原理及结原始处理结果分本文利用IVS的2009至2010年NGS格式的数据⑤，使用4.14.1日长变化结果及精度统计（添加图例4.1(d)为日长变化值及其精度的线描绘。OCCAM6.24.1表 日长变化精度统计（单位（，2010）[14]的日长变化解算结果85%的数据解算精度在0.02ms以内，精度为0.0826ms。由此也可以说明本文的数据处理结果是有效的。从图4.1（c）⑤数据地址中可以明显看到日长变化的时间序列有明显的波峰和波谷的缺失，因此在进行下一步数据分析前需要对原始处理数据进行插值。插值方法的选择与原始数据插由于选取的VBI与计算，这样对于日长变化周期性和与月地距相关性的研究会产生较大的影响。需要对日长变化时间序列进行插值。插值方法的选IERS0910年以一天为时间间隔的日长变化（记为LOD）210210本文初选的插值方法有线性插值、三次多项式插值、三次样条插值。不同插值4.2所示。4.2不同差值方法效果对比（图d4.2(a)4.2(b)为三次样条插值法4.2(c)为三次多项式插值法得到的结果与原值较差序列；图4.2(d)为三种插值结果的对比图。从图4.2可以推断三次多项式插值较优(意！！！)，为了更清楚的分析不同差值方法的差别，将不同插值方法插值后与原4.24.2不同插值方法插值后与原值的差值（单位从图4.2和表4.2中可以看出线性插值法（Linear）得到的结果，三次多项优于0.2ms，也达到0.26ms可以满足精度要求（为什么？？？）。原始数据插值（目的4.2.1OCCAM6.2的计算结果进行插4.3所示。4.3OCCAM6.2_ΔLOD趋势项去趋势项是振动周期比信号采样长度大的频率成分，趋势项的存在会使相关函数、功率谱函数处理中出现变形，甚至可能使谱估计失真[18（JonhM.Whr，1988）[19总2ms203045s的周期波动。356,400km到406,700km之间变化，而且月球有逐渐远离地球的趋势。（J.S.Bandat和A..Piersol，199）[20对最小二乘4.2节中4.54.4图4.4(a)为月地距时间序列图；图4.4(b)为月地距趋势曲线；图4.4(c)为去趋势项4.4(d)为月地距去趋势项前后对比。从图4.4中可以看到在09至10年日地距离整体上有先增加再减小的过程，对于这两年月地距离总的来说是在增加的，去除趋势月地距离的数据在0值附近波动。图4.5(a)为三次样条插值后的OCCAM6.2处理的日长变化数据；图4.5(b)为日长变化趋势曲线；图4.5(c)去趋势项之前的日长变化数据与日长变化趋势的关系；图4.5(d)4.50910年日长变长的速度有所放缓，去除趋势日长变化的数据在0值附近波动。4.5日长变化趋势项提取和消除（图例周期项分在去除日长变化和月地距离的趋势，为了分析它们之间的相关性，先对二者的周期进行分析。本文采取频谱分析和三角曲线拟合两种周期分析方式（为什频谱分分析能够清楚揭示数据的频谱结构，通过时间序列的能量密度频谱，对4.6所示。图4.6(a)去趋势项前月地距的频谱图；图4.6(b)为去趋势月地距的时间序列图；图4.6(c)去趋势月地距的频谱图。从图4.6(a和图4.6(b的对比可以看出去除趋势项前月地距的周期并不能被发现，在去除趋势月地距变化的频谱图主频4.3中去趋势项对月地距周期分析的有效性。4.6从图4.6(c)中可读出能量峰值对应的频率值约为0.036023(1/天)，由于用于频谱分析的数据采样率为1天，所以可以得到09年至10年月地距变化周期为1/0.036023=27.76天。图4.7(a)为去趋势项前日长变化的频谱图；图4.7(b)为去趋势项后日长变化的时间序列图；图4.7(c)去趋势日长变化的频谱图。4.7图4.7(a)和图4.7(c)的对比可以看出在去趋势日长变化变化的频谱图主频更4.314.34.3频率(1/天5.7637*10-2.8818*10-周期(天091027.7827.3天是比较符合的。091013.61天、天、173.50347.0127.780910年月地距的周期在数值上吻合的很好，13.61天也接近于半月周期。曲线拟在4.1.1中采用了频谱分析的方法求周期，由于日长变化和月地距变化都有明显的正弦特性，所以本文利用CFtool工具箱对二者进行正弦曲线近。以得到不同方法求周期的结果。CFtool工具箱中的正弦曲线近，有8种类型，基础型是a1*sin(b1*x+c1)。曲线拟合效果的评价主要有SSE（残差平方和）、R-Square（可决系数）、AdjustedR-Square（调整的可决系数）、RMSE(均误差)四个指标[21]。SSESSERMSEniSSE(i

f(x

RMSE

nnTSSnTSS(yi2

R2R21时模型与样本观测值完全拟合，R21，模型的拟合优度越高。R21

（4）AdjustedAdjustedR2R21SSE(nkTSS(n

以式（4.1）至式（4.5）yif(xi)代表拟合值，n代表样本平均值n-k-1代表SSE的自由度。不同类型的正弦曲线对去除趋势日长变化的近结果评价指标统计如表4.3所示，*表示的为拟合最优的情况。4.3D12345678说明：根据表4.3中示的结果，去除趋势的日长变化选取第8种正弦拟（a1*in（b1*x1）+...+a8sin（b8*xc8））拟合，去除趋势的月地距采用第7种正弦拟合（a1*sin（b1*x+c1）+...+a8*sin（b7*x+c7））拟合最优。拟4.84.4图4.8（a）和（b）分别为去趋势的日长边和和月地距的正弦拟合。表4.4中bi项的第二列为由bi的值求得的周期，在ai取值过小时bi项中不对其对应周期进4.4376.24天、180.29天、13.6313.64天、27.57天五个明显的周期。取13.63天和13.64天平均为13.64天。4.4.1中频谱分析的结果有一定的差异，日长变化的月周期和半月周期与频谱分析的结果差异较小，尤其是半周期项。存在这种差异的原因与二者的基本原理有关频谱分析基本原理是变换，正弦曲线拟合的基本原理是最小二乘，此外采样率不高也是导致分析结果有所偏差的原因。4.84.41234-5/-6-7/8/-D1-2/3/4/5-/6/7/不论是何种求周期的方法都得出了日长变化存在与月地距离变化很相符的月周期和半月周期项，基于此推断日长变化的月周期分量和半月周期分量与月地距离存在相关性。相关性分4.4析，由于考虑到可能存在时滞的相关性，本节采用互相关分析，得到最大相关系数的同时获取时滞信息。日长变化周期分量提小波分析是一种窗口大小不变形状可变的时频局域化分析方法，可以将不同频率的高低频信号区分出来。因此可以用于周期分量的提取，然而小波变换提取周期分量会有混频的现象，因此要根据不同小波基的特点进行选择。经过不断的测试调6demy4.9。4.9demy法读出D3层的周期约为1/0.0735=13.61天；D4层的周期约为1/0.036=27.78天。6dmey的。由于月地距的周期项很明显，为了减少在剔除高频项对月地距产生的负影响，对于月地距离不做上述周期项提取。互相关分互相关函数可以描述两列时间序列X(t),Y()在任意两个不同时刻t，t2间的相关程度。相对于相关函数在本文的相关分析中更适用。查相关系数显著性检验表可以查到在样本量为7000.13499.9%的。图4.10各层小波时频分析（注意纵坐标标值匹配性日长变化月周期分量与去趋势月地距的相关在图4.11(b)中我们也可以明显看到，随着月地距的增加，日长变化月周期分量化的月周期分量（27.78天）与去趋势月地距用进行互相关分析，分析结果显示，在月地距离(D)滞后日长变化的月周期分量(MLOD)0天时有最大负相关-0.822，在D滞后MLOD一天时有次大负相关-0.803，因此日长变化的月周期分量与去趋势的99.9%的置信水平下是显著负相关的。日长变化半月周期分量与去趋势月地距的相关直接将日长变化半月周期分量(SMLOD)与去趋势月地距(D)进行互相关分析4.114.11(a)为日长变化月周期分量（MLOD）与月地距（D）时间序列；图4.11(b)为图4.11(a)的细节放大；图4.11(c)为日长变化半月周期分量与月地距（D）时间序列；图4.11(d)为图4.11(c)的细节放大；4.11(e)为日长变化半月周期分量（MLOD）与月地距变化速率绝对值（VS）时间序列；4.11(f)4.11(e)的细节放大。从图4.11()中可以看到日长变化的半月周期分量大致存在着在日地距离变化速率最快时出现周期内最大值的情况。为了验证这一推论，我们对日长变化求变化速率的绝对值，将其与日长变化的半月周期项进行互相关分析。地月距离变化绝对速30.564.势的月地距的变化速率绝对值在99.9%的置信水平下时显著正相关的。日长变的月周期分量与去趋势的月地距的变化速率相关性没有日长变化月周期分量与去趋势月地距的相关性高。然而儒略日54836到54999有时滞1天的最大相关数.765550005520020.738系数的值产生了影响。笔者认为，因为日长变化对月地距变化的反应存在延时，而这种延时又与地球的内外部环境相关，针对不同的互相关量在不同的时期有不同的时滞。（赵凯华，1996）提出“大量地球物理观测表明，地球对力的响应并是纯弹性的，而是滞弹性的，即应变稍有延迟。"[22]这可以在一定程度解释日长变化与月地距而变化的时滞现象。结

5论与展OCCAM6.3OCCAM6.2相同条件下处理结果方差因子的对比，选择处理结果较优的OCCAM6.2（值的真实精度外比较？？？）。通过不同模型组合处理