重力正反演程序

function g=forward(x,T) %二度板状体异常子程序
x0=T(1);z0=T(2);b=T(3);l=T(4);a0=T(5);Q=T(6);%二度板状体各个参数
a=(180-a0)*pi/180;%转化成台阶的角度
f=6.67*10^(-11);%引力常数
h=z0-l*sin(a0*pi/180);%上顶面深度
H=z0+l*sin(a0*pi/180);%下底面深度
for i=1:length(x)
%左台阶异常
X(i)=x(i)-(z0*cot(a)+(x0-b));%坐标转换
p1(i)=pi*(H-h)+2*H*atan((X(i)+H*cot(a))/H);
p2(i)=-2*h*atan((X(i)+h*cot(a))/h);
p31(i)=(H+X(i)*sin(a)*cos(a))^2+X(i)^2*(sin(a))^4;
p32(i)=(h+X(i)*sin(a)*cos(a))^2+X(i)^2*(sin(a))^4;
p3(i)=X(i)*(sin(a))^2*log(p31(i)/p32(i));
p41(i)=X(i)*(H-h)*(sin(a))^2;
p42(i)=X(i)^2*(sin(a))^2+(H+h)*X(i)*sin(a)*cos(a)+H*h;
p4(i)=-2*X(i)*sin(a)*cos(a)*atan(p41(i)/p42(i));
g1(i)=f*Q*(p1(i)+p2(i)+p3(i)+p4(i));%理论计算公式
%右台阶异常
X(i)=x(i)-(z0*cot(a)+(x0+b));%坐标转换
p1(i)=pi*(H-h)+2*H*atan((X(i)+H*cot(a))/H);
p2(i)=-2*h*atan((X(i)+h*cot(a))/h);
p31(i)=(H+X(i)*sin(a)*cos(a))^2+X(i)^2*(sin(a))^4;
p32(i)=(h+X(i)*sin(a)*cos(a))^2+X(i)^2*(sin(a))^4;
p3(i)=X(i)*(sin(a))^2*log(p31(i)/p32(i));
p41(i)=X(i)*(H-h)*(sin(a))^2;
p42(i)=X(i)^2*(sin(a))^2+(H+h)*X(i)*sin(a)*cos(a)+H*h;
p4(i)=-2*X(i)*sin(a)*cos(a)*atan(p41(i)/p42(i));
g2(i)=f*Q*(p1(i)+p2(i)+p3(i)+p4(i));%理论计算公式
g(i)=g1(i)-g2(i);%二度板状体异常（左台阶与右台阶异常差）
g(i)=10^6*g(i);%将重力异常的单位转化成高斯
end



clear all
close all
x=-200:5:200; %观测点横坐标
z=0; %观测点纵坐标
T=[0,50,20,10,45,500];%给出正演模型参数T
g=forward(x,T);%调用重力异常子程序
figure(1)
plot(x,g,'.-g') %画出重力异常
ylabel('\Deltag/g.u.','FontSize',12)
xlabel('测点横坐标/m','FontSize',12)


clear all
close all
x=-200:5:200; %测点横坐标
z=0; %观测点纵坐标
T=[0,50,20,10,45,500];%给出模型参数T
g=forward(x,T);%调用重力异常子程序
gn=g;
n=5;%含n%的随机误差
noise=2*rand(1,length(x))-1;
gn=gn+gn*n.*noise/100;
figure(2)
plot(x,g,'g')
hold on
plot(x,gn,'.-r')
ylabel('\Deltag/g.u.','FontSize',12)
xlabel('测点横坐标/m','FontSize',12)
legend('无误差重力异常','含5%随机误差重力异常');



=[0,50,20,10,45,500];%实际模型参数
T=0.5*(t+[2,30,18,12,35,520]);%反演初始模型参数
x=-200:5:200;%测点横坐标
z=0;%观测点纵坐标
m=1;%阻尼因子
v=100;%阻尼因子调整系数
smax=20;%最大迭代次数
for s=1:smax %循环迭代
s=s
gm=forward(x,T);%由模型参数计算出的重力异常
dg=gn-gm;%观测数据与模型重力异常差值
for j=1:length(T) %求出雅可比矩阵P
dT=10^(-2)*T(j);
T(j)=T(j)+dT;
aa=forward(x,T);
T(j)=T(j)-dT;
bb=forward(x,T);
P(:,j)=(aa-bb)/dT;
end
A=P'*P;q=P'*dg';D=diag(diag(A));%线性方程组（(A+mD)PT=q）的

系数矩阵
[U,S,V]=svd(A+m*D);%奇异值分解法(svd)
for i=1:size(A,1) %限定S矩阵的奇异值范围
if (S(i,i)<10^(-12))
S(i,i)=10^(-12);
end
if (S(i,i)>10^(12))
S(i,i)=10^(12);
end
end
PT=inv(V')*inv(S)*inv(U)*q;%奇异值分解法(svd)求解线性方程组
T=T+PT';%改正模型参数
gm=forward(x,T);
F(s)=0;%初始化反演重力异常总误差
for i=1:length(x)
F(s)=F(s)+abs(gn(i)-gm(i))/abs(gm(i));%求反演重力异常总误差
end
F(s)=F(s)/length(x);%求反演重力异常平均误差
if s==smax%超过设定最大迭代次数，跳出循环
break;
end
if s>=3&&F(s)>=F(s-2)&&F(s-1)>=F(s-2)%反演重力异常平均误差不收敛，跳出循环
break;
end
if F(s)<=5/100 %达到反演重力异常平均误差容许值，跳出循环
break;
end
if s>1&&F(s)<=F(s-1)%调整阻尼因子
m=m/v;
elseif s>1&&F(s)>F(s-1)
m=m*v;
end
end
figure(3)
semilogy(1:s,F(1:s)*100);
ylabel('反演重力异常的平均相对误差/%','FontSize',10);
xlabel('迭代次数s','FontSize',10);
figure(4)
plot(x,gn);
hold on
plot(x,gm,'.-r');
ylabel('\Deltag/g.u.','FontSize',15);
xlabel('测点横坐标/m','FontSize',15);
legend('实测重力异常','反演模型正演重力异常');










%二度异常向上延拓(利用剖面上的异常值)matlab7.4计算程序
%Design by 赵玉岩
%最后修改日期：2008年4月21日23:03:19
%Δg原始数据为excel文件，内容格式为：
%点号 Δg
%-2 0.00
%-1 0.00
%0 0.00
%1 0.00
%2 0.00
%Bi原始数据为，内excel文件容格式为：
%Bi
%0.00
%0.00
%0.00
%##########################################################################
clear;
clc;%清空窗体和内存内容
%以下读入相关参数的对话框
%以下是读入点个数的对话框
cydj=input('输入采样点距') %采样点距
ytgd=input('输入延拓高度') %延拓高度
if ytgdmsgbox '延拓高度小于采样点距，程序无法工作'
break
end
if rem(ytgd,cydj)~=0
msgbox '延拓高度不是采样点距的整数倍，程序无法工作'
break
end
[datafile,pathname]=uigetfile({'*.xls';'*.*'},'打开Δg原始数据文件：');%打开Δg原始数据文件
[ydata headtext]= xlsread([pathname datafile]);
subplot(2,1,1),plot(ydata(:,1),ydata(:,2))
hold on
ydx=size(ydata);%矩阵大小
yhs=ydx(1);%矩阵行数
ydh=ydata(:,1);
ydg=ydata(:,2);
[datafile,pathname]=uigetfile({'*.xls';'*.*'},'打开系数原始数据文件：');%打开系数原始数据文件
[bdata headtext]= xlsread([pathname datafile]);
bdx=size(bdata);%矩阵大小
bhs=bdx(1)-1;%矩阵行数
if 2*bhs>yhs
msgbox '系数值个数大于测点的1/2，程序无法工作'
break
end
%

对系数进行处理
bn=[];
for i=0:bhs-1
bn=[bn bdata(bhs-i)];
end
bn=[bn bdata'];
%ydg进行扩边处理
nydg=td(ydg',ytgd/cydj*bhs);
ydg=nydg';
yhs=yhs+ytgd/cydj*bhs*2;
%以下开始计算,不包括边缘数据
ytdgjz=[];
for i=(ytgd/cydj*bhs+1):(yhs-ytgd/cydj*bhs)
ndg=[];
for j=(i-ytgd/cydj*bhs):(ytgd/cydj):(i+ytgd/cydj*bhs)
ndg=[ndg ydg(j)];
end
ytdg=ndg*(bn');
ytdgjz=[ytdgjz ytdg];
end
[filename2, pathname2] = uiputfile( ...
{'*.xls', '电子表格(*.xls)'}, ...
'二度异常向上延拓计算结果：');
a1={'点号' '向上延拓值'};
xlswrite([pathname2 filename2],a1,'sheet1','A1');
xlswrite([pathname2 filename2], ydh ,'sheet1','A2');
xlswrite([pathname2 filename2], ytdgjz' ,'sheet1','B2');
subplot(2,1,2),plot(ytdgjz)
msgbox '计算完成！'
hold off