本文将地质体目标离散成一些矩形二度体的组合，对于复杂情况，只需定义密度异常非零的二度体单元，这样会大大节省内存。
1数据结构
1.1定义以下的数据结构用于存储模型参数
{
模型中矩形二度体的个数
矩形二度体的水平位置
矩形二度体的垂直位置
矩形二度体的宽
矩形二度体的高
矩形二度体的密度
}
matlab代码：
function
[Model] =
modelDef(Size,XArray,ZArray,dXArray,dZArray,RhoArray)
Model.Size=Size;
Model.X=XArray;
Model.Z=ZArray;
Model.Dx=dXArray;
Model.Dz=dZArray;
Model.Rho=RhoArray;
end
1.2定义如下结构体用于描述观测系统
{
测线的观测点数
测点的水平位置
测点的垂直位置
测点的重力值
重力值的单位
}
matlab代码：
function
[ObsSystem]=ObsSystemDef(Size,XArray,ZArray,VArray,Unit)
ObsSystem.Size=Size;
ObsSystem.X=XArray;
ObsSystem.Z=ZArray;
ObsSystem.Value=VArray;
ObsSystem.Unit=Unit;
end
2重力正演计算
由地质模型计算观测系统上的重力异常
function
ObsSystem=Model2ObsSys(model,ObsSystem)
G=6.67e-11*1e5;
ObsSystem.Value(1:ObsSystem.Size)=0;
for
(j=1:ObsSystem.Size)
for
(i=1:model.Size)
x1=model.X(i)-model.Dx(i)/2-ObsSystem.X(j);
x2=model.X(i)+model.Dx(i)/2-ObsSystem.X(j);
z1=model.Z(i)-model.Dz(i)/2-ObsSystem.Z(j);
z2=model.Z(i)+model.Dz(i)/2-ObsSystem.Z(j);
if
(model.Rho(i)~=0)
gg=x1*log(x1^2+z1^2)+2*z1*atan(x1/z1)...
-x1*log(x1^2+z2^2)-2*z2*atan(x1/z2)...
-x2*log(x2^2+z1^2)-2*z1*atan(x2/z1)...
+x2*log(x2^2+z2^2)+2*z2*atan(x2/z2);
ObsSystem.Value(j)=ObsSystem.Value(j)+gg*model.Rho(i);
end
end
end
ObsSystem.Value=ObsSystem.Value*G;
ObsSystem.Unit='mgal';
end
matlab绘图程序：
function
plotModel(model)
hold
on;
modelXmin=min(model.X)-max(model.Dx);
modelXmax=max(model.X)+max(model.Dx);
modelZmin=min(model.Z)-max(model.Dz);
modelZmax=max(model.Z)+max(model.Dz);
modelW=modelXmax-modelXmin;
modelH=modelZmax-modelZmin;
modelX=modelXmin;
modelZ=modelZmin;
modelRmin=min(model.Rho);
modelR=max(model.Rho)-min(model.Rho);
rectangle('Position',[modelX,modelZ,modelW,modelH],...
'LineStyle','--');
for(i=1:model.Size)
if
(model.Rho(i)~=0)
rectangle('Position',[model.X(i)-model.Dx(i)/2,model.Z(i)-model.Dz(i)/2,model.Dx(i),model.Dz(i)],...
'FaceColor',[(model.Rho(i)-modelRmin)/modelR
1-(model.Rho(i)-modelRmin)/modelR 1]);
end
end
xlim([modelXmin modelXmax]);
ylim([modelZmin modelZmax]);
xlabel('X
(m)');
ylabel('Z
(m)');
set(gca,'ydir','reverse');
hold
off;
end
function
plotModelJet16(model,cmin,cmax)
hold
on;
modelXmin=min(model.X)-max(model.Dx);
modelXmax=max(model.X)+max(model.Dx);
modelZmin=min(model.Z)-max(model.Dz);
modelZmax=max(model.Z)+max(model.Dz);
modelW=modelXmax-modelXmin;
modelH=modelZmax-modelZmin;
modelX=modelXmin;
modelZ=modelZmin;
modelRmin=min(model.Rho);
modelR=max(model.Rho)-min(model.Rho);
rectangle('Position',[modelX,modelZ,modelW,modelH],...
'LineStyle','--');
colorindex=jet(16);
colormap(jet(16));
rho=model.Rho;
rho=(rho-cmin)/(cmax-cmin)*16;
rho(rho<1)=1;
rho(rho>16)=16;
rho=round(rho);
for(i=1:model.Size)
if
(model.Rho(i)~=0)
rectangle('Position',[model.X(i)-model.Dx(i)/2,model.Z(i)-model.Dz(i)/2,model.Dx(i),model.Dz(i)],...
'FaceColor',colorindex(rho(i),:));
end
end
�xis([cmin cmax]);
%colorbar;
xlim([modelXmin modelXmax]);
ylim([modelZmin modelZmax]);
xlabel('X
(m)');
ylabel('Z
(m)');
set(gca,'ydir','reverse');
hold
off;
end
function
plotObsSystem(ObsSystem)
hold
on;
Xmin=min(ObsSystem.X);
Xmax=max(ObsSystem.X);
Vmin=min(ObsSystem.Value);
Vmax=max(ObsSystem.Value);
scatter(ObsSystem.X,ObsSystem.Value);
xlim([Xmin Xmax]);
if
(Vmin~=Vmax)
ylim([Vmin Vmax]);
end
xlabel('X
(m)');
ylabel('gravity
anomaly (mgal)');
hold
off;
end
3例子二度圆柱体正演
clear
all;
clf;
%定义二度圆柱体
lx=100;
dx=2e3;
lz=100;
dz=2e3;
xx=(0:lx-1)*dx;
zz=(1:lz)*dz;
[xx,zz]=meshgrid(xx,zz);
%均匀密度圆柱体
drho_model=zeros(lx,lz);
%半径
radius=22.5e3;
%埋深
depth=30e3;
%圆心距离原点的水平距离
x0=45e3;
%密度差
drho=1000;
%生成密度模型
drho_model(((xx-x0).^2+(zz-depth).^2)
%将地质模型转换为结构体
%也可以将密度非零的二度体挑出来定义，将会节省内存，本文仅仅简单转换，包含了冗余的二度体
Size=lx*lz;
XArray=reshape(xx,1,[]);
ZArray=reshape(zz,1,[]);
dXArray=XArray*0+dx;
dZArray=ZArray*0+dz;
RhoArray=reshape(drho_model,1,[]);
model=modelDef(Size,XArray,ZArray,dXArray,dZArray,RhoArray);
%绘制模型
%subplot(2,2,2);
%plotModel(model);
%定义观测系统
Size=lx;
XArray=(0:lx-1)*dx;
ZArray=(0:lx-1)*0;
VArray=(0:lx-1)*0;
Unit='mgal';
ObsSystem=ObsSystemDef(Size,XArray,ZArray,VArray,Unit);
%正演计算圆柱体模型的重力异常值
ObsSystem=Model2ObsSys(model,ObsSystem);
%绘出重力异常值
%subplot(2,2,1);
%plotObsSystem(ObsSystem);