《应用地球物理学》课程作业
基于MATLAB的球体重力正演程序实验报告
1
一 程序简介
本程序基于MATLAB软件的GUI模块编写，旨在实现球体重力正演结果的可视化分析。MATLAB是一个高级的编程语言，其矩阵思想方便了地球物理的编程工作。随着该语言和相应软件的发展，其内部也集成了许多模块，如该实验用到的GUI模块。在该模块中，可以通过窗口、按键和赋值框等基本元素的组合，编写出可视化的应用程序，再配合MATLAB强大的作图功能，可以实现正演结果的展示与分析。
该程序应包含以下内容：
1. 可以自由输入参数，如球体半径，埋深和剩余密度。 2. 可以计算出Δg、VZZ、VXZ和VZZZ这四种重力异常及其导数的对应值。
3. 可以绘制剖面图及平面图两种图像。
二 源程序
由于GUI程序的头文件均大同小异，这里只列出赋值框及绘图按键的程序代码。
% --- Executes on selection change in popupmenu1. function popupmenu1_Callback(hObject, eventdata, handles) % hObject handle to popupmenu1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)
% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu1
contents as cell array
% contents{get(hObject,'Value')} returns selected item from popupmenu1
s=get(hObject,'value'); handles.s = s;
guidata(hObject, handles);
function edit1_Callback(hObject, eventdata, handles) % hObject handle to edit1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of edit1 as text % str2double(get(hObject,'String')) returns contents of edit1 as a double
d = str2double(get(hObject,'string')); handles.d = d;
guidata(hObject, handles);
function edit2_Callback(hObject, eventdata, handles) % hObject handle to edit2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of edit2 as text % str2double(get(hObject,'String')) returns contents of edit2 as a double
r = str2double(get(hObject,'string')); handles.r = r;
guidata(hObject, handles);
function edit3_Callback(hObject, eventdata, handles) % hObject handle to edit3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of edit3 as text % str2double(get(hObject,'String')) returns contents of edit3 as a double
ro = str2double(get(hObject,'string')); handles.ro = ro;
guidata(hObject, handles);
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)
ro = handles.ro;r = handles.r;d = handles.d; x=-2*d:2*d;
G=6.67e-11;pi=3.14159; switch handles.s case 2
z=4*ro*pi*G*r^3*d./(3*(x.^2+d^2).^1.5);
plot(x,z.*1e6);xlabel('X/m');ylabel('\\Deltag/g.u.'); case 3
z=4*ro*pi*G*r^3*(2*d^2-x.^2)./(3*(x.^2+d^2).^2.5); plot(x,z.*1e9);xlabel('X/m');ylabel('Vzz/E'); case 4
z=-4*ro*pi*G*r^3*(d.*x)./(x.^2+d^2).^2.5; plot(x,z.*1e9);xlabel('X/m');ylabel('Vxz/E'); case 5
z=4*ro*pi*G*r^3*(2*d^2-3.*x.^2)./(x.^2+d^2).^3.5; plot(x,z.*1e9);xlabel('X/m');ylabel('Vzzz/nMKS'); end
% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles) % hObject handle to pushbutton3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) ro = handles.ro;r = handles.r;d = handles.d; x=-2*d:2*d;y=-2*d:2*d; [X,Y]=meshgrid(x,y); G=6.67e-11;pi=3.14159; switch handles.s case 2
c=4*ro*pi*G*r^3*d./(3*(X.^2+Y.^2+d^2).^1.5);
contour(X,Y,c.*1e6,'showtext','on');colormap(winter);xlabel('X/m');ylabel('Y/m'); case 3
c=4*ro*pi*G*r^3*(2*d^2-X.^2-Y.^2)./(3*(X.^2+Y.^2+d^2).^2.5);
contour(X,Y,c.*1e9,'showtext','on');colormap(winter);xlabel('X/m');ylabel('Y/m');
case 4
c=-4*ro*pi*G*r^3*(d.*X.*Y)./(X.^2+Y.^2+d^2).^2.5;
contour(X,Y,c.*1e9,'showtext','on');colormap(winter);xlabel('X/m');ylabel('Y/m'); case 5
c=4*ro*pi*G*r^3*(2*d^2-3.*X.^2-Y.^2)./(X.^2+Y.^2+d^2).^3.5;
contour(X,Y,c.*1e9,'showtext','on');colormap(winter);xlabel('X/m');ylabel('Y/m'); end
三 程序运行结果
假设埋深为100m，球体半径为10m，剩余密度为2kg/m3. Δg的剖面图和平面图如下：
Vzz的剖面图和平面图如下：