MRST绘制三维网格图
MRST绘制三维网格图
Matlab储层模拟工具箱(The MATLAB Reservoir Simulation Toolbox, MRST)是一款用于储层建模的免费、开源的软件，主要由 SINTEF Digital 数学与控制论系的计算地球科学小组开发。
更多介绍MRST (sintef.no)
在下载完工具箱之后，将其加入库路径，随后开始绘图。
首先是针对规则网格的绘制与展示：
1 规则网格
Step1 用cartGrid可以创建一个空的笛卡尔网格，如下图(a)所示；
G = cartGrid([3, 3, 3], [1, 1, 1].*meter); %创建一个3*3*3的网格
plotGrid(G), view([-40, 32]), xlabel('x'), ylabel('y'), zlabel('Depth')
Step2 对空网格进行赋值；
在这一步中我们对这3×3×3=27个网格进行赋值，然后成图，如图中(b)
G = cartGrid([3, 3, 3], [1, 1, 1].*meter);
Values = linspace(1,27,27).'; %将数据转成与对应网格数相同的n×1的序列即可
plotCellData(G, Values) %通过plotCellData显示
view([-40, 32]), xlabel('x'), ylabel('y'), zlabel('Depth'), colorbar
如果我们需要的网格就是这类规则的方格，那成图过程简单概括也就是“建网格，贴数值”。
2 不规则网格
Step1 用cartGrid可以创建一个空的笛卡尔网格，抠除掉一个无效的网格，使之成为透明的，如下图(a)所示；
%同样还是先创建一个3*3*3的网格
% 创建 3x3x3 的笛卡尔网格
G = cartGrid([3, 3, 3]);
actnum = ones(G.cells.num, 1);
actnum(21) = 0;
% 移除不可见的网格
G = removeCells(G, actnum == 0);
plotGrid(G), view([-40, 32]), xlabel('x'), ylabel('y'), zlabel('Depth')
Step2 对空网格进行赋值；
%% 抠网格贴颜色
%同样还是先创建一个3*3*3的网格
% 创建 3x3x3 的笛卡尔网格
G = cartGrid([3, 3, 3]);
actnum = ones(G.cells.num, 1);
actnum(21) = 0;
% 移除不可见的网格
G = removeCells(G, actnum == 0);
% 为剩余网格生成随机数值
values = rand(G.cells.num, 1);
plotCellData(G, values); %通过plotCellData显示
view([-40, 32]), xlabel('x'), ylabel('y'), zlabel('Depth'), colorbar
将得到这样一个结果：
3 任意网格
当我们需要对一个任意的三维数据，进行3D网格展示时，执行的操作可概括为，创建一个范围足够的规则网格，随后从规则的网格里抠除无值部分。具体到如切片数据时，仅需依据Line、CDP、Time范围对三维网格赋值、重排，坐标排序格式如下，包含X,Y,Z和指定位置处的Value，有值时Actnum为1，否则为0：
完整代码如下：
filename = 'Demo.slice'; %切片文件
[line,cdp,X,Y,Time,value]=readslicedata(filename); %读取切片函数
showslice3D(line,cdp,Time,value); %切片三维可视化
其中包含两个子函数
（1）读取切片函数readslicedata
function [line,cdp,X,Y,Time,value]=readslicedata(filename)
fid=fopen(filename,'rt');
if fid<0
warndlg('打开文件失败!');
return;
end
FormatString=repmat('%f',1,7);
dataoutput=textscan(fid,FormatString,'HeaderLines',1);
line=dataoutput{2};
cdp=dataoutput{3};
X=dataoutput{4};
Y=dataoutput{5};
Time=dataoutput{6};
value=dataoutput{7};
fclose(fid);
end
（2）切片三维可视化函数showslice3D
function showslice3D(line,cdp,Time,value)
%汇总切片信息
num = size(line,1);
sliceInfo = zeros(num,5);
sliceInfo(:,1)=line;
sliceInfo(:,2)=cdp;
sliceInfo(:,3)=int64(Time);
sliceInfo(:,4)=value;
sliceInfo(:,5)=1;
%获取线道号、时间最大、最小范围
minLine=min(line);
maxLine=max(line);
minCDP=min(cdp);
maxCDP=max(cdp);
minTime=min(Time);
maxTime=max(Time);
%重组切片坐标
[X, Y, Z] = meshgrid(minLine:maxLine, minCDP:maxCDP, minTime:maxTime);
COO = [X(:), Y(:), Z(:)];
% 数组转表格
COOtable = array2table(COO, 'VariableNames', {'X', 'Y', 'Z'});
Slicetable = array2table(sliceInfo, 'VariableNames', {'X', 'Y', 'Z', 'Value','Valid'});
% 连接两表
joinedTable = outerjoin(COOtable, Slicetable,'Keys', {'X', 'Y', 'Z'});
selectedTable = joinedTable(:, {'X_COOtable', 'Y_COOtable', 'Z_COOtable', 'Value','Valid'});
selectedTable = sortrows(selectedTable, [3, 2, 1]); %按z->y->x的顺序重新排序
data=table2array(selectedTable);
values = data(~isnan(data(:,4)),4); %只取有值的部分
data(isnan(data)) = 0; %对空值填0
actnums = data(:,5);
Xnum = maxLine-minLine+1;
Ynum = maxCDP-minCDP+1;
Znum = maxTime-minTime+1;
G = cartGrid([Xnum, Ynum, Znum]);
% 移除不可见的网格
G = removeCells(G, actnums == 0);
plotCellData(G, values,'EdgeColor', 'none'); %通过plotCellData显示
view([-40, 32]), xlabel('Line'), ylabel('CDP'), zlabel('Time'), colorbar
end
最后结果如下：
