CN102609560B - 一种3d任意粗糙表面的数字化模拟方法 - Google Patents

Abstract

本发明公开了一种3D任意粗糙表面的数字化模拟方法，提高了模拟效率与统计参数的模拟精度。通过白噪声序列的反傅里叶变换获得随机初始相角序列，利用初始相角序列获得新的白噪声及其傅里叶变换；通过指定的自相关函数，进行离散与傅里叶变换等处理，获得高斯粗糙表面高度序列的功率谱密度与传递函数；利用频域点乘并求反傅里叶变换的方法完成高斯表面高度序列的模拟；在此基础上，通过指定偏斜度与峰度等高度分布统计特征参数，利用Pearson与Johnson非高斯转换***相结合，生成非高斯粗糙表面；若偏斜度与峰度模拟精度不合格，则更新相角序列与白噪声的傅里叶变换，重新进行高斯滤波与非高斯转换，直到满足给定精度要求。

Description

一种3D任意粗糙表面的数字化模拟方法

技术领域

本发明涉及一种3D微观任意粗糙表面的数字模拟方法，更具体地说，本发明，根据连接表面点的高度分布前四阶矩统计特征与水平方向自相关函数统计规律，可以高效准确地构建具有相同统计特征的粗糙表面。

背景技术

高速、精密加工技术的发展及产品可靠度要求的不断提高，推动着粗糙表面的加工制造、检测、数字化模拟与性能预测等领域的快速发展。

近年来，国内外在粗糙表面的数字模拟领域的研究，主要集中在针对粗糙高度符合高斯分布的粗糙表面方面，在非高斯粗糙表面的数字化模拟方面，已能成功将高斯分布的数据转化为非高斯数据，出现了基于Johnson、Pearson等非高斯转换***。但在非高斯粗糙表面的数字化模拟当中，经非高斯转换后的非高斯序列还需要开展滤波，给模拟序列的偏斜度、峰度等统计特征参数引入了新的误差，同时也改变了偏斜度-峰度的对应的模拟区域。从而导致有些统计特征的粗糙表面不能通过Johnson、Pearson等转换***转换与滤波生成相应的非高斯序列，即使能转换，也会在滤波阶段引入新的误差，影响滤波的非高斯表面数字化模拟的精度。

考虑到滤波（滤波系数为h(k，l)）的影响，Johnson、Pearson等非高斯转换***中生成的非高斯序列粗糙表面的偏斜度S_kη、峰度K_uη与滤波后的粗糙表面的偏斜度S_kZ、峰度K_uZ存在如下关系：

$\left\{\begin{array}{c}{S}_{\mathrm{kz}}=\frac{{\mathrm{\Σ}}_{i=1}^{\mathrm{mn}}{\mathrm{\θ}}_i^3}{{\left({\mathrm{\Σ}}_{i=1}^{\mathrm{mn}}{\mathrm{\θ}}_i^2\right)}^{3/2}}{S}_{\mathrm{k\η}}\\ K_{\mathrm{uz}}=\frac{K_{\mathrm{u\η}}{\mathrm{\Σ}}_{i=1}^{\mathrm{mn}}{\mathrm{\θ}}_i^4+6{\mathrm{\Σ}}_{i=1}^{\mathrm{mn}-1}{\mathrm{\Σ}}_{j=i+1}^{\mathrm{mn}}{\mathrm{\θ}}_i^2{\mathrm{\θ}}_j^2}{{\left({\mathrm{\Σ}}_{i=1}^{\mathrm{mn}}{\mathrm{\θ}}_i^2\right)}^2}\end{array}\begin{array}{c}{\mathrm{\θ}}_i=h(k,l),k=\mathrm{1,2}...,m\\ l=1,...,n,i=(k-1)m+l\end{array}\right.$

即使目标偏斜度与峰度满足非高斯转换必要条件：但在考虑滤波系数与自相关函数的相互影响后，S_kη、K_uη有可能不再满足非高斯转换的必要条件：从而导致不能通过Johnson或Pearson等非高斯转换***生成非高斯序列；即使能转换，也会在滤波阶段引入新的误差，影响滤波的非高斯表面数字化模拟的精度。

现实中，不同的加工方法会产生具有不同统计特征的粗糙表面，粗糙表面的数字化模拟模型作为结合面性能分析的载体，其模拟的准确性直接影响相应结合面的接触性能预测。本研究利用一定的相角序列生成白噪声，通过白噪声的高斯滤波与高斯表面高度序列的非高斯转换***相结合，形成满足高度分布统计特征与自相关函数的非高斯高度序列，高斯或非高斯粗糙表面的数字化模拟方法为结合面性能研究奠定了基础。

发明内容

本发明的主要目的是提供一种高效高精度实现任意粗糙表面的数字化模拟方法，减少了粗糙表面数字化模拟中的滤波误差环节，在满足自相关函数、高度分布前四阶矩等统计特征的前提下，提高了模拟粗糙表面高度分布统计特征的精度与数字化模拟的效率。

本发明通过以下技术方案实现，主要步骤包括：

一种3D任意粗糙表面的数字化模拟方法，其特征在于，包含下述步骤：

1)当任意粗糙表面的偏斜度值约等于3且峰度值约等于0时，其表面形貌符合高斯分布，否则，其表面形貌符合非高斯分布；

2)如果步骤1）中粗糙表面形貌为高斯分布，利用三维高斯粗糙表面的模拟方法获得粗糙表面高度序列z，完成粗糙表面的数字化模拟；

3)如果步骤1）中粗糙表面形貌为非高斯分布，利用三维非高斯粗糙表面的模拟方法获得粗糙表面高度序列z，完成粗糙表面的数字化模拟。

所述的三维高斯粗糙表面的模拟方法，具体包含以下步骤：

1)通过给定的自相关函数生成离散的自相关函数矩阵R(m,n)，其中m、n分别表示x、y方向取点的数目；

2)对离散的自相关函数矩阵R(m,n)进行快速傅里叶变换，获得高斯粗糙表面的功率谱密度P(m,n)；

$P(1+k,1+l)=\left|\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}R(r+1,s+1)e^{-2\mathrm{\πi}(\mathrm{kr}/m+\mathrm{ls}/n)}\right|\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

3)假设白噪声的功率谱密度常数C等于1，利用获得的功率谱密度，计算传递函数H(m,n)；

$H=\sqrt{P}$

4)利用白噪声生成器，获得白噪声序列η(m,n)，利用反傅里叶变换生成相角序列Φ；

$\mathrm{\Φ}(k+1,l+1)={\tan }^{-1}\left(\frac{-\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}\mathrm{\η}(r+1,s+1)\sin (2\mathrm{\πkr}/m+2\mathrm{\πls}/n)}{-\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}\mathrm{\η}(r+1,s+1)\cos (2\mathrm{\πkr}/m+2\mathrm{\πls}/n)}\right)\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

5)以步骤4）获得的相角序列Φ，计算相角相关序列exp(iΦ)，并对其进行反傅里叶变换，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)；

$\mathrm{\η}(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}\exp (\mathrm{i\Φ}+\frac{2\mathrm{\πikr}}{m}+\frac{2\mathrm{\πils}}{n})\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

$A(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}\mathrm{\η}(r+1,s+1)e^{-2\mathrm{\πi}(\mathrm{kr}/m+\mathrm{ls}/n)}\≈\mathrm{mn}\·e^{\mathrm{i\Φ}}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

6)以步骤3）获得的传递函数H(m,n)，步骤5）获得的白噪声傅里叶变换序列A(m,n)，求序列H与序列A的点积运算，并对点积运算结果进行反傅里叶变换，获得高斯表面的高度序列z₀(m,n)；

$z_0(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}(A.*H)e^{2\mathrm{\πi}(\mathrm{kr}/m+\mathrm{ls}/n)}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

7)以步骤6）生成的高斯高度序列z₀，对其进行高度方向平移与缩放，获得满足均值与标准偏差的高斯高度序列z，完成高斯粗糙表面的数字化模拟。

所述的三维非高斯粗糙表面的模拟方法，具体包含以下步骤：

1)以获得的高斯高度序列z，利用高斯序列的标准化方法，在高度方向进行平移与缩放，获得标准化正态分布序列Z₀；

2)以步骤1）中获得的标准化正态分布序列Z₀与高度序列的前四阶矩统计特征作为输入，利用Johnson或Pearson转换***进行非高斯变换，获得非高斯序列Z₁；

3)利用步骤2）生成的非高斯序列Z₁，对其进行平移与缩放，获得满足均值与标准偏差高度统计特征非高斯序列Z，同时，计算其偏斜度与峰度高度统计特征，判断其精度是否满足要求；如满足精度要求，则完成非高斯粗糙表面的数字化模拟，否则，继续执行后面步骤；

4)以步骤3）获得的非高斯序列Z₁，利用反傅里叶变换生成相角序列Φ，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)；

5)以步骤5）获得序列A与权利要求1所获得的传递函数H(m,n)，通过求序列A与H的点积运算，并对点积运算结果进行反傅里叶变换，获得新的高斯表面的高度序列z₀(m,n)；

$z_0(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{s=0}{\overset{n-1}{\mathrm{\Σ}}}(A.*H)e^{2\mathrm{\πi}(\mathrm{kr}/m+\mathrm{ls}/n)}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

6)以步骤5）生成的高斯高度序列z₀，利用高斯序列的标准化方法，获得满足均值与标准偏差分别为0与1的新标准化正态分布序列Z₀；

7)重复步骤2）～3），直到非高斯序列Z满足偏斜度与峰度的精度要求为止。

所述对z₀高度方向进行平移与缩放的方法，具体包含以下步骤：

1)对高度序列z₀进行平移，使其均值等于给定μ；

z₀=z₀-mean2(z₀)+μ

2)对经平移后的高度序列z₀在高度方向进行缩放，获得标准偏差为σ的高斯高度序列z。

z=σ·z₀/std2(z₀)

所述高斯序列的标准化方法，具体包含以下步骤：

1)对高度序列z或z₀进行平移，使其均值等于0；

z₀=z₀-mean2(z₀)或z＝z-mean2(z)

2)对经平移后的高度序列z或z₀在高度方向进行缩放，获得标准偏差1的标准高斯序列。

Z₀=σ·z₀/std2(z₀)或Z₀=σ·z/std2(z)

所述的利用Johnson或Pearson转换***进行非高斯变换，具体步骤如下：

1)优先采用Johnson非高斯转换***；

对数正态拟合方式S_L：Z₁=ξ+λexp[(Z₀-γ)δ]  (ξ<η)

有界***拟合方式S_B：z=ξ+λe^(η-γ)δ/[1+e^(η-γ)/δ]  (ξ<η<ξ+λ)

无界***拟合方式S_U：z=ξ+λsinh[(η-γ)/δ]

γ、δ为形状参数，λ、ξ分别为比例参数与位置参数，γ、δ、λ与ξ均由均值、标准偏差、偏斜度与峰度等高度统计特征输入参数确定。

2)如果有界***S_B拟合方式无法收敛，则采用Pearson非高斯转换***，以提高非高斯转换***的精度：

利用概率密度函数P满足以下微分方程，根据均值、方差、偏斜度与峰度等高度统计参数获得概率密度函数与高度分布函数。

$\frac{d\left(\log p\left(x\right)\right)}{\mathrm{dx}}=-\frac{a+x}{c_0+c_{1}x+c_2{x}^2}$

式中a、c₀、c₁与c₂是由均值、标准偏差、偏斜度与峰度等高度统计特征输入参数确定。

所述的利用反傅里叶变换生成相角序列Φ，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)，具体包含以下步骤：

1)以权利要求3所获得的非高斯序列Z₁，通过反傅里叶变换，生成新的白噪声相角序列Φ(m,n)；

$\mathrm{\&Phi;}(k+1,l+1)={\tan }^{-1}\left(\frac{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{Z}_1(r+1,s+1)\sin (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{Z}_1(r+1,s+1)\cos (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}\right)\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

2)以步骤1）获得的相角序列Φ(m,n)，计算相角相关序列exp(iΦ)，并对其进行反傅里叶变换，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)；

$\mathrm{\&eta;}(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\exp (\mathrm{i\&Phi;}+\frac{2\mathrm{\&pi;ikr}}{m}+\frac{2\mathrm{\&pi;ils}}{n})\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

3)更新高斯粗糙表面数字化模拟中的白噪声序列η(m,n)的傅里叶变换序列A(m,n)。

$A(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\&ap;\mathrm{mn}\&CenterDot;e^{\mathrm{i\&Phi;}}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

本发明提出一种任意3D粗糙表面的数字化模拟方法，从高斯滤波、高精度非高斯变换、变相角循环迭代技术三个方面入手，克服了现有技术瓶颈，提高了3D微观粗糙表面的模拟精度和模拟效率，在一定程度上保证了整个偏斜度-峰度平面对应高度统计特征的任意非高斯表面的数字化模拟的可行性、准确性与高效性。可为结合面的接触性能研究提供高精度统计特征的形貌模型。

附图说明

图1各向同性高斯表面及其高度分布（β_x=β_y=5，μ=0，σ=1.6，点数128×128）

图2各向异性高斯表面及其高度分布（β_x=5,β_y=50，μ=0，σ＝1.6，点数512×512）

图3各向同性非高斯粗糙表面及其高度分布（β_x=β_y=5，μ=0，σ=1.6，S_kz=-1，K_uz=4，点数128×128）

图4各向同性粗糙表面的x、y方向平均自相关函数与理论自相关函数的对比：（a）x方向自相关函数（β_x=5），（b）y方向自相关函数（β_x=5）

图5各向异性非高斯粗糙表面及其高度分布（β_x=5,β_y=50,μ=0,σ＝1.6,Skz=1,K_uz=5点数512×512）

图6各向异性粗糙表面的x、y方向平均自相关函数与理论自相关函数的对比：（a）x方向自相关函数（β_x=5），（b）y方向自相关函数（β_y=50）

具体实施方式

以满足一定的自相关函数、均值、标准偏差、偏斜度与峰度等统计特征的3D高斯或非高斯粗糙表面的数字化模拟为例，对本发明的任意3D粗糙表面的数字化模拟方法进行说明。

1.高斯粗糙表面的数字化模拟：

1)利用给定自相关函数生成相关函数离散矩阵R(m,n)

a)假定自相关函数为f(x,y)（σ为高斯序列的标准偏差，β_x、β_y分别为X、Y方向的自相关长度）；

$f(x,y)={\mathrm{\&sigma;}}^2\exp (-2.3\sqrt{{\left(\frac{x}{{\mathrm{\&beta;}}_x}\right)}^2+{\left(\frac{y}{{\mathrm{\&beta;}}_y}\right)}^2})$

b)指定X、Y方向的离散间距均为1μm，并在-m/2≤x≤m/2与-n/2≤y≤n/2区域内对自相关函数进行离散，获得自相关函数的离散矩阵R(m+1,n+1)；

$R(k+\frac{m}{2}+1,l+\frac{n}{2}+1)=f(k,l),k=-\frac{m}{2},-\frac{m}{2}+1...\frac{m}{2},l=-\frac{n}{2},-\frac{n}{2}+1...\frac{n}{2}$

c)所得R(m+1,n+1)序列中的第m/2+1，m/2+2行重复，删除其中一行；第n/2+1与n/2+2列重复，删除其中一列；获得与模拟表面高度序列具有相同行列数的自相关函数矩阵R(m,n)。

2)根据自相关函数求高斯表面的功率谱密度与传递函数；

a)对自相关函数离散矩阵进行傅里叶变换，并求解高斯粗糙表面的功率谱密度P；

$P(1+k,1+l)=\left|\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}R(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\right|\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

b)根据白噪声的功率谱密度为常数C，假定C=1，则传递函数为；

$H=\sqrt{P}$

3)白噪声序列的生成与傅里叶变换。

a)利用白噪声生成器randn(m,n)产生白噪声序列η(m,n)；

b)求白噪声序列η的傅里叶变换A(m,n)：

$A(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

4)利用频域点乘（A.*H）取反傅里叶变换的方法获得高斯表面的高度序列；

$z_0(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}(A.*H)e^{2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

5)根据给定高度均值μ、标准偏差σ与获得序列z₀的实际均值mean2(z₀)（mean2(z₀)求序列均值的函数）与标准偏差std2(z)（std2为求序列标准偏差函数）的关系，对高斯高度序列z₀进行高度方向平移与缩放，消除白噪声功率谱密度长度C不等于1对高度序列统计特征的影响，获得满足均值与标准偏差分布为μ与σ的高斯高度序列z，完成高斯粗糙表面的数字化模拟。

z₀=z₀-mean2(z₀)+μ  z=σ·z₀/std2(z₀)

利用上述给定的自相关函数，x、y方向的自相关函数均取5μm，模拟表面的x，y方向取点为128×128、表面的均值为0，标准偏差为1.6μm，模拟获得的各向同性高斯表面形貌及其高度分布如图1所示，其前四阶矩高度统计特征误差如表1所示。x、y方向的自相关函数分别取5μm、50μm，模拟获得的高斯表面形貌表现出各向异性，具有一定的纹理。模拟获得的各向同性高斯表面形貌及其高度分布如图2所示，其前四阶矩高度统计特征误差如表2所示。

表1 各向同性高斯粗糙表面（β_x=β_y=5，μ=0，σ＝1.6，点数128×128）的高度统计前四阶矩仿真精度

表2 各向异性高斯粗糙表面（β_x=5，β_y=50，μ=0，σ＝1.6，点数512×512）的高度统计前四阶矩仿真精度计算

2.非高斯粗糙表面的数字化模拟：

6)对上述模拟的高斯高度序列z平移与缩放，获得均值与标准偏差分别为0与1的标准化正态分布序列Z₀；

z＝z-mean2(z)  Z₀=z/std2(z)

7)采用Johnson或Pearson非高斯转换***进行非高斯转换，生成非高斯序列Z₁；

a)优先采用Johnson非高斯转换***；

对数正态拟合方式S_L：Z₁=ξ+λexp[(Z₀-γ)δ]  (ξ<η)

有界***拟合方式S_B： $Z_1=\mathrm{\&xi;}+{\mathrm{\&lambda;e}}^{(Z_0-\mathrm{\&gamma;})/\mathrm{\&delta;}}/[1+e^{(Z_0-\mathrm{\&gamma;})/\mathrm{\&delta;}}],(\mathrm{\&xi;}<\mathrm{\&eta;}<\mathrm{\&xi;}+\mathrm{\&lambda;})$

无界***拟合方式S_U：Z₁=ξ+λsinh[(Z₀-γ)/δ]

γ、δ为形状参数，λ、ξ分别为比例参数与位置参数，γ、δ、λ与ξ均由均值、标准偏差、偏斜度与峰度等高度统计特征输入参数确定。

b)如果有界***S_B拟合方式无法收敛，则采用Pearson非高斯转换***，以提高非高斯转换***的精度：

利用概率密度函数P满足以下微分方程，根据均值、方差、偏斜度与峰度等高度统计参数获得概率密度函数与高度分布函数。

$\frac{d\left(\log p\left(x\right)\right)}{\mathrm{dx}}=-\frac{a+x}{c_0+c_{1}x+c_2{x}^2}$

式中a、c₀、c₁与c₂是由均值、标准偏差、偏斜度与峰度等高度统计特征输入参数确定。

8)非高斯高度序列偏斜度与峰度精度控制与循环迭代策略：

a)对非高斯高度序列Z₁进行平移与缩放，获得均值与标准偏差分别为μ与σ的非高斯高度序列Z；

Z₁=Z₁-mean2(Z₁)+μ  Z=σ·Z₁/std2(Z₁)

b)计算非高斯序列Z的偏斜度S_kz与峰度K_uz，并判断其精度能否满足要求，若满足精度要求，则完成非高斯粗糙表面的模拟，否则继续执行后面的步骤；

$S_{\mathrm{kz}}=\frac{\sqrt{\mathrm{mn}}\&CenterDot;\underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}{(Z-\stackrel{\&OverBar;}{Z})}^3}{{\left(\underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}{(Z-\stackrel{\&OverBar;}{Z})}^2\right)}^{1.5}}$ $K_{\mathrm{uz}}=\frac{\mathrm{mn}\&CenterDot;\underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}{(Z-\stackrel{\&OverBar;}{Z})}^4}{{\left(\underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}{(Z-\stackrel{\&OverBar;}{Z})}^2\right)}^2}$

c)利用非高斯序列Z₁的快速反傅里叶变换，更新的相角序列Φ(m,n)；

$\mathrm{\&Phi;}(k+1,l+1)={\tan }^{-1}\left(\frac{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{Z}_1(r+1,s+1)\sin (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{Z}_1(r+1,s+1)\cos (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}\right)\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

d)利用相角序列，计算相角相关序列exp(iΦ)，对其进行快速反傅里叶变换更新高斯粗糙表面数字化模拟中的白噪声序列η，则更新的白噪声序列的傅里叶变换A约等于mn·exp(iΦ)；

$\mathrm{\&eta;}(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\exp (\mathrm{i\&Phi;}+\frac{2\mathrm{\&pi;ikr}}{m}+\frac{2\mathrm{\&pi;ils}}{n})\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

$A(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\&ap;\mathrm{mn}\&CenterDot;e^{\mathrm{i\&Phi;}}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

e)重复步骤4)～8)，直到非高斯序列Z满足偏斜度与峰度的精度要求。

利用模拟的各向同性高斯粗糙表面高度序列z，其x、y方向的自相关函数均取5μm，x、y方向取点为128×128、表面的均值为0，标准偏差为1.6μm，通过非高斯模拟获得偏斜度S_kz与峰度K_uz分别为-1与4的各向同性非高斯序列，对应的非高斯表面形貌及其高度分布如图3所示。产生非高斯粗糙表面的迭代循环次数为2，模拟时间小于0.5秒。在高度统计特征方面：其前四阶矩高度统计特征误差如表3所示，其误差小于1.00%；在自相关函数统计特征方面：如图4所示，产生的高斯粗糙粗糙表面与非高斯粗糙表面对应的x方向与y方向平均自相关函数（已经标准化，最大值为1）接近重合，且与相应理论自相关函数基本吻合，在下部稍有偏离，可以通过增加模拟的点数来提高其吻合程度。

利用模拟的各向异性高斯粗糙表面高度序列z，其x、y方向的自相关函数分别为5μm与50μm，x、y方向取点为512×512、表面的均值为0，标准偏差为1.6μm，通过非高斯模拟获得偏斜度S_kz与峰度K_uz分别为1与5的各向异性非高斯序列，对应的非高斯表面形貌及其高度分布如图5所示。产生非高斯粗糙表面的迭代循环次数为3，模拟时间小于0.5秒。在高度统计特征方面：其前四阶矩高度统计特征误差如表4所示，其误差小于1.00%；在自相关函数统计特征方面：如图6所示，由于x方向模拟点数512（反映模拟长度，原来点数为128）与x方向自相关长度5μm的比值明显增大，模拟产生的高斯粗糙粗糙表面与非高斯粗糙表面的x方向平均自相关函数（已经标准化，最大值为1）重合，几乎与相应理论自相关函数曲线重合；由于y方向模拟点数512（反映模拟长度，原来点数为128）与y方向自相关长度50μm的比值稍有减小，导致在自相关函数值较小的情况下，模拟产生的高斯粗糙粗糙表面与非高斯粗糙表面的y方向平均自相关函数曲线与相应的理论自相关曲线有一定的偏离，这是由于模拟点数反映模拟长度，某方向的模拟长度与相应自相关长度比值的减小，会导致该方向的相似波数目变少，导致该方向的统计特征偏离。如要提高自相关函数的统计精度，可以通过增加点数与减小自相关长度。

表3 非高斯表面（β_x＝β_y=5，μ=0，σ=1.6，S_kz=-1，K_uz=4，点数128×128）的高度统计前四阶矩仿真精度计算

表4 非高斯表面（β_x=5,β_y=50,μ=0,σ=1.6,S_kz=1,K_uz=5点数512×512）的高度统计前四阶矩仿真精度计算

Claims (6)

1.一种3D任意粗糙表面的数字化模拟方法，其特征在于，包含下述步骤：

1)当任意粗糙表面的偏斜度值等于0且峰度值约等于3时，其表面形貌符合高斯分布，否则，其表面形貌符合非高斯分布；

2)如果步骤1)中粗糙表面形貌为高斯分布，利用三维高斯粗糙表面的模拟方法获得粗糙表面高度序列z，完成粗糙表面的数字化模拟；

3)如果步骤1)中粗糙表面形貌为非高斯分布，利用三维非高斯粗糙表面的模拟方法获得粗糙表面高度序列Z，完成粗糙表面的数字化模拟；

其中，步骤2)中三维高斯粗糙表面的模拟方法，具体包含以下步骤：

(1)通过给定的自相关函数生成离散的自相关函数矩阵R(m,n)，其中m、n分别表示x、y方向取点的数目；

(2)对离散的自相关函数矩阵R(m,n)进行快速傅里叶变换，获得高斯粗糙表面的功率谱密度P(m,n)；

$P(1+k,1+l)=\left|\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}R(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\right|\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

(3)假设白噪声的功率谱密度常数C等于1，利用获得的功率谱密度，计算传递函数H(m,n)；

$H=\sqrt{P}$

(4)利用白噪声生成器，获得白噪声序列η(m,n)，利用反傅里叶变换生成相角序列Φ；

$\mathrm{\&Phi;}(k+1,l+1)={\tan }^{-1}\left(\frac{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1s+1)\sin (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1,s+1)\cos (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}\right)\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

(5)以步骤(4)获得的相角序列Φ，计算相角相关序列exp(iΦ)，并对其进行反傅里叶变换，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)；

$\mathrm{\&eta;}(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\exp (\mathrm{i\&Phi;}+\frac{2\mathrm{\&pi;ikr}}{m}+\frac{2\mathrm{\&pi;ils}}{n})\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

$A(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\&ap;\mathrm{mn}\&CenterDot;e^{\mathrm{i\&Phi;}}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

(6)以步骤(3)获得的传递函数H(m,n)，步骤5)获得的白噪声傅里叶变换序列A(m,n)，求序列H与序列A的点积运算，并对点积运算结果进行反傅里叶变换，获得高斯表面的高度序列z₀(m,n)；

$z_0(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}(A.*H)e^{2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

(7)以步骤(6)生成的高斯高度序列z₀，对其进行高度方向平移与缩放，获得满足均值与标准偏差的高斯高度序列z，完成高斯粗糙表面的数字化模拟。

2.根据权利要求1所述的一种3D任意粗糙表面的数字化模拟方法，所述对z₀高度方向进行平移与缩放的方法，具体包含以下步骤：

1)对高度序列z₀进行平移，使其均值等于给定μ；

z₀＝z₀-mean2(z₀)+μ

其中mean2(z₀)为求序列z₀均值的函数，μ为高度均值；

2)对经平移后的高度序列z₀在高度方向进行缩放，获得标准偏差为σ的高斯高度序

列z，

z＝σ·z₀/std2(z₀)

其中std2(z₀)为求序列z₀的标准偏差函数。

3.根据权利要求1所述的一种3D任意粗糙表面的数字化模拟方法，所述的三维非高斯粗糙表面的模拟方法，具体包含以下步骤：

1)以三维高斯粗糙表面的模拟方法中所获得的高斯高度序列z，利用高斯序列的标准化方法，获得标准化正态分布序列Z₀；

2)以步骤1)中获得的标准化正态分布序列Z₀与高度序列的前四阶矩统计特征作为输入，利用Johnson或Pearson转换***进行非高斯变换，获得非高斯序列Z₁；

3)利用步骤2)生成的非高斯序列Z₁，对其进行平移与缩放，获得满足均值与标准偏差高度统计特征非高斯序列Z，同时，计算其偏斜度与峰度高度统计特征，判断其精度是否满足要求；如满足精度要求，则完成非高斯粗糙表面的数字化模拟，否则，继续执行后面步骤；

4)以步骤3)获得的非高斯序列Z₁，利用反傅里叶变换生成相角序列Φ，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)；

5)以步骤5)获得序列A与权利要求1所获得的传递函数H(m,n)，通过求序列A与H的点积运算，并对点积运算结果进行反傅里叶变换，获得新的高斯表面的高度序列z₀(m,n)；

$z_0(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}(A.*H)e^{2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

6)以步骤5)生成的高斯高度序列z₀，利用高斯序列的标准化方法，获得满足均值与标准偏差分别为0与1的新标准化正态分布序列Z₀；

7)重复步骤2)～3)，直到非高斯序列Z满足偏斜度与峰度的精度要求为止。

4.根据权利要求3所述的一种3D任意粗糙表面的数字化模拟方法，所述高斯序列的标准化方法，具体包含以下步骤：

1)对高度序列z₀进行平移，使其均值等于0；

z₀＝z₀-mean2(z₀)

2)对经平移后的高度序列z₀在高度方向进行缩放，获得标准偏差1的标准高斯序列，Z₀＝z₀/std2(z₀)。

5.根据权利要求3所述的一种3D任意粗糙表面的数字化模拟方法，所述的利用Johnson或Pearson转换***进行非高斯变换，具体步骤如下：

1)优先采用Johnson非高斯转换***；

对数正态拟合方式S_L：Z₁＝ξ+λexp[(Z₀-γ)/δ]   (ξ<η)

有界***拟合方式S_B：z＝ξ+λe^(η-γ)/δ/[1+e^(η-γ)/δ]   (ξ<η<ξ+λ)

无界***拟合方式S_U：z＝ξ+λsinh[(η-γ)/δ]

γ、δ为形状参数，λ、ξ分别为比例参数与位置参数，γ、δ、λ与ξ均由均值、标准偏差、偏斜度与峰度等高度统计特征输入参数确定；

2)如果有界***S_B拟合方式无法收敛，则采用Pearson非高斯转换***，以提高非高斯转换***的精度：

利用概率密度函数P满足以下微分方程，根据均值、方差、偏斜度与峰度等高度统计参数获得概率密度函数与高度分布函数：

$\frac{d\left(\log p\left(x\right)\right)}{\mathrm{dx}}=-\frac{a+x}{c_0+c_{1}x+c_2{x}^2}$

式中a、c₀、c₁与c₂是由均值、标准偏差、偏斜度与峰度等高度统计特征输入参数确定。

6.根据权利要求3所述的一种3D任意粗糙表面的数字化模拟方法，所述的利用反傅里叶变换生成相角序列Φ，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)，具体包含以下步骤：

1)以三维非高斯粗糙表面的模拟方法中所获得的非高斯序列Z₁，通过反傅里叶变换，生成新的白噪声相角序列Φ(m,n)；

$\mathrm{\&Phi;}(k+1,l+1)={\tan }^{-1}\left(\frac{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{Z}_1(r+1,s+1)\sin (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}{-\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{Z}_1(r+1,s+1)\cos (2\mathrm{\&pi;kr}/m+2\mathrm{\&pi;ls}/n)}\right)\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

2)以步骤1)获得的相角序列Φ(m,n)，计算相角相关序列exp(iΦ)，并对其进行反傅里叶变换，更新白噪声序列η(m,n)及其傅里叶变换A(m,n)；

$\mathrm{\&eta;}(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\exp (\mathrm{i\&Phi;}+\frac{2\mathrm{\&pi;ikr}}{m}+\frac{2\mathrm{\&pi;ils}}{n})\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}$

3)更新高斯粗糙表面数字化模拟中的白噪声序列η(m,n)的傅里叶变换序列A(m,n)；

$A(k+1,l+1)=\underset{r=0}{\overset{m-1}{\mathrm{\&Sigma;}}}\underset{s=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\mathrm{\&eta;}(r+1,s+1)e^{-2\mathrm{\&pi;i}(\mathrm{kr}/m+\mathrm{ls}/n)}\&ap;\mathrm{mn}\&CenterDot;e^{\mathrm{i\&Phi;}}\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}.$