理论固体潮_matlab
- 引言
固体潮是指在日、月引潮力的作用下,固体地球产生的周期形变的现象。固体潮的存在说明固体地球具有一定的弹性,固体潮就是弹性地球在日月引力作用下发生的弹性变形。引潮力不但可以引起地球表面流体的潮汐(如海潮、大气潮),还能引起地球固体部分的周期性形变。受固体潮的影响,地面不停的变形,影响到各种测量数据的精确度。精密大地测量结果应加入相应的修正。重力观测精度已达到2~5微伽的量级,而重力潮汐变化影响的最大幅度可达±130 微伽 ,故须加入改正。卫星激光测距精度达到3厘米 ,而地面测站的垂直潮汐形变达到30~40厘米的幅度,必须加以改正。固体潮的变化对卫星轨道也有摄动作用,所以在卫星轨道设计中必须顾及这一影响。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-UnAHE8aE-1630728376695)(https://dspstack.com/uploads/images/201901/10/68/CdSsgBTmNS.jpg)]
- 理论固体潮计算
固体潮是一种可以预先精确计算得到的物理量,基本计算过程可分为四步。
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;
- 参考维基百科、百度百科、郜晓亮等重力固体潮理论计算文章