• bwoah
    了解作者
  • matlab
    开发工具
  • 2KB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 0
    下载次数
  • 2020-11-11 10:10
    上传日期
三维规则网格剖分下的重磁正演计算,得到重力异常分布
新建文件夹.zip
  • 新建文件夹
  • forward_gm3d.m
    10KB
内容介绍
function [sens_g,sens_m]=forward_gm3d(Nx,Ny,Nz,deltx,delty,deltz) %重磁三维正演,单位分别为mGal、nT (g/cm3 SI)正演磁化率 %Nx、Ny、Nz分别为x、y、z三个方向上网格剖分的数量,deltx、delty、deltz分别为x、y、z三个方向上网格剖分的间隔 %sens_g为重力灵敏度矩阵sens_g(Nx*Ny,Nx*Ny*Nz) sens_m为磁法灵敏度矩阵sens_m(Nx*Ny,Nx*Ny*Nz) %排列顺序station:x1,y1 x1,y2 x1,y3 ... x1,yn x2,y1 ... %model:x1,y1,z1 x1,y2,z1 x1,y3,z1 ... x1,yn,z1 x2,y1,z1 ... xn,yn,z1 x1,y1,z2 ... Stax=0.5*deltx:deltx:(Nx-0.5)*deltx; Stay=0.5*delty:delty:(Ny-0.5)*delty; [Sta_x,Sta_y]=meshgrid(Stax,Stay); Stah=-2;%观测点高程 Sta=[reshape(Sta_x,[Nx*Ny,1]),reshape(Sta_y,[Nx*Ny,1]),zeros(Nx*Ny,1)+Stah];%观测点坐标 Gc=6.6732e-3;%万有引力常数 T=50000;%地磁场强度,unit nT I=degtorad(90);%地磁倾角 A1=degtorad(0);%测线方位角,测线方向x轴与磁北夹角 A=0.5*pi-A1; sens_g=zeros(Nx*Ny,Nx*Ny*Nz); Hax=zeros(Nx*Ny,Nx*Ny*Nz); Hay=zeros(Nx*Ny,Nx*Ny*Nz); Za=zeros(Nx*Ny,Nx*Ny*Nz); sens_m=zeros(Nx*Ny,Nx*Ny*Nz); I=-I; Trax=cos(I)*sin(A); Tray=cos(I)*cos(A); Traz=sin(I); for i=1:Nx*Ny k=1; for Zm=1:Nz for Xm=1:Nx for Ym=1:Ny Xm1=(Xm-1)*deltx; Xm2=Xm*deltx; Ym1=(Ym-1)*delty; Ym2=Ym*delty; Zm1=(Zm-1)*deltz; Zm2=Zm*deltz; R1=norm(Sta(i,:)-[Xm2 Ym2 Zm2]); R2=norm(Sta(i,:)-[Xm1 Ym2 Zm2]); R3=norm(Sta(i,:)-[Xm2 Ym1 Zm2]); R4=norm(Sta(i,:)-[Xm1 Ym1 Zm2]); R5=norm(Sta(i,:)-[Xm2 Ym2 Zm1]); R6=norm(Sta(i,:)-[Xm1 Ym2 Zm1]); R7=norm(Sta(i,:)-[Xm2 Ym1 Zm1]); R8=norm(Sta(i,:)-[Xm1 Ym1 Zm1]); % R1=norm(Sta(i,:)-[Xm2-0.5*deltx Ym2-0.5*delty Zm2-0.5*deltz]); % R2=R1; % R3=R1; % R4=R1; % R5=R1; % R6=R1; % R7=R1; % R8=R1; if Xm1==0 Xm1=0.00001; end if Ym1==0 Ym1=0.00001; end % if Zm1==0 % Zm1=0.0000001; % end sens_g(i,k)=-Gc*( ((Xm2-Sta(i,1))*log((Ym2-Sta(i,2))+R1)+(Ym2-Sta(i,2))*log((Xm2-Sta(i,1))+R1)-(Zm2-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Ym2-Sta(i,2)))/((Zm2-Sta(i,3))*R1))) ... -((Xm1-Sta(i,1))*log((Ym2-Sta(i,2))+R2)+(Ym2-Sta(i,2))*log((Xm1-Sta(i,1))+R2)-(Zm2-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Ym2-Sta(i,2)))/((Zm2-Sta(i,3))*R2))) ... -((Xm2-Sta(i,1))*log((Ym1-Sta(i,2))+R3)+(Ym1-Sta(i,2))*log((Xm2-Sta(i,1))+R3)-(Zm2-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Ym1-Sta(i,2)))/((Zm2-Sta(i,3))*R3))) ... +((Xm1-Sta(i,1))*log((Ym1-Sta(i,2))+R4)+(Ym1-Sta(i,2))*log((Xm1-Sta(i,1))+R4)-(Zm2-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Ym1-Sta(i,2)))/((Zm2-Sta(i,3))*R4))) ... -((Xm2-Sta(i,1))*log((Ym2-Sta(i,2))+R5)+(Ym2-Sta(i,2))*log((Xm2-Sta(i,1))+R5)-(Zm1-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Ym2-Sta(i,2)))/((Zm1-Sta(i,3))*R5))) ... +((Xm1-Sta(i,1))*log((Ym2-Sta(i,2))+R6)+(Ym2-Sta(i,2))*log((Xm1-Sta(i,1))+R6)-(Zm1-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Ym2-Sta(i,2)))/((Zm1-Sta(i,3))*R6))) ... +((Xm2-Sta(i,1))*log((Ym1-Sta(i,2))+R7)+(Ym1-Sta(i,2))*log((Xm2-Sta(i,1))+R7)-(Zm1-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Ym1-Sta(i,2)))/((Zm1-Sta(i,3))*R7))) ... -((Xm1-Sta(i,1))*log((Ym1-Sta(i,2))+R8)+(Ym1-Sta(i,2))*log((Xm1-Sta(i,1))+R8)-(Zm1-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Ym1-Sta(i,2)))/((Zm1-Sta(i,3))*R8))) ); % -(Xm1*log(Ym2+R2)+Ym2*log(Xm1+R2)-Zm2*atan((Zm2*R2)/(Xm1*Ym2))) ... % -(Xm2*log(Ym1+R3)+Ym1*log(Xm2+R3)-Zm2*atan((Zm2*R3)/(Xm2*Ym1)))+(Xm1*log(Ym1+R4)+Ym1*log(Xm1+R4)-Zm2*atan((Zm2*R4)/(Xm1*Ym1))) ... % -(Xm2*log(Ym2+R5)+Ym2*log(Xm2+R5)-Zm1*atan((Zm1*R5)/(Xm2*Ym2)))+(Xm1*log(Ym2+R6)+Ym2*log(Xm1+R6)-Zm1*atan((Zm1*R6)/(Xm1*Ym2))) ... % +(Xm2*log(Ym1+R7)+Ym1*log(Xm2+R7)-Zm1*atan((Zm1*R7)/(Xm2*Ym1)))-(Xm1*log(Ym1+R8)+Ym1*log(Xm1+R8)-Zm1*atan((Zm1*R8)/(Xm1*Ym1))) ); % sens_g(i,k)=-Gc*( ( (Xm2-Sta(i,1))*log(Ym2-Sta(i,2)+R1)+(Ym2-Sta(i,2))*log(Xm2-Sta(i,1)+R1)+(Zm2-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Zm2-Sta(i,3)))/((Ym2-Sta(i,2))^2+(Zm2-Sta(i,3))^2+(Ym2-Sta(i,2))*R1))) ... % -( (Xm1-Sta(i,1))*log(Ym2-Sta(i,2)+R2)+(Ym2-Sta(i,2))*log(Xm1-Sta(i,1)+R2)+(Zm2-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Zm2-Sta(i,3)))/((Ym2-Sta(i,2))^2+(Zm2-Sta(i,3))^2+(Ym2-Sta(i,2))*R2))) ... % -( (Xm2-Sta(i,1))*log(Ym1-Sta(i,2)+R3)+(Ym1-Sta(i,2))*log(Xm2-Sta(i,1)+R3)+(Zm2-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Zm2-Sta(i,3)))/((Ym1-Sta(i,2))^2+(Zm2-Sta(i,3))^2+(Ym1-Sta(i,2))*R3))) ... % +( (Xm1-Sta(i,1))*log(Ym1-Sta(i,2)+R4)+(Ym1-Sta(i,2))*log(Xm1-Sta(i,1)+R4)+(Zm2-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Zm2-Sta(i,3)))/((Ym1-Sta(i,2))^2+(Zm2-Sta(i,3))^2+(Ym1-Sta(i,2))*R4))) ... % -( (Xm2-Sta(i,1))*log(Ym2-Sta(i,2)+R5)+(Ym2-Sta(i,2))*log(Xm2-Sta(i,1)+R5)+(Zm1-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Zm1-Sta(i,3)))/((Ym2-Sta(i,2))^2+(Zm1-Sta(i,3))^2+(Ym2-Sta(i,2))*R5))) ... % +( (Xm1-Sta(i,1))*log(Ym2-Sta(i,2)+R6)+(Ym2-Sta(i,2))*log(Xm1-Sta(i,1)+R6)+(Zm1-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Zm1-Sta(i,3)))/((Ym2-Sta(i,2))^2+(Zm1-Sta(i,3))^2+(Ym2-Sta(i,2))*R6))) ... % +( (Xm2-Sta(i,1))*log(Ym1-Sta(i,2)+R7)+(Ym1-Sta(i,2))*log(Xm2-Sta(i,1)+R7)+(Zm1-Sta(i,3))*atan(((Xm2-Sta(i,1))*(Zm1-Sta(i,3)))/((Ym1-Sta(i,2))^2+(Zm1-Sta(i,3))^2+(Ym1-Sta(i,2))*R7))) ... % -( (Xm1-Sta(i,1))*log(Ym1-Sta(i,2)+R8)+(Ym1-Sta(i,2))*log(Xm1-Sta(i,1)+R8)+(Zm1-Sta(i,3))*atan(((Xm1-Sta(i,1))*(Zm1-Sta(i,3)))/((Ym1-Sta(i,2))^2+(Zm1-Sta(i,3))^2+(Ym1-Sta(i,2))*R8))) ); Hax(i,k)=(-Trax*atan((Xm2-Sta(i,1))*(Ym2-Sta(i,2))/((Xm2-Sta(i,1))^2+(Zm2-Sta(i,3))*R1+(Zm2-Sta(i,3))^2))+Tray*log(R1+(Zm2-Sta(i,3)))+Traz*log(R1+(Ym2-Sta(i,2)))) ... -(-Trax*atan((Xm1-Sta(i,1))*(Ym2-Sta(i,2))/((Xm1-Sta(i,1))^2+(Zm2-Sta(i,3))*R2+(Zm2-Sta(i,3))^2))+Tray*log(R2+(Zm2-Sta(i,3)))+Traz*log(R2+(Ym2-Sta(i,2)))) ... -(-Trax*atan((Xm2-Sta(i,1))*(Ym1-Sta(i,2))/((Xm2-Sta(i,1))^2+(Zm2-Sta(i,3))*R3+(Zm2-Sta(i,3))^2))+Tray*log(R3+(Zm2-Sta(i,3)))+Traz*log(R3+(Ym1-Sta(i,2)))) ... +(-Trax*atan((Xm1-Sta(i,1))*(Ym1-Sta(i,2))/((Xm1-Sta(i,1))^2+(Zm2-Sta(i,3))*R4+(Zm2-Sta(i,3))^2))+Tray*log(R4+(Zm2-Sta(i,3)))+Traz*log(R4+(Ym1-Sta(i,2)))) ... -(-Trax*atan((Xm2-Sta(i,1))*(Ym2-Sta(i,2))/((Xm2-Sta(i,1))^2+(Zm1-Sta(i,3))*R5+(Zm1-Sta(i,3))^2))+Tray*log(R5+(Zm1-Sta(i,3)))+Traz*log(R5+(Ym2-Sta(i,2)))) ... +(-Trax*atan((Xm1-Sta(i,1))*(Ym2-Sta(i,2))/((Xm1-Sta(i,1))^2+(Zm1-Sta(i,3))*R6+(Zm1-Sta(i,3))^2))+Tray*log(R6+(Zm1-Sta(i,3)))+Traz*log(R6+(Ym2-Sta(i,2)))) ... +(-Trax*atan((Xm2-Sta(i,1))*(Ym1-Sta(i,2))/((Xm2-Sta(i,1))^2+(Zm1-Sta(i,3))*R7+(Zm1-Sta(i,3))^2))+Tray*log(R7+(Zm1-Sta(i,3)))+Traz*log(R7+(Ym1-Sta(i,2)))) ... -(-Trax*atan((Xm1-Sta(i,1))*(Ym1-Sta(i,2))/((Xm1-Sta(i,1))^2+(Zm1-Sta(i,3))*R8+(Zm1-Sta(i,3))^2))+Tray*log(R8+(Zm1-Sta(i,3)))+Traz*log(R8+(Ym1-Sta(i,2)))) ; Hay(i,k)=(Trax*log(R1+(Zm2-Sta(i,3)))-Tray*atan((Xm2-Sta(i,1))*(Ym2-Sta(i,2))/((Ym2-Sta(i,2))^2+R1*(Zm2-Sta(i,3))+(Zm2-Sta(i,3))^2))+Traz*log(R1+(Xm2-Sta(i,1)))) ... -(Trax*log(R2+(Zm2-Sta(i,3)))-Tray*atan((Xm1-Sta(i,1))*(Ym2-Sta(i,2))/((Ym2-Sta(i,2))^2+R2*(Zm2-Sta(i,3))+(Zm2-Sta(i,3))^2))+Traz*log(R2+(Xm1-Sta(i,1)))) ... -(Trax*log(R3+(Zm2-Sta(i,3)))-Tray*atan((Xm2-Sta(i,1))*(Ym1-Sta(i,2))/((Ym1-Sta(i,2))^2+R3*(Zm2-Sta(i,3))+(Zm2-Sta(i,3))^2))+Traz*log(R3+(Xm2-Sta(i,1)))) ... +(Trax*log(R4+(Zm2-Sta(i,3)))-Tray*atan((Xm1-
评论
    相关推荐