下面给出一个思路示例，展示如何在 MATLAB 中基于“椭圆掩码 + 二维正态分布”来模拟楼梯在长时间踩踏下的磨损过程。思路分为以下几个步骤：

---

1. 建立楼梯几何及网格

假设楼梯是一个矩形立方体，尺寸为：

• x 方向长度 1 m
• y 方向宽度 0.4 m
• 高度（z 方向）0.3 m

在数值模拟中，我们往往只关注“楼梯表面的磨损”即 z 方向高度 $h(x,y)$ 随时间的演化，因此可以在平面 $(x,y)$ 上离散出一个网格，用一个矩阵来表示每个网格点的当前高度。

示例：在 x ∈ [0,1]，y ∈ [0,0.4] 的区域做均匀网格划分

matlab
% 1. 建立网格
Nx = 101;  % x方向网格数
Ny = 41;   % y方向网格数
xVals = linspace(0, 1, Nx);
yVals = linspace(0, 0.4, Ny);
[XX, YY] = meshgrid(xVals, yVals);

% 2. 定义楼梯表面高度 h(x,y)，初始均为 0.3 m
h = 0.3 * ones(Ny, Nx);

这样就得到一个大小为 (Ny x Nx) 的矩阵 h 表示楼梯表面的高度。

---

2. 定义踩踏分布（概率分布 + 椭圆掩码）

2.1 二维正态分布

用户给出“踩踏点”在楼梯上的分布服从二维正态分布，均值与标准差分别为：

• 均值 $\mu_x = 0.5$，$\mu_y = 0.15$
• 标准差 $\sigma_x = 0.3$，$\sigma_y = 0.2$

可根据需要做截断，或者直接在楼梯区域外的概率视为 0。

在 MATLAB 中，可直接用 mvnrnd 函数一次性生成随机样本点，或者自己编写二维正态函数并采样：

matlab
Nwalk = 1e4;   % 生成 1万 个踩踏点
mu    = [0.5, 0.15];
Sigma = [0.3^2, 0; 0, 0.2^2];  % 协方差矩阵：对角为 sigma_x^2, sigma_y^2
pts   = mvnrnd(mu, Sigma, Nwalk);

% 只取落在楼梯范围 0<=x<=1, 0<=y<=0.4 的点
validMask = (pts(:,1)>=0 & pts(:,1)<=1 & pts(:,2)>=0 & pts(:,2)<=0.4);
pts = pts(validMask, :);

经过这一步，就得到了分布在楼梯表面的若干“踩踏点”。

2.2 椭圆掩码

假设每次踩踏在楼梯表面形成的接触区域近似为椭圆，给定：

matlab
ellipseRx = 0.02;  % 椭圆 x方向半径
ellipseRy = 0.01;  % 椭圆 y方向半径

如果有一点 $(x_0,y_0)$ 是踩踏中心，则它对应的“接触区域”可以表示为满足

$$\frac{(x-x_0)^2}{\text{ellipseRx}^2} \;+\; \frac{(y-y_0)^2}{\text{ellipseRy}^2} \;\le\; 1$$

的点的集合。

---

3. 定义磨损模型（离散形式）

3.1 磨损速率方程

你给出的连续形式为

$$\frac{\partial h}{\partial t}(x,y,t) = -\,\alpha_{\text{bulk}} \int_{\Omega} K_{\text{bulk}}(x-\xi, y-\eta) \,P_{\text{bulk}}(\xi,\eta,t) \,\mathrm{d}\xi\,\mathrm{d}\eta \;-\; \alpha_{\text{edge}} \int_{\Omega} K_{\text{edge}}(x-\xi, y-\eta) \,P_{\text{edge}}(\xi,\eta,t) \,\mathrm{d}\xi\,\mathrm{d}\eta.$$

在简化情形下，可以把每次踩踏（出现一次接触区域）看作对楼梯表面某些网格点做高度扣减（或增加“损耗值”）。若不细分 bulk 与 edge 两种情形，先做最直接的 离散化：

1. 每出现一脚，产生一个椭圆区域。
2. 对该区域内的网格点，按照某个常数或者空间核函数 $K(\Delta x,\Delta y)$ 进行磨损扣减。

这样我们就不必直接实现上面带有积分卷积形式的方程，而是用蒙特卡洛式的“每脚一扣减”来逼近。后面若需要区分 bulk/edge，可在接触区域外围加一圈“edge”核函数。

---

3.2 用蒙特卡洛迭代实现

假设一天有 $N_{\text{day}}$ 次踩踏（或者统一按总次数），每出现一次随机踩踏中心 $(x_0, y_0)$，对落在椭圆范围内的网格点高度进行扣减。示例：

matlab
% 给定每脚对椭圆内网格点的磨损深度(可视为常数或包含某种分布)
wearDepthPerStep = 5e-5;  % 例如一次踩踏的局部磨损深度

for i = 1 : size(pts,1)
    x0 = pts(i,1);
    y0 = pts(i,2);

    % 找出网格中哪些点落入椭圆掩码
    % 椭圆方程 ( (x - x0)/Rx )^2 + ( (y - y0)/Ry )^2 <= 1
    distEllipse = ((XX - x0)/ellipseRx).^2 + ((YY - y0)/ellipseRy).^2;
    inMask = distEllipse <= 1;

    % 对掩码区域对应的高度做磨损扣减
    h(inMask) = h(inMask) - wearDepthPerStep;
end

注：上面是一种最简单的离散蒙特卡洛思路，每一脚都“确定性扣减”同样深度。若你需要引入更复杂的材料学或者核函数 $K_{\text{bulk}}, K_{\text{edge}}$ 等，也可以在每次扣减时使用不同的加权系数。

---

4. 考虑核函数 (Bulk/Edge) 的实现

如果要更贴近你的公式 (B)，即：

$$\int K_{\text{bulk}}\!(x-\xi,y-\eta)\,P_{\text{bulk}}(\xi,\eta,t)\,\mathrm{d}\xi\,\mathrm{d}\eta,$$

可以将该核函数离散化为一个“卷积核”矩阵，比如在 x-y 网格上定义一个矩阵 $K_{\text{bulk}}[\Delta i,\Delta j]$。当某处出现一次踏点，就相当于在那个踏点中心上，对整张楼梯表面的“磨损场”加上一张卷积核的“叠加”。代码示例如下：

matlab
% 1) 先构造K矩阵, 假设空间范围在 [-3*Rx, 3*Rx] x [-3*Ry, 3*Ry]
RxGrid = 3 * ellipseRx; 
RyGrid = 3 * ellipseRy;
dx = xVals(2) - xVals(1);
dy = yVals(2) - yVals(1);
ixRange = -floor(RxGrid/dx) : floor(RxGrid/dx);
iyRange = -floor(RyGrid/dy) : floor(RyGrid/dy);

[KX, KY] = meshgrid(ixRange*dx, iyRange*dy);

% 一个示例：K_bulk 为 2D 高斯核或别的函数
sigmaKx = ellipseRx / 2;
sigmaKy = ellipseRy / 2;
K_bulk = exp(-((KX).^2/(2*sigmaKx^2) + (KY).^2/(2*sigmaKy^2)));
K_bulk = K_bulk / sum(K_bulk(:));  % 归一化

% 2) 在每个踏点 (x0,y0) 的网格索引附近，对 h 做叠加扣减
for i = 1 : size(pts,1)
    x0 = pts(i,1);
    y0 = pts(i,2);
    % 找到 (x0,y0) 对应到网格的索引中心 (ix0, iy0)
    [~, ix0] = min(abs(xVals - x0));
    [~, iy0] = min(abs(yVals - y0));

    % 计算本次K_bulk应加到 h 的哪一片区
    % 这里需要注意边界检查
    for ii = 1:length(ixRange)
        for jj = 1:length(iyRange)
            gx = ix0 + ixRange(ii);
            gy = iy0 + iyRange(jj);
            if gx >=1 && gx <= Nx && gy >=1 && gy <= Ny
                h(gy,gx) = h(gy,gx) - wearDepthPerStep * K_bulk(jj, ii);
            end
        end
    end
end

其中：

• K_bulk 只是一个示例，你可以根据材料学、脚掌形状、速度、重力等细节自定义更合理的核函数。
• 如果要区分 K_edge 与 K_bulk，只需额外再做一个 K_edge 并在某些区域或条件下叠加相应磨损。

---

5. 观察和可视化结果

在迭代完成后，我们可用 surf、imagesc 等函数查看楼梯表面高度分布的演变。例如：

matlab
figure; 
surf(xVals, yVals, h); 
shading interp; colorbar;
title('楼梯磨损后的表面高度'); 
xlabel('x'); ylabel('y'); zlabel('h(x,y)');

若需要做动画，可在每次扣减后或每隔若干步后更新图像。

---

6. 总结与提示

1. 宏观方程 (B) 与离散实现
   你列出的方程 (B) 可以视为一种“空间卷积 + 时间累积”形式；在数值上，通常可以借助卷积(convolution) 或者蒙特卡洛随机叠加来实现。
2. 随机分布 & 椭圆掩码
   • 先用二维正态分布产生随机踩踏点。
   • 再用椭圆掩码或核函数对楼梯表面做局部扣减，逼近实际的磨损过程。
3. 区分 Bulk / Edge
   如果要把“台阶边缘处的额外磨损”提出来，可以令 $K_{\text{edge}}$ 对应楼梯边缘处出现时，磨损更大；也可以单独统计踏点离楼梯边缘的距离等信息，通过额外系数 $\alpha_{\text{edge}}$ 来加大磨损强度。
4. 单位和尺度
   注意检查单位 (m, s, 次/天, 磨损深度量级等)。通常在实际模拟中，一次踩踏的磨损可能非常小，需要很多次叠加才显现显著的高度变化。
5. 加速与边界
   如果网格过细、踩踏次数很多，直接在循环中做卷积可能效率较低，可以考虑用 FFT 做 2D 卷积加速（MATLAB 中的 conv2 或 fft2 等），或者用更高效的数据结构处理。

通过以上步骤，就可以在 MATLAB 里初步搭建一个基于“椭圆掩码 + 二维正态分布”来模拟楼梯踩踏磨损的数值演示模型。你也可以不断迭代、细化材料属性和分布假设，使之更贴近实际的物理情形。

下面给出一个自包含的 MATLAB 脚本示例，演示如何利用「二维正态分布 + 椭圆掩码」来模拟楼梯表面因踩踏而产生的磨损。所有可调整的参数都放在脚本前面，便于修改和调试。请将以下内容复制到一个名为 simulateStairWear.m 的文件中，然后在 MATLAB 中运行。

---

matlab
function simulateStairWear()
% SIMULATESTAIRWEAR - 基于“椭圆掩码 + 二维正态分布”模拟楼梯磨损
%
%   作者示例：完整代码 + 注释
%
%----------------------------------------------------------------------
%% =============== 0. 参数定义（可以自由调整） =========================
%
% 楼梯尺寸和网格相关
stairLength    = 1.0;    % 楼梯 x 方向长度 (米)
stairWidth     = 0.4;    % 楼梯 y 方向宽度 (米)
initialHeight  = 0.3;    % 初始楼梯高度 (米)
Nx             = 101;    % x方向网格数
Ny             = 41;     % y方向网格数
%
% 踩踏点(人流分布)相关
Nwalk          = 1e4;          % 生成多少个踩踏点(总人次)
muX            = 0.5;          % 踩踏点在 x 方向的均值
muY            = 0.15;         % 踩踏点在 y 方向的均值
sigmaX         = 0.3;          % x方向标准差
sigmaY         = 0.2;          % y方向标准差
%
% 每个踩踏的“椭圆掩码”形状
ellipseRx      = 0.02;         % 椭圆 x方向半径 (米)
ellipseRy      = 0.01;         % 椭圆 y方向半径 (米)
%
% 磨损深度相关
wearDepthPerStep = 5e-5;       % 每次踩踏对楼梯产生的磨损深度(米)
%
%----------------------------------------------------------------------
%% =============== 1. 建立楼梯表面网格 =============================
% 在 x ∈ [0, stairLength], y ∈ [0, stairWidth] 上做均匀网格

xVals = linspace(0, stairLength, Nx);  % x方向坐标
yVals = linspace(0, stairWidth,  Ny);  % y方向坐标
[XX, YY] = meshgrid(xVals, yVals);     % 生成二维网格
%
% 初始化楼梯表面高度矩阵 h，大小为 Ny x Nx
% (注：MATLAB 中默认 h(row, col)，所以 y 对应行索引, x 对应列索引)
h = initialHeight * ones(Ny, Nx);  

%----------------------------------------------------------------------
%% =============== 2. 生成踩踏点 (二维正态随机采样) ================
% 假设踩踏点 (x0, y0) ~ N( [muX, muY], [sigmaX^2, 0; 0, sigmaY^2] )

mu    = [muX, muY];
Sigma = [sigmaX^2, 0; 0, sigmaY^2];  % 协方差矩阵
pts   = mvnrnd(mu, Sigma, Nwalk);    % 随机产生 Nwalk 个点 (x,y)

% 只保留落在楼梯顶面范围 [0, stairLength] x [0, stairWidth] 的点
validMask = (pts(:,1) >= 0 & pts(:,1) <= stairLength & ...
             pts(:,2) >= 0 & pts(:,2) <= stairWidth);
pts = pts(validMask, :);

%----------------------------------------------------------------------
%% =============== 3. 蒙特卡洛式迭代：每脚产生一次磨损 =============
%
% 思路：逐一读取随机踩踏点 (x0, y0)，构造椭圆掩码 
%       对 h(x,y) 进行磨损扣减
%
nPoints = size(pts,1);  % 有效踩踏点总数

for i = 1 : nPoints
    
    % 读取本次踩踏点中心 (x0, y0)
    x0 = pts(i,1);
    y0 = pts(i,2);
    
    % 计算网格中哪些点落入椭圆掩码
    % 椭圆方程： ((x - x0)/ellipseRx)^2 + ((y - y0)/ellipseRy)^2 <= 1
    distEllipse = ((XX - x0)/ellipseRx).^2 + ((YY - y0)/ellipseRy).^2;
    inMask = (distEllipse <= 1);
    
    % 对掩码范围内的网格点做磨损：减去 wearDepthPerStep
    h(inMask) = h(inMask) - wearDepthPerStep;
    
end

% 在长时间(多次脚步)累加后，楼梯表面的高度分布 h 已经变化

%----------------------------------------------------------------------
%% =============== 4. 可视化楼梯表面高度 ===========================
%
% 使用 surf 函数查看最终结果，也可用 imagesc 等

figure('Name','楼梯磨损模拟结果','Color','w');

% 注意 surf 默认 x->列索引, y->行索引：这里需将 xVals 与 yVals 顺序传入
surf(xVals, yVals, h, 'EdgeColor','none');
xlabel('x (m)'); ylabel('y (m)'); zlabel('Height (m)');
title('磨损后的楼梯表面高度分布');
axis tight; colorbar; shading interp;
view(45,45); % 旋转一下视角

%----------------------------------------------------------------------
%% =============== 5. 统计与显示分布信息(可选) ======================
%
% 例如，查看最终高度的平均值或最小值:
minH = min(h(:));
meanH = mean(h(:));
disp('---------------------------------------------');
disp(['最终楼梯最小高度: ', num2str(minH), ' m']);
disp(['最终楼梯平均高度: ', num2str(meanH), ' m']);
disp(['踩踏次数(有效点数): ', num2str(nPoints)]);

end

---

代码说明

1. 可调参数
   代码开头部分（%% =============== 0. 参数定义 =========================）列出所有可调参数，例如网格大小、楼梯尺寸、二维正态分布参数（$\mu$ 与 $\sigma$）、踩踏次数、椭圆半径、磨损深度等。

2. 楼梯表面网格 (h)

   • 在 $[0, \text{stairLength}] \times [0, \text{stairWidth}]$ 的 2D 平面内，构造大小为 Ny x Nx 的网格并初始化高度 h = initialHeight。
   • 变量 XX, YY 分别存储网格点的 x, y 坐标。

3. 踩踏点分布

   • 用 mvnrnd 函数生成二维正态分布的随机采样点 (x, y)。
   • 参数 muX, muY, sigmaX, sigmaY 控制分布中心与离散程度。
   • 用 validMask 筛掉落在楼梯区域外的点，保留最终分布。

4. 蒙特卡洛式磨损叠加

   • 对每一个踩踏点 $(x_0, y_0)$，构建椭圆掩码： $$\frac{(x - x_0)^2}{\text{ellipseRx}^2} + \frac{(y - y_0)^2}{\text{ellipseRy}^2} \le 1.$$
   • 找到满足该条件的网格点后，对应 h 矩阵做一次固定深度（wearDepthPerStep）的扣减。
   • 这样经过所有踩踏点累加，就得到最后磨损形态。

5. 可视化

   • 通过 surf 命令显示结果的 3D 视图，也可用 imagesc+contour 等方式查看。
   • 最后打印出最小高度 minH、平均高度 meanH 等简单统计信息。

如果想要模拟“核函数卷积形式”($K_{\text{bulk}}$或$K_{\text{edge}}$)，可以在 step 3 中改为根据距离进行不同权重的扣减，或者利用 conv2 函数做基于卷积的快速叠加。

---

运行该脚本后，即可获得一个简单的楼梯磨损数值模拟示例。若需要更高精度或更长时间模拟，可以加大网格密度或增大踩踏次数；如果关心真实物理单位，就需要根据实际人流、材料、时间尺度来调整 wearDepthPerStep 和 Nwalk 等参数。

下面给出一个完整可运行的 MATLAB 脚本示例，演示如何在“以千年为单位”的迭代过程中，实时刷新并生成三个图像：

1. 平面热力图（2D 颜色图）展示楼梯高度分布。
2. 立体图（3D surf）可旋转查看楼梯表面的磨损形态。
3. 掩码生成图：把每次踩踏对应的椭圆区域（掩码）都画在同一张 2D 图上，随着迭代进行不断累加展示。

为演示方便，脚本中将“总踩踏次数”拆分为“每千年脚步数”与“千年数”相乘，从而实现“以千年为单位的迭代”。在实际应用中，可根据真实人流量、时间尺度等进行调整。请将以下代码保存为 simulateStairWear_3Figures.m，在 MATLAB 中运行。

---

matlab
function simulateStairWear_3Figures()
% SIMULATESTAIRWEAR_3FIGURES - 以千年为单位进行楼梯磨损模拟，并在迭代过程中
%                               动态绘制三个图：
%                               1) 平面热力图
%                               2) 三维立体图
%                               3) 踩踏椭圆掩码分布图（所有踩踏都叠加）
%
% 说明：
%  - 每个“脚步”都从预先生成的二维正态分布中抽样出 (x,y)，
%  - 并在楼梯表面用椭圆掩码做磨损扣减。
%  - 将脚步数拆分为 [每千年脚步数] x [千年数]，并在迭代中每隔 n 步刷新绘图。
%
% 作者示例代码，欢迎根据需要自行修改参数。

    %% =========== 0. 可调参数：楼梯尺寸 & 网格 & 时间 & 分布等 ===========
    % 楼梯尺寸和网格相关
    stairLength      = 1.0;    % 楼梯 x 方向长度 (米)
    stairWidth       = 0.4;    % 楼梯 y 方向宽度 (米)
    initialHeight    = 0.3;    % 初始楼梯高度 (米)
    Nx               = 101;    % x方向网格数
    Ny               = 41;     % y方向网格数

    % 时间与脚步：以“千年”为单位的迭代次数
    totalMillennia       = 2;         % 总共模拟多少个“千年”
    footStepsPerMillenni = 500;       % 每一千年有多少“脚步”
    totalSteps           = totalMillennia * footStepsPerMillenni;
    
    % 每隔多少步刷新/绘制一次图像（避免每步都画图过慢）
    plotInterval = 100;              

    % 二维正态分布参数
    muX            = 0.5;    % x方向均值
    muY            = 0.15;   % y方向均值
    sigmaX         = 0.3;    % x方向标准差
    sigmaY         = 0.2;    % y方向标准差

    % 椭圆掩码参数
    ellipseRx      = 0.02;   % 椭圆 x方向半径 (米)
    ellipseRy      = 0.01;   % 椭圆 y方向半径 (米)

    % 每次踩踏对楼梯的磨损深度
    wearDepthPerStep = 5e-5; % (米)

    % 为了在第三张图上绘制所有椭圆边界，需要采样多少个点来画椭圆
    ellipsePlotResolution = 36;  % 较大一些，可使椭圆平滑

    %% =========== 1. 建立楼梯表面网格并初始化高度 ====================
    xVals = linspace(0, stairLength, Nx);  % x方向坐标序列
    yVals = linspace(0, stairWidth,  Ny);  % y方向坐标序列
    [XX, YY] = meshgrid(xVals, yVals);     % 生成二维网格 (Ny x Nx)
    h = initialHeight * ones(Ny, Nx);      % 楼梯表面高度矩阵

    %% =========== 2. 生成所有脚步的随机分布 (二维正态) ================
    % 总共需要 totalSteps 个踩踏点
    mu    = [muX, muY];
    Sigma = [sigmaX^2, 0; 0, sigmaY^2];  
    pts   = mvnrnd(mu, Sigma, totalSteps);  % (totalSteps x 2) 的矩阵

    % 只保留落在 [0, stairLength] x [0, stairWidth] 范围内的点，若超出则重新采样
    % 为简单起见，这里直接用一个循环来保证生成够 totalSteps 个有效点
    nValid = 0;
    validPts = zeros(totalSteps, 2);
    iAttempt = 0;
    while nValid < totalSteps
        iAttempt = iAttempt + 1;
        xCandidate = muX + sigmaX*randn;
        yCandidate = muY + sigmaY*randn;
        if xCandidate>=0 && xCandidate<=stairLength && ...
           yCandidate>=0 && yCandidate<=stairWidth
            nValid = nValid + 1;
            validPts(nValid, :) = [xCandidate, yCandidate];
        end
        % 若想避免无限循环，可加一个上限，但一般不会持续很久
    end
    % 用 validPts 替代最初的 pts
    pts = validPts;

    %% =========== 3. 预先计算每个踩踏点的椭圆边界(用于第三个图) =======
    % 由于第三张图要把所有踩踏点的椭圆都画上去，需要获得每个点的椭圆轮廓。
    % 椭圆参数方程：xEllipse = x0 + Rx*cos(t), yEllipse = y0 + Ry*sin(t)
    tVals = linspace(0, 2*pi, ellipsePlotResolution);
    allEllipseX = cell(totalSteps, 1);
    allEllipseY = cell(totalSteps, 1);
    for i = 1 : totalSteps
        x0 = pts(i,1);
        y0 = pts(i,2);
        xBoundary = x0 + ellipseRx*cos(tVals);
        yBoundary = y0 + ellipseRy*sin(tVals);
        allEllipseX{i} = xBoundary;
        allEllipseY{i} = yBoundary;
    end

    %% =========== 4. 创建三个 figure，用于在迭代中实时刷新 ==========

    % Figure 1: 平面热力图(2D)
    fig1 = figure('Name','平面热力图','Color','w');
    ax1  = axes('Parent',fig1);
    imagesc(ax1, xVals, yVals, h); 
    set(ax1,'YDir','normal');  % 修正 y 方向
    axis(ax1, 'equal'); 
    colorbar; 
    title(ax1, '楼梯高度平面分布（初始）');
    xlabel(ax1, 'x (m)'); ylabel(ax1, 'y (m)');

    % Figure 2: 三维立体图
    fig2 = figure('Name','楼梯3D磨损图','Color','w');
    ax2  = axes('Parent',fig2);
    surf(ax2, xVals, yVals, h, 'EdgeColor','none');
    axis(ax2, 'tight');
    shading(ax2, 'interp');
    colorbar(ax2);
    title(ax2, '楼梯表面高度（初始）');
    xlabel(ax2, 'x (m)'); ylabel(ax2, 'y (m)'); zlabel(ax2, 'Height (m)');
    view(ax2, 45,45);

    % Figure 3: 椭圆掩码分布图
    fig3 = figure('Name','椭圆掩码累计图','Color','w');
    ax3  = axes('Parent',fig3);
    hold(ax3, 'on');
    axis(ax3, [0 stairLength 0 stairWidth]);
    axis(ax3, 'equal');
    title(ax3, '已踩踏的所有椭圆掩码边界');
    xlabel(ax3, 'x (m)'); ylabel(ax3, 'y (m)');

    %% =========== 5. 迭代：每“脚”一次，扣减楼梯高度，定期刷新图 ======
    % 此处用一个 for 循环模拟从第1千年 ~ 第 totalMillennia 千年，总共 totalSteps 步

    for iStep = 1 : totalSteps

        % 读取本次踩踏点
        x0 = pts(iStep,1);
        y0 = pts(iStep,2);

        % 构建椭圆掩码，对楼梯表面进行磨损
        distEllipse = ((XX - x0)/ellipseRx).^2 + ((YY - y0)/ellipseRy).^2;
        inMask = (distEllipse <= 1);
        h(inMask) = h(inMask) - wearDepthPerStep;  % 扣减

        % 判断是否需要刷新绘图
        if (mod(iStep, plotInterval) == 0) || (iStep == 1) || (iStep == totalSteps)

            % ============ 更新 Figure 1：平面热力图 ============
            figure(fig1);
            cla(ax1); % 清除旧图
            imagesc(ax1, xVals, yVals, h);
            set(ax1,'YDir','normal');
            axis(ax1, 'equal');
            colorbar; 
            title(ax1, sprintf('楼梯高度平面分布：第 %d 步 (%.1f 千年)', ...
                   iStep, iStep/footStepsPerMillenni));
            xlabel(ax1, 'x (m)'); ylabel(ax1, 'y (m)');

            % ============ 更新 Figure 2：3D 表面图 ============
            figure(fig2);
            cla(ax2);
            surf(ax2, xVals, yVals, h, 'EdgeColor','none');
            shading(ax2, 'interp');
            colorbar(ax2);
            title(ax2, sprintf('楼梯3D高度：第 %d 步 (%.1f 千年)', ...
                   iStep, iStep/footStepsPerMillenni));
            xlabel(ax2, 'x (m)'); ylabel(ax2, 'y (m)'); zlabel(ax2, 'Height (m)');
            view(ax2, 45,45);

            % ============ 更新 Figure 3：叠加所有椭圆边界 ============
            figure(fig3);
            cla(ax3);
            hold(ax3, 'on');
            % 将前 iStep 次的椭圆都画出来
            for j = 1 : iStep
                plot(ax3, allEllipseX{j}, allEllipseY{j}, 'r-', 'LineWidth', 0.5);
            end
            axis(ax3, [0 stairLength 0 stairWidth]);
            axis(ax3, 'equal');
            title(ax3, sprintf('椭圆掩码累计：前 %d 步 (%.1f 千年)', ...
                   iStep, iStep/footStepsPerMillenni));
            xlabel(ax3, 'x (m)'); ylabel(ax3, 'y (m)');

            drawnow;  % 刷新图像
        end

    end

    %% =========== 6. 输出最终高度统计并结束 =========================
    disp('--------------------------------------------------');
    disp(['总迭代步数(总脚步): ', num2str(totalSteps)]);
    disp(['最终楼梯最小高度: ', num2str(min(h(:))), ' m']);
    disp(['最终楼梯平均高度: ', num2str(mean(h(:))), ' m']);
    disp('模拟结束.');
end

---

代码功能与说明

1. 参数定义

   • stairLength, stairWidth, initialHeight, Nx, Ny 等用于控制楼梯平面范围及网格精度。
   • totalMillennia（总千年数） 与 footStepsPerMillenni（每千年脚步数）相乘得到 totalSteps（总脚步数）。在大规模模拟时可酌情增大。
   • plotInterval 控制每隔多少步刷新/绘图一次，以免频繁绘图导致程序速度过慢。

2. 生成踩踏点 (二维正态分布)

   • 通过 mvnrnd 或者手动 randn 采样出候选点 (x, y)，并判断是否在 [0, stairLength] x [0, stairWidth] 区域内。
   • 这里写了一个循环保底，保证一定能采样到所需数量的点数。

3. 椭圆边界坐标

   • 为在第三张图里绘制所有踩踏掩码，需要提前计算好每个踩踏点对应的椭圆轮廓坐标（allEllipseX{i}, allEllipseY{i}）。
   • 椭圆边界方程： $$x_{\text{ellipse}} = x_0 + R_x\cos(t), \quad y_{\text{ellipse}} = y_0 + R_y\sin(t).$$

4. 主迭代过程

   • 循环 iStep = 1 : totalSteps，每一步取出该脚步的中心 (x0,y0)，构造椭圆掩码 inMask，并对楼梯高度矩阵 h 做磨损扣减。
   • 每当 mod(iStep, plotInterval)==0（或刚开始、或结束）时，刷新所有图：
     1. Figure 1：平面热力图（imagesc）。
     2. Figure 2：3D 表面图（surf）。
     3. Figure 3：将前 iStep 所有椭圆掩码的边界都绘制到同一个坐标系里。

5. 运行速度与可视化

   • 若脚步数 totalSteps 很大，并且 plotInterval 很小，则会大量绘图，程序可能会非常慢。可根据需要增大 plotInterval、或在模拟结束后再统一绘制结果。
   • 第三张图要求把每一脚的椭圆都叠加在同一图上，若脚步数过多，此图会非常密集，也会增加绘制时间。可视化时应酌情取较小的脚步数或更大的刷新间隔。

6. 时间单位/千年理解

   • 脚本中将脚步数拆分为 [每千年脚步数] x [千年数]。可以自行理解：如 footStepsPerMillenni=500 表示“每千年”发生 500 次踩踏，这仅是一个演示用的数字。
   • 在实际物理场景中，可以根据人流量估计“每天多少人”“一年多少人”，再乘以 1000 即可得到“千年”总踩踏次数，从而替换本代码的参数。

---

运行后即可在 MATLAB 中弹出三个窗口分别显示：

1. 平面热力图：楼梯不同网格点的高度颜色分布；
2. 三维立体图：可用鼠标旋转查看楼梯磨损表面的 3D 形态；
3. 掩码生成图：所有已经发生的踩踏椭圆边界叠加在一起的示意。

如需进一步扩展为 “bulk / edge 两种核函数” 或其他更复杂的摩擦磨损模型，可在本示例的基础上自行修改磨损公式与可视化方式。祝模拟顺利!

下面给出一个单窗口多子图的完整示例，演示如何在 MATLAB 中进行“以千年为单位”迭代的楼梯磨损模拟，并将2D 热力图、3D 表面图、以及所有踩踏椭圆掩码分布三个视图统一绘制在同一个 figure 的不同子图 (subplots) 上，随着迭代更新：

---

matlab
function simulateStairWear_SingleFigure()
% SIMULATESTAIRWEAR_SINGLEFIGURE
%
% 在同一个窗口内的 3 个子图实时更新：
%   (1) 楼梯表面热力图 (2D)
%   (2) 楼梯表面三维图 (3D)
%   (3) 已踩踏的所有椭圆掩码叠加图
%
% 其中“以千年为单位”的迭代，通过 [footStepsPerMillenni x totalMillennia]
% 控制总踩踏次数，每隔 plotInterval 步刷新一下图像。
%
% 如果一次踩踏总数很大，建议调大 plotInterval 以减少绘图开销。

    %% ================== 0. 可调参数 ======================
    % --- 楼梯与网格 ---
    stairLength       = 1.0;    % 楼梯 x 方向长度 (米)
    stairWidth        = 0.4;    % 楼梯 y 方向宽度 (米)
    initialHeight     = 0.3;    % 初始楼梯高度 (米)
    Nx                = 101;    % x方向网格数
    Ny                = 41;     % y方向网格数

    % --- 时间与脚步 (以千年为单位) ---
    totalMillennia       = 2;        % 总计多少个"千年"
    footStepsPerMillenni = 200;      % 每千年发生多少"脚步"
    totalSteps           = totalMillennia * footStepsPerMillenni;
    
    % --- 刷新绘图步长 ---
    plotInterval = 50;   % 每隔多少步刷新一次图像

    % --- 二维正态分布参数 (踩踏中心分布) ---
    muX    = 0.5;    % x方向均值
    muY    = 0.15;   % y方向均值
    sigmaX = 0.3;    % x方向标准差
    sigmaY = 0.2;    % y方向标准差

    % --- 椭圆掩码参数 ---
    ellipseRx = 0.02;   % 椭圆 x方向半径 (米)
    ellipseRy = 0.01;   % 椭圆 y方向半径 (米)

    % --- 每次踩踏产生的磨损深度 ---
    wearDepthPerStep = 5e-5;  % (米)

    % --- 绘制椭圆边界的精细度(用于子图3) ---
    ellipsePlotResolution = 36; 

    %% ============ 1. 建立楼梯表面网格并初始化高度 h =============

    xVals = linspace(0, stairLength, Nx);  % x方向坐标
    yVals = linspace(0, stairWidth,  Ny);  % y方向坐标
    [XX, YY] = meshgrid(xVals, yVals);     % (Ny x Nx)
    h = initialHeight * ones(Ny, Nx);      % 初始楼梯高度矩阵

    %% ============ 2. 生成随机踩踏点 (二维正态) =============
    %   需要 totalSteps 个点，每个点位于楼梯区域内

    mu    = [muX, muY];
    Sigma = [sigmaX^2, 0; 0, sigmaY^2];

    % 为确保能得到足够多在 [0,stairLength]x[0,stairWidth] 范围内的点，写一个循环
    pts = zeros(totalSteps,2);
    countValid = 0;
    while countValid < totalSteps
        candidate = mvnrnd(mu, Sigma);
        xC = candidate(1);
        yC = candidate(2);
        if (xC>=0 && xC<=stairLength && yC>=0 && yC<=stairWidth)
            countValid = countValid + 1;
            pts(countValid,:) = [xC,yC];
        end
    end

    %% ============ 3. 事先计算所有踩踏椭圆边界(用于子图3) =========
    %  这样每次刷新时，可将当前之前所有椭圆一并画上去

    tVals = linspace(0, 2*pi, ellipsePlotResolution);
    ellipseBordersX = cell(totalSteps,1);
    ellipseBordersY = cell(totalSteps,1);
    for i = 1 : totalSteps
        x0 = pts(i,1);
        y0 = pts(i,2);
        xBoundary = x0 + ellipseRx * cos(tVals);
        yBoundary = y0 + ellipseRy * sin(tVals);
        ellipseBordersX{i} = xBoundary;
        ellipseBordersY{i} = yBoundary;
    end

    %% ============ 4. 准备一个 figure，3 个 subplot =============
    mainFig = figure('Name','楼梯磨损模拟','Color','w','Position',[100,100,1400,450]);

    % subplot(1,3,1): 平面热力图
    ax1 = subplot(1,3,1,'Parent',mainFig);
    imagesc(ax1, xVals, yVals, h);
    set(ax1,'YDir','normal');    % 让 y 轴随增大向上
    axis(ax1,'equal'); axis(ax1,'tight');
    colorbar(ax1);
    xlabel(ax1,'x (m)'); ylabel(ax1,'y (m)');
    title(ax1,'初始楼梯高度 - 2D');
    
    % subplot(1,3,2): 三维立体图
    ax2 = subplot(1,3,2,'Parent',mainFig);
    surf(ax2, xVals, yVals, h, 'EdgeColor','none');
    shading(ax2,'interp');
    colorbar(ax2);
    xlabel(ax2,'x (m)'); ylabel(ax2,'y (m)'); zlabel(ax2,'Height');
    title(ax2,'初始楼梯 - 3D');
    view(ax2,45,45);

    % subplot(1,3,3): 踩踏椭圆分布图
    ax3 = subplot(1,3,3,'Parent',mainFig);
    hold(ax3,'on'); 
    axis(ax3,[0 stairLength 0 stairWidth]); 
    axis(ax3,'equal');
    xlabel(ax3,'x (m)'); ylabel(ax3,'y (m)');
    title(ax3,'椭圆掩码分布 (初始)');

    drawnow;

    %% ============ 5. 进行迭代，每步扣减楼梯表面 =============
    for iStep = 1 : totalSteps

        % 当前脚的中心
        x0 = pts(iStep,1);
        y0 = pts(iStep,2);

        % 构造椭圆掩码
        distEllipse = ((XX - x0)/ellipseRx).^2 + ((YY - y0)/ellipseRy).^2;
        inMask = (distEllipse <= 1);

        % 扣减高度
        h(inMask) = h(inMask) - wearDepthPerStep;

        % 每隔 plotInterval 步 或最后一步 刷新一次可视化
        if mod(iStep, plotInterval)==0 || iStep==totalSteps || iStep==1
            
            % ---------- 更新子图1：平面热力图 ----------
            subplot(1,3,1);
            cla(ax1);
            imagesc(ax1, xVals, yVals, h);
            set(ax1,'YDir','normal');
            axis(ax1,'equal'); axis(ax1,'tight');
            colorbar(ax1);
            xlabel(ax1,'x (m)'); ylabel(ax1,'y (m)');
            title(ax1, sprintf('楼梯高度 - 2D (第 %d 步, %.2f 千年)', ...
                  iStep, iStep/footStepsPerMillenni));

            % ---------- 更新子图2：3D 表面图 ----------
            subplot(1,3,2);
            cla(ax2);
            surf(ax2, xVals, yVals, h, 'EdgeColor','none');
            shading(ax2,'interp');
            colorbar(ax2);
            view(ax2,45,45);
            xlabel(ax2,'x (m)'); ylabel(ax2,'y (m)'); zlabel(ax2,'Height');
            title(ax2, sprintf('楼梯 - 3D (第 %d 步, %.2f 千年)', ...
                  iStep, iStep/footStepsPerMillenni));

            % ---------- 更新子图3：椭圆掩码累计图 ----------
            subplot(1,3,3);
            cla(ax3);
            hold(ax3,'on');
            for j = 1 : iStep
                plot(ax3, ellipseBordersX{j}, ellipseBordersY{j}, 'r-', 'LineWidth', 0.5);
            end
            axis(ax3,[0 stairLength 0 stairWidth]); 
            axis(ax3,'equal');
            xlabel(ax3,'x (m)'); ylabel(ax3,'y (m)');
            title(ax3, sprintf('已踩踏椭圆 (前 %d 步)', iStep));

            drawnow;  % 刷新到屏幕
        end

    end

    %% ============ 6. 输出一些统计信息 =============
    minH = min(h(:));
    avgH = mean(h(:));
    disp('---------------------------------------------');
    disp(['模拟结束，共有 ', num2str(totalSteps), ' 次踩踏.']);
    disp(['最终楼梯最小高度: ', num2str(minH), ' m']);
    disp(['最终楼梯平均高度: ', num2str(avgH), ' m']);
end

---

代码说明

1. 单窗口 + 三个子图
   使用一个 figure（mainFig），在其中 subplot(1,3,1)、subplot(1,3,2)、subplot(1,3,3) 分别放置 2D 热力图、3D 曲面图 和 椭圆掩码累计图。在每次刷新时，依次 cla（清除旧图）后重新绘制。

2. 以千年为单位

   • totalMillennia 表示要模拟的“千年”数，例如 2 即 2 个千年。
   • footStepsPerMillenni 表示每千年发生的“脚步”数。
   • 二者相乘得到 totalSteps，即本次模拟的总脚步数。迭代结束后，也可根据真实需求进行缩放。

3. plotInterval

   • 为避免每一步都画图（速度慢），仅在 plotInterval 的倍数步数以及最后一步时进行可视化更新。
   • 如果想更加频繁更新，可将 plotInterval 调小；若想加速模拟（少刷新），可将其调大。

4. 第三个子图（踩踏掩码累计）

   • 由于用户要求“第三个图要显示所有已踩踏的椭圆”，因此在刷新时，会将当前脚之前的所有椭圆边界一并绘制。
   • 如果 totalSteps 很大，此步会比较耗时，可根据需要缩小可视化频率或只显示最近的一部分踩踏点。

5. 其他

   • 如果想让脚步更加随机地分布在楼梯平面，可以改用更高的 “脚步采样总数 + 过滤” 或用截断正态分布等方式。
   • 可在迭代结束后再一次性绘制最终图，以避免重复绘图的性能消耗。

这样就实现了一个单窗口、分 3 个子图的动态更新楼梯磨损模拟，涵盖2D热力图、3D表面图和所有踩踏椭圆掩码叠加图。

下面给出一个思路示例，展示如何将势函数 $\phi$ 纳入到人流踩踏分布中，从而让人们“更倾向踩到好踩的地方”。核心思想是：

1. 给定楼梯表面高度 $h(x,y,t)$ 后，根据某种规则计算势函数 $\phi(x,y,t)$。
2. 用指数形式 $\exp[-\,\gamma\,\phi]$ 来得到“未归一化”的加权分布 $\widetilde{P}$。
3. 归一化到期望的总流量 $\Lambda_{\text{bulk}}(t)$，或在离散网格上直接归一化为一个概率分布，再用随机抽样来确定每一步的踩踏点。
4. 在数值上，每一次迭代都可以（可选）根据最新的楼梯表面高度 $h$（及其斜度等）来更新 $\phi$ 和分布，从而模拟“楼梯磨损 -> 人们踩踏方式改变 -> 进一步磨损”的循环。

下面给出一个完整的 MATLAB 演示脚本，与之前的 simulateStairWear 类似，只是把踩踏点的产生从“简单正态分布”改为“基于势函数的加权分布”。并在单个窗口里绘制三幅子图实时更新。

注意

• 这只是示例，供你理解如何将 $\phi$ 融合到踩踏分布中。实际中，$\phi$ 的定义会更复杂，比如依赖楼梯表面的斜度、凹凸度、湿滑程度等。
• 代码示例里演示了“每一步都根据最新 $h$ 重新计算势函数并做抽样”，如果计算量大，可以视需求改为“隔几步才更新势函数”、“在外层迭代做更高层次更新”等。
• 若要在连续坐标 $(x,y)$ 中采样，需要将离散网格概率转换为一系列小单元的累积分布，这里演示的是“在网格中心附近”进行踩踏抽样的离散化方法。

---

matlab
function simulateStairWear_withPotential()
% SIMULATESTAIRWEAR_WITHPOTENTIAL
%
% 演示如何将势函数 phi 纳入到人群踩踏分布中，使人更倾向踩到“势函数较低”的地方。
% 
% 1) 每次根据当前楼梯表面高度 h(x,y)，计算 phi(i,j)；
% 2) 用 exp(-gamma * phi) 得到加权分布后离散采样，
% 3) 在被采样的网格附近产生一次脚步，然后用椭圆掩码扣减楼梯高度，
% 4) 循环迭代并在单窗口的三个子图里实时更新。
%
% 说明：phi 的定义可自由扩展，这里仅以“楼梯局部斜度”为例进行简化示意。

    %% =========== 0. 可调参数 ===========

    % 楼梯尺寸和网格
    stairLength   = 1.0;    % 楼梯 x 方向长度 (米)
    stairWidth    = 0.4;    % 楼梯 y 方向宽度 (米)
    initialHeight = 0.3;    % 初始楼梯高度 (米)
    Nx            = 101;    % x方向网格数
    Ny            = 41;     % y方向网格数

    % 时间/迭代步骤
    totalMillennia       = 2;         % 总共模拟多少个“千年”
    footStepsPerMillenni = 200;       % 每千年踩踏步数
    totalSteps           = totalMillennia * footStepsPerMillenni;
    plotInterval         = 50;        % 每隔多少步刷新一次图像

    % 每次踩踏的椭圆掩码半径
    ellipseRx = 0.02;
    ellipseRy = 0.01;

    % 每次踩踏的磨损深度
    wearDepthPerStep = 5e-5;

    % 势函数 phi 相关
    gamma = 10;  % 人对势函数的敏感度参数（越大代表越不愿踩“势”高的地方）

    % 椭圆可视化分辨率（第三张图）
    ellipsePlotResolution = 36;

    %% =========== 1. 建立楼梯表面网格及高度矩阵 ===========

    xVals = linspace(0, stairLength, Nx);
    yVals = linspace(0, stairWidth,  Ny);
    [XX, YY] = meshgrid(xVals, yVals);   % (Ny x Nx)
    h = initialHeight * ones(Ny, Nx);    % 初始楼梯高度

    %% =========== 2. 为绘制椭圆边界做准备（先不定踩踏点） ===========
    % 因为本例中，踩踏点要“每一步根据势函数分布随机生成”，
    % 所以只有在生成点之后，才能得知椭圆位置。
    % 这里先预留一个 cell，用来保存每一步的椭圆坐标，以便画在第三幅图中。

    ellipseBordersX = cell(totalSteps,1);
    ellipseBordersY = cell(totalSteps,1);

    %% =========== 3. 创建 figure 及 3 个子图 ===========
    mainFig = figure('Name','楼梯磨损(势函数分布)','Color','w',...
                     'Position',[100,100,1400,450]);

    ax1 = subplot(1,3,1,'Parent',mainFig); % 2D热力图
    ax2 = subplot(1,3,2,'Parent',mainFig); % 3D表面图
    ax3 = subplot(1,3,3,'Parent',mainFig); % 椭圆掩码分布

    %% =========== 4. 迭代过程：每一步都重新计算phi，并离散采样 ===========

    for iStep = 1 : totalSteps

        % ---- 4.1 根据当前楼梯高度 h 计算势函数 phi ----
        %    这里给出一个简化示例：phi(i,j) = 局部斜度(梯度幅值)。
        %    在实际中，你可以用更复杂的公式：phi = f( h, grad(h), 湿滑度, 等 ).
        %
        %    简单做法：用有限差分近似梯度:
        %       dh/dx(i,j) ~ [ h(i,j+1) - h(i,j-1) ] / (2*dx)
        %       dh/dy(i,j) ~ [ h(i+1,j) - h(i-1,j) ] / (2*dy)
        %    注意边界可采用向前/向后差分，或假设边界斜度为0，仅作演示。
        
        % 网格步长
        dx = xVals(2) - xVals(1);
        dy = yVals(2) - yVals(1);
        
        % 初始化势函数矩阵
        phi = zeros(Ny, Nx);

        for iy = 2 : (Ny-1)
            for ix = 2 : (Nx-1)
                dhdx = (h(iy, ix+1) - h(iy, ix-1)) / (2*dx);
                dhdy = (h(iy+1, ix) - h(iy-1, ix)) / (2*dy);
                slopeMag = sqrt(dhdx^2 + dhdy^2);
                % 在此处可以组合更多因子，如: phi(iy,ix) = slopeMag + alpha*(...)
                phi(iy, ix) = slopeMag;
            end
        end
        % 边界简化：令 phi=0 或复制相邻值
        phi(1,:)   = phi(2,:);
        phi(end,:) = phi(end-1,:);
        phi(:,1)   = phi(:,2);
        phi(:,end) = phi(:,end-1);

        % ---- 4.2 构造加权分布 w_ij = exp( - gamma * phi_ij ) ----
        wMat = exp(-gamma * phi);

        % ---- 4.3 归一化：把 wMat 当作“离散概率” ----
        %   pMat(i,j) = wMat(i,j) / sum(wMat(:))
        sumW = sum(wMat(:));
        pMat = wMat / sumW;

        % ---- 4.4 从离散分布 pMat 中采样一个网格点 (ix, iy) 作为踩踏中心 ----
        %   方法：将 pMat(:) 展平成一维，再用 randsample 按权重抽样
        flatP = pMat(:);
        idxChosen = randsample(1 : (Nx*Ny), 1, true, flatP);
        % 将 idxChosen 转回 (iy, ix)
        [iyChosen, ixChosen] = ind2sub([Ny, Nx], idxChosen);

        %   把这个网格中心当成 (x0, y0)，若要更精细，也可加一点 [0,dx],[0,dy] 的随机
        x0 = xVals(ixChosen);
        y0 = yVals(iyChosen);

        % ---- 4.5 在 (x0,y0) 处构建椭圆掩码并扣减楼梯高度 ----
        distEllipse = ((XX - x0)/ellipseRx).^2 + ((YY - y0)/ellipseRy).^2;
        inMask = (distEllipse <= 1);
        h(inMask) = h(inMask) - wearDepthPerStep;

        % ---- 4.6 保存椭圆边界坐标（用于第三张子图可视化） ----
        tVals = linspace(0,2*pi, ellipsePlotResolution);
        ellX = x0 + ellipseRx*cos(tVals);
        ellY = y0 + ellipseRy*sin(tVals);
        ellipseBordersX{iStep} = ellX;
        ellipseBordersY{iStep} = ellY;

        % ---- 4.7 每隔若干步刷新图像 ----
        if (mod(iStep, plotInterval)==0) || (iStep==1) || (iStep==totalSteps)

            % ==== 子图1：平面热力图 ====
            subplot(ax1);  cla(ax1);
            imagesc(ax1, xVals, yVals, h);
            set(ax1,'YDir','normal');
            axis(ax1,'equal'); axis(ax1,'tight');
            colorbar(ax1);
            xlabel(ax1,'x (m)'); ylabel(ax1,'y (m)');
            title(ax1, sprintf('楼梯高度 - 2D (第 %d 步)', iStep));

            % ==== 子图2：3D 表面图 ====
            subplot(ax2);  cla(ax2);
            surf(ax2, xVals, yVals, h, 'EdgeColor','none');
            shading(ax2,'interp');
            colorbar(ax2);
            view(ax2,45,45);
            xlabel(ax2,'x (m)'); ylabel(ax2,'y (m)'); zlabel(ax2,'Height (m)');
            title(ax2, sprintf('楼梯3D (第 %d 步)', iStep));

            % ==== 子图3：迄今为止所有椭圆掩码边界 ====
            subplot(ax3);  cla(ax3);
            hold(ax3,'on');
            for k = 1 : iStep
                plot(ax3, ellipseBordersX{k}, ellipseBordersY{k}, 'r-', 'LineWidth',0.5);
            end
            axis(ax3,[0 stairLength 0 stairWidth]);
            axis(ax3,'equal');
            xlabel(ax3,'x (m)'); ylabel(ax3,'y (m)');
            title(ax3, sprintf('已踩踏椭圆 (前 %d 步)', iStep));
            drawnow;

        end

    end

    %% =========== 5. 输出一些统计信息 =============
    disp('----------------------------------');
    disp(['共完成 ', num2str(totalSteps), ' 次踩踏.']);
    disp(['最终楼梯最小高度: ', num2str(min(h(:))), ' m']);
    disp(['最终楼梯平均高度: ', num2str(mean(h(:))), ' m']);
    disp('模拟结束！');

end

---

主要改动要点

1. 势函数 $\phi$ 的计算

   • 在示例中，用了一个简单的“楼梯表面斜度”来定义 $\phi$，即局部斜度越大，势函数越高，人就越不愿踩。
   • 具体做法：用有限差分计算 $\frac{\partial h}{\partial x}$ 和 $\frac{\partial h}{\partial y}$，令 $\phi=\sqrt{(\partial h/\partial x)^2 + (\partial h/\partial y)^2}$。
   • 边界处做了简单处理（复制邻近值）。
   • 在实际应用中，你可以写一个更复杂的函数 phi = potentialFunction(h, ...)，结合台阶磨损、湿滑、材料等指标。

2. 加权分布 $\exp[-\gamma\,\phi]$

   • 先用 wMat = exp(- gamma* phi) 得到未归一化的权重矩阵。
   • 然后 pMat = wMat / sum(wMat(:)) 得到离散概率矩阵。

3. 离散采样

   • 通过 randsample(1:Nx*Ny, 1, true, flatP) 在整张楼梯网格上按概率抽取一个点 (ixChosen, iyChosen)。
   • 将其坐标转换回实际楼梯面上的 (x0, y0)。
   • 也可以在离散点基础上再加一点小随机量，使踩踏位置在当前网格内更加连续。

4. 递回迭代

   • 由于楼梯高度 h 在每一步踩踏后都会变化，势函数 $\phi$ 也会随之改变，因此示例中在每一步都重新计算。
   • 如果计算量大，可以改为“每隔 N 步更新一次势函数”或类似策略。

5. 单窗口 3 子图

   • 与之前示例的多窗口区别在于，这里用 subplot(1,3,X) 放置2D热力图、3D表面图和所有椭圆掩码累计图在同一窗口内，每隔 plotInterval 步更新一次。

---

通过这种方式，就将“楼梯表面斜度(或其他物理因素)”映射为“势函数 $\phi$”，再用 $\exp[-\gamma\,\phi]$ 来偏移人流分布，让人们更倾向于踩到平缓或“好踩”的区域，从而可以模拟在磨损过程中的“人流自适应选择”。

下面给出一个思路示例，展示如何在踩踏椭圆掩码下，检测“局部较高的尖峰”并给予额外磨损扣减，以模拟“脚步踩到局部突起（尖峰）时会削弱更多”的效果。

基本思路

1. 识别踩踏区域：和之前一样，通过 $(x_0,y_0)$ 及椭圆掩码判断哪些网格点 $(i,j)$ 处于踩踏区域。
2. 检测尖峰：取该踩踏区域内的所有高度 h(inMask) 进行统计。例如：
   • 最大高度与区域平均高度之差（$\text{peakDiff} = \max(h) - \text{mean}(h)$)。
   • 或者直接检测 $\max(h)$ 是否大于某阈值。
3. 根据尖峰大小，调整磨损深度：如果发现“尖峰”超过了给定阈值，就额外加大（或乘以某因子）本次磨损深度，达到“多削弱一些尖峰”的效果。

下面是一个示例化的简单实现，基于之前的 simulateStairWear_withPotential 框架，只需在“椭圆掩码”部分增加尖峰检测与额外磨损逻辑即可。

---

matlab
function simulateStairWear_spikeReduction()
% SIMULATESTAIRWEAR_SPIKEREDUCTION
%
% 演示如何在踩踏椭圆掩码中检测局部的高度尖峰(peak)，若明显凸起，令磨损更强。
% 基础逻辑：
%   1) 每步我们依然根据势函数来采样踩踏点 (x0, y0)；
%   2) 找出椭圆掩码 inMask 下的所有点；
%   3) 计算尖峰指标 (如 max(h区域) - mean(h区域))，若超过阈值则乘以额外因子；
%   4) 对该区域扣减高度(磨损)。
%
% 与之前类似，也在单窗口3个子图中实时更新。

    %% =========== 0. 可调参数 ===========
    % 楼梯尺寸与网格
    stairLength   = 1.0;
    stairWidth    = 0.4;
    initialHeight = 0.3;
    Nx            = 101;
    Ny            = 41;

    % 时间/步骤
    totalMillennia       = 2;
    footStepsPerMillenni = 200;
    totalSteps           = totalMillennia * footStepsPerMillenni;
    plotInterval         = 50;

    % 每次踩踏的椭圆掩码
    ellipseRx = 0.02;
    ellipseRy = 0.01;

    % 基础磨损深度
    wearDepthPerStep = 5e-5;

    % 尖峰检测参数 (示例)
    peakThreshold = 0.002;   % 若 (maxHeight - avgHeight) 超过2mm则视为尖峰
    spikeFactor   = 2.0;     % 出现尖峰时，磨损加倍

    % 势函数 phi 相关(示例:对楼梯斜度敏感)
    gamma = 10;  

    % 椭圆可视化分辨率(第三张图)
    ellipsePlotResolution = 36;

    %% =========== 1. 网格+高度初始化 ===========
    xVals = linspace(0, stairLength, Nx);
    yVals = linspace(0, stairWidth,  Ny);
    [XX, YY] = meshgrid(xVals, yVals);
    h = initialHeight * ones(Ny, Nx);

    %% =========== 2. 预留椭圆边界存储(用于可视化) ===========
    ellipseBordersX = cell(totalSteps,1);
    ellipseBordersY = cell(totalSteps,1);

    %% =========== 3. 创建 figure 并分三个子图 ===========
    mainFig = figure('Name','楼梯磨损(检测尖峰)','Color','w',...
        'Position',[100,100,1400,450]);
    ax1 = subplot(1,3,1,'Parent',mainFig); % 2D图
    ax2 = subplot(1,3,2,'Parent',mainFig); % 3D图
    ax3 = subplot(1,3,3,'Parent',mainFig); % 椭圆掩码分布

    %% =========== 4. 迭代踩踏过程 ===========
    for iStep = 1 : totalSteps

        % ---- 4.1 计算当前楼梯势函数 phi (基于斜度的简化例) ----
        dx = xVals(2)-xVals(1);
        dy = yVals(2)-yVals(1);

        phi = zeros(Ny,Nx);
        for iy=2:Ny-1
            for ix=2:Nx-1
                dhdx = (h(iy, ix+1)-h(iy, ix-1))/(2*dx);
                dhdy = (h(iy+1, ix)-h(iy-1, ix))/(2*dy);
                phi(iy, ix) = sqrt(dhdx^2 + dhdy^2);
            end
        end
        % 边界处理(简化)
        phi(1,:) = phi(2,:);
        phi(end,:)=phi(end-1,:);
        phi(:,1)=phi(:,2);
        phi(:,end)=phi(:,end-1);

        % ---- 4.2 构造加权分布 exp(- gamma*phi)，离散化采样 ----
        wMat = exp(-gamma * phi);
        sumW = sum(wMat(:));
        pMat = wMat / sumW;
        flatP = pMat(:);
        idxChosen = randsample(1:(Nx*Ny),1,true,flatP);
        [iyChosen, ixChosen] = ind2sub([Ny,Nx], idxChosen);

        x0 = xVals(ixChosen);
        y0 = yVals(iyChosen);

        % ---- 4.3 椭圆掩码 + 检测区域尖峰 ----
        distEllipse = ((XX - x0)/ellipseRx).^2 + ((YY - y0)/ellipseRy).^2;
        inMask = (distEllipse <= 1);
        
        localHeights = h(inMask);
        hMax = max(localHeights);
        hAvg = mean(localHeights);
        peakDiff = hMax - hAvg;  % 尖峰指标

        % 若尖峰差值超过阈值, 则加大磨损
        if peakDiff > peakThreshold
            wearFactor = spikeFactor;  % 例如加倍
        else
            wearFactor = 1.0;
        end

        % ---- 4.4 磨损扣减 ----
        h(inMask) = h(inMask) - (wearDepthPerStep * wearFactor);

        % ---- 4.5 记录椭圆边界(用于第三张子图) ----
        tVals = linspace(0,2*pi,ellipsePlotResolution);
        ellX = x0 + ellipseRx*cos(tVals);
        ellY = y0 + ellipseRy*sin(tVals);
        ellipseBordersX{iStep} = ellX;
        ellipseBordersY{iStep} = ellY;

        % ---- 4.6 每隔 plotInterval 步刷新图像 ----
        if (mod(iStep,plotInterval)==0) || (iStep==1) || (iStep==totalSteps)
            % ========== 子图1: 2D 热力图 ==========
            subplot(ax1); cla(ax1);
            imagesc(ax1, xVals, yVals, h);
            set(ax1,'YDir','normal');
            axis(ax1,'equal'); axis(ax1,'tight');
            colorbar(ax1);
            xlabel(ax1,'x (m)'); ylabel(ax1,'y (m)');
            title(ax1, sprintf('楼梯高度-2D (步 %d)', iStep));

            % ========== 子图2: 3D 表面图 ==========
            subplot(ax2); cla(ax2);
            surf(ax2, xVals, yVals, h, 'EdgeColor','none');
            shading(ax2,'interp');
            colorbar(ax2);
            view(ax2,45,45);
            xlabel(ax2,'x (m)'); ylabel(ax2,'y (m)'); zlabel(ax2,'Height');
            title(ax2, sprintf('楼梯3D (步 %d)', iStep));

            % ========== 子图3: 椭圆掩码累计图 ==========
            subplot(ax3); cla(ax3);
            hold(ax3,'on');
            for k=1:iStep
                plot(ax3, ellipseBordersX{k}, ellipseBordersY{k}, 'r-','LineWidth',0.5);
            end
            axis(ax3,[0 stairLength 0 stairWidth]);
            axis(ax3,'equal');
            xlabel(ax3,'x (m)'); ylabel(ax3,'y (m)');
            title(ax3, sprintf('已踩踏椭圆 (前 %d 步)', iStep));

            drawnow;
        end

    end

    %% =========== 5. 输出信息 ===========
    disp('-------------------------------------');
    disp(['完成 ', num2str(totalSteps), ' 次踩踏.']);
    disp(['最终最小高度: ', num2str(min(h(:)))]);
    disp(['最终平均高度: ', num2str(mean(h(:)))]);
    disp('模拟结束！');
end

---

关键新增逻辑

1. 尖峰检测

   matlab
   localHeights = h(inMask);  % 椭圆区域内的高度
   hMax = max(localHeights);
   hAvg = mean(localHeights);
   peakDiff = hMax - hAvg;

   • 这里用 $\max - \text{mean}$ 代表“尖峰程度”，当然也可定义其他指标，例如 $\max(h)$ 与某绝对阈值比较。
2. 阈值与加倍因子
   • 设定 peakThreshold（例如 0.002 m，即 2 mm），当 peakDiff > peakThreshold 时，就把磨损深度乘以 spikeFactor（例如 2.0）。
   • 可以根据材料、实际需求调整。也可以用一个更细腻的连续函数，而非简单地二分判断。
3. 额外扣减

   matlab
   if peakDiff > peakThreshold
       wearFactor = spikeFactor;
   else
       wearFactor = 1.0;
   end
   h(inMask) = h(inMask) - (wearDepthPerStep * wearFactor);

   • 当检测到尖峰，一次踩踏就会额外多扣减一些高度。这可以模拟“局部突出的台阶棱角容易被踩得更快磨损”。

---

可扩展性

1. 连续削弱因子
   • 目前的示例是如果超过阈值则直接加倍；也可以让“峰值差”越大，对应磨损因子越高。例如： $$\text{wearFactor} = 1 \;+\; k \times (\text{peakDiff} - \text{threshold})^+$$
   其中 $(x)^+ = \max(x,0)$ 表示大于0才生效，$k$ 为某个系数。
2. 本地极大检测
   • 如果想更准确识别“尖峰”，可以在椭圆区域内再对某些网格点做更精细的局部极大值检测，而不是单纯看 $\max - \text{mean}$。
   • 例如只对距离中心较近的地方检查，或者看梯度大小等。
3. 多重阈值
   • 也可以设置多个区间，如“尖峰 < 1 mm”不加成，“1~3 mm”加成1.5倍，“>3 mm”加成2倍等。
4. 与势函数配合
   • 在示例中，我们先按照势函数($\phi$)采样踩踏点，再进行“尖峰检查”来加大磨损。你也可以把“是否有凸起”这个信息直接融入势函数，让高凸起区域拥有更高的势函数，降低踩踏概率，从而形成双向耦合： 凸起 -> 人更不想踩 -> 踩踏更少 -> 磨损却可能更慢…… 等多种行为。

通过上述方法，即可在脚步踩踏的椭圆区域内，对“局部尖峰”做加强磨损处理，模拟实际使用中“突起位置更容易被迅速磨平”的现象。