理论固体潮计算_matlab

原作者: [db:作者] 来自: [db:来源] 收藏邀请

理论固体潮_matlab

引言
固体潮是指在日、月引潮力的作用下，固体地球产生的周期形变的现象。固体潮的存在说明固体地球具有一定的弹性，固体潮就是弹性地球在日月引力作用下发生的弹性变形。引潮力不但可以引起地球表面流体的潮汐（如海潮、大气潮），还能引起地球固体部分的周期性形变。受固体潮的影响，地面不停的变形，影响到各种测量数据的精确度。精密大地测量结果应加入相应的修正。重力观测精度已达到2～5微伽的量级，而重力潮汐变化影响的最大幅度可达±130 微伽，故须加入改正。卫星激光测距精度达到3厘米，而地面测站的垂直潮汐形变达到30～40厘米的幅度，必须加以改正。固体潮的变化对卫星轨道也有摄动作用，所以在卫星轨道设计中必须顾及这一影响。

理论固体潮计算
固体潮是一种可以预先精确计算得到的物理量，基本计算过程可分为四步。
1、地理纬度转换为地心纬度，弧度制。
2、求儒略世纪数T
3、由儒略数计算6个天文参数（月球平黄经S、月球深交点黄经N、月球近地点黄经P、太阳平黄经h、太阳近地点黄经Ps、黄赤交角sigma）
4、根据固体潮公式带人计算。
之前从甚宽频带地震计中提取固体潮信息的工作时，简单整理写了一下计算理论固体潮的脚本，如下：
计算固体潮函数gravity_earthtide

function deltag=gravity_earthtide(ymdh,Longtitude,latitude)
% 理论固体计算函数
% output 
% 计算时间[2000,1,1,1,0,0];
% 经度
% 纬度
% 2015/10/22
% timevec = datevec(daynum);
y = ymdh(1);
m = ymdh(2);
d = ymdh(3);
t = ymdh(4);
Phi = latitude;     
%换算
Phip = Phi-0.192424*sin(2*Phi*pi/180.0);   
Longtitude = Longtitude*pi/180.0;
Latitude = latitude*pi/180.0;
Phi = Latitude;
Phip = Phip*pi/180.0;
To = Julian(y,m,d);
T = (To-2415020+(t-8)/24)./(36525);
% 天文学参数
S = 270.43416+481267.89057*T-0.001133*T*T+0.000002*T*T*T;
H = 279.696678+36000.76892*T+0.0003025*T*T;
P = 334.3295556+4069.034033*T+0.0001*T*T*T;
N = 259.183275-1934.142008*T+0.002077*T*T+0.000002*T*T*T;
Ps = 281.2208333+1.719175*T+0.0004527*T*T+0.000003*T*T*T;
sigma = 23.452294-0.0130125*T-0.000002*T*T;
S = S*pi/180;
h = H*pi/180;
p = P*pi/180;
N = N*pi/180;
ps = Ps*pi/180;
e = sigma*pi/180;
crm = 1+0.0545*cos(S-p)+0.0030*cos(2*(S-p))+0.01*cos(S-2*h+p)+0.0082*cos(2*(S-h))...
+0.0006*cos(2*S-3*h+ps)+0.0009*cos(3*S-2*h-p);
Lambdam = S+0.0222*sin(S-2*h+p)+0.1098*sin(S-p)+0.0115*sin(2*S-2*h)+0.0037*sin(2*S-2*p)...
        -0.0032*sin(h-ps)-0.001*sin(2*h-2*p)+0.001*sin(S-3*h+p+ps)+0.0007*sin(S-h-p+ps)...
        -0.0006*sin(S-h)-0.0005*sin(S+h-p-ps)+0.0008*sin(2*S-3*h+ps)-0.002*sin(2*S-2*N)...
        +0.0009*sin(3*S-2*h-p);
Beltam = 0.003*sin(S-2*h+N)+0.0895*sin(S-N)+0.0049*sin(2*S-p-N)-0.0048*sin(p-N)...
    -0.0008*sin(2*h-p-N)+0.001*sin(2*S-2*h+p-N)+0.0006*sin(3*S-2*h-N);
Delta = sin(e)*sin(Lambdam)*cos(Beltam)+cos(e)*sin(Beltam);
Theta = ((t-8)*(15*pi/180.0))+h+Longtitude-pi;
H = cos(Beltam)*cos(Lambdam)*cos(Theta)+sin(Theta)*(cos(e)*cos(Beltam)*sin(Lambdam)-sin(e)*sin(Beltam));
Zm = sin(Phip)*Delta+cos(Phip)*H; 
crs = 1+0.0168*cos(h-ps)+0.0003*cos(2*h-2*ps);
Lambdas = h+0.0335*sin(h-ps)+0.0004*sin(2*h-2*ps);
Beltas = 0.0;
Zs = sin(Phip)*sin(e)*sin(Lambdas)+cos(Phip)*(cos(Lambdas)*cos(Theta)+sin(Theta)*cos(e)*sin(Lambdas)); 
F = 0.998327+0.001676*cos(2*Phi);
Gt = -165.17*F*crm*crm*crm*(Zm*Zm-1.0/3.0)-1.3708*F*F*crm*crm*crm*crm*Zm*(5*Zm*Zm-3)-76.08*F*crs*crs*crs*(Zs*Zs-1.0/3.0);
deltag = Gt*1.16;
end

儒略日计算Julian.m

function R=Julian(y,m,d)
% 儒略数
    yy = y-1900;
    mm = m-1;
    dd = d;
    w = floor(yy/4);
    if(y == 4*w & mm<2)
        dd=dd-1;
    end
    if mm == 0
        d1 = yy*365+w-0.5+mm+dd;
    else if (mm == 1)
            mm = 31;
            d1 = yy*365 + w - 0.5 + mm + dd;
        else
            mm = floor(mm*365/12)-floor(10/(4+mm));
            d1 = yy*365 + w - 0.5 + mm + dd;
        end
    end
    R=d1+2415020.0;
end

计算例子main.m,计算2019.1.1-2019.1.10

clc
clear;
Longtitude=104;   Latitude=30.6;
date1 = [2019,1,1];
date2 = [2019,1,10];
amend= [];
kk = 1;
for i = 1:datenum(date2)-datenum(date1)
    date = datenum(date1) + i;
    for j = 0:23
        ymdh = datevec(date);
        ymdh(4) = j;
        ymdhnum(kk) = datenum(ymdh);
        amend = [amend,gravity_earthtide(ymdh,Longtitude,Latitude)];
        kk = kk+1;
    end
end
disp('========= complete =========');
plot(ymdhnum,amend)
datetick;