526互联

MRST绘制三维网格图

发布时间 2023-06-06 00:09:37作者: GeoFXR

MRST绘制三维网格图

Matlab储层模拟工具箱(The MATLAB Reservoir Simulation Toolbox, MRST)是一款用于储层建模的免费、开源的软件，主要由 SINTEF Digital 数学与控制论系的计算地球科学小组开发。

更多介绍MRST (sintef.no)

下载地址Download (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')
image-20220531095049587

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

将得到这样一个结果：

image-20230605234517809

3 任意网格

当我们需要对一个任意的三维数据，进行3D网格展示时，执行的操作可概括为，创建一个范围足够的规则网格，随后从规则的网格里抠除无值部分。具体到如切片数据时，仅需依据Line、CDP、Time范围对三维网格赋值、重排，坐标排序格式如下，包含X,Y,Z和指定位置处的Value，有值时Actnum为1，否则为0：

image-20220531133416676

完整代码如下：

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

最后结果如下：

切片2