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")

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;  %对空值填0actnums = 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"), colorbarend

最后结果如下：