CN102609560B - Digitalized simulation method of any three-dimensional (3D) rough surface - Google Patents

Abstract

The invention discloses a digitalized simulation method of any three-dimensional (3D) rough surface, which improves simulation efficiency and simulation accuracy of statistical parameters. Random initial phase angle sequences are obtained through inverse Fourier transform of white noise sequences, and novel white noise and Fourier transform are obtained by utilizing the initial phase angle sequences. Processing including dispersion and Fourier transform and the like is performed through appointed autocorrelation functions to obtain power spectral density and transfer functions of Gaussian rough surface height sequences. The simulation of the Gaussian rough surface height sequences is finished by utilizing frequency domain dot product and a method for obtaining the inverse Fourier transform. On the basis, non-Gaussian rough surfaces are generated through height distribution statistical characteristic parameters including appointed skewness degree, kurtosis and the like and by utilizing combined Pearson and Johnson non-Gaussian transition systems. If the simulation accuracy of the skewness degree and the kurtosis is below standard, the phase angle sequences and the Fourier transform of the white noise are updated to re-perform Gaussian filtering and non-Gaussian transition till the given accuracy requirement is met.

Description

The digitized simulation method of any rough surface of a kind of 3D

Technical field

The present invention relates to the digital simulation method of any rough surface of a kind of 3D microcosmic, more particularly, the present invention, according to Fourth-order moment statistical nature and horizontal direction autocorrelation function statistical law before the height distribution of connection table millet cake, the rough surface with identical statistical nature can be built efficiently and accurately.

Background technology

At a high speed, improving constantly of requiring of the development of precision processing technology and production reliability, promote the fast development in the fields such as the processing and manufacturing of rough surface, detection, digitized simulation and performance prediction.

In recent years, research both at home and abroad in the digital simulation field of rough surface, mainly concentrate on the rough surface aspect meeting Gaussian distribution for roughness height, in the digitized simulation of non-gaussian rough surface, can be successfully non-gaussian data by the data transformations of Gaussian distribution, occur based on non-gaussian converting systems such as Johnson, Pearson.But in the middle of the digitized simulation of non-gaussian rough surface, non-gaussian sequence after non-gaussian conversion also needs to carry out filtering, introduce new error to statistical nature parameters such as the measure of skewness of simulated series, kurtosis, have also been changed the simulated domain of the correspondence of measure of skewness-kurtosis simultaneously.Thus cause the rough surface of some statistical nature not generate corresponding non-gaussian sequence by converting system conversions such as Johnson, Pearson to filtering, i.e. enable transition, also can introduce new error at filtering stage, affect the precision of the non-gaussian surface digitizing simulation of filtering.

Consider the impact of filtering (filter factor is h (k, l)), the measure of skewness S of the non-gaussian sequence rough surface generated in the non-gaussian converting systems such as Johnson, Pearson _{k η}, kurtosis K _{u η}with the measure of skewness S of filtered rough surface _kZ, kurtosis K _uZthere is following relation:

$\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.$

Even if target deflection degree and kurtosis meet non-gaussian change necessary condition: but after considering the influencing each other of filter factor and autocorrelation function, S _{k η}, K _{u η}likely no longer meet the necessary condition of non-gaussian conversion: thus cause not generating non-gaussian sequence by non-gaussian converting systems such as Johnson or Pearson; I.e. enable transition, also can introduce new error at filtering stage, affects the precision of the non-gaussian surface digitizing simulation of filtering.

In reality, different job operations can produce the rough surface with different statistical nature, the carrier that the digitized simulation model of rough surface is analyzed as interface properties, and accuracy of its simulation directly affects the contact performance prediction of corresponding faying face.The sequence of phase angles generation white noise that this research and utilization is certain, combined with the non-gaussian converting system of Gaussian surface high degree of sequence by the gaussian filtering of white noise, form the non-gaussian high degree of sequence meeting height distribution statistics feature and autocorrelation function, the digitized simulation method of Gauss or non-gaussian rough surface is that interface properties research is laid a good foundation.

Summary of the invention

Fundamental purpose of the present invention is to provide a kind of digitized simulation method that high-efficiency high-accuracy realizes any rough surface, decrease the filtering error link in rough surface digitized simulation, before meeting autocorrelation function, highly distribution the statistical nature such as Fourth-order moment prerequisite under, improve and simulate the precision of rough surface height distribution statistics feature and the efficiency of digitized simulation.

The present invention is achieved through the following technical solutions, and key step comprises:

A digitized simulation method for any rough surface of 3D, is characterized in that, comprises following step:

1) when the deflection angle value of any rough surface approximates 3 and kurtosis value approximates 0, its surface topography meets Gaussian distribution, otherwise its surface topography meets non-gaussian distribution;

2) if rough surface morphology is Gaussian distribution in step 1), utilize the analogy method of three-dimensional Gaussian rough surface to obtain rough surface high degree of sequence z, complete the digitized simulation of rough surface;

3) if rough surface morphology is non-gaussian distribution in step 1), utilize the analogy method of three-dimensional non-gaussian rough surface to obtain rough surface high degree of sequence z, complete the digitized simulation of rough surface.

The analogy method of described three-dimensional Gaussian rough surface, specifically comprises following steps:

1) generate discrete autocorrelation function matrix R (m, n) by given autocorrelation function, wherein m, n represent that number is a little got in x, y direction respectively;

2) Fast Fourier Transform (FFT) is carried out to discrete autocorrelation function matrix R (m, n), obtain the power spectrum density P (m, n) of Gauss's rough surface;

$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) suppose that the power spectrum density constant C of white noise equals 1, utilize the power spectrum density obtained, calculation of transfer function H (m, n);

$H=\sqrt{P}$

4) utilize white noise maker, obtain white noise sequence η (m, n), utilize inversefouriertransform to generate sequence of phase angles Φ;

$\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) with the sequence of phase angles Φ that step 4) obtains, calculate phase angle correlated series exp (i Φ), and inversefouriertransform is carried out to it, upgrade white noise sequence η (m, n) and Fourier transform A (m, n) thereof;

$\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) with the transfer function H (m, n) that step 3) obtains, the white noise Fourier transform sequence A (m that step 5) obtains, n), ask the dot-product operation of sequence H and sequence A, and inversefouriertransform is carried out to dot-product operation result, obtain the high degree of sequence z of Gaussian surface ₀(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) with Gauss's high degree of sequence z that step 6) generates ₀, short transverse translation and convergent-divergent are carried out to it, obtain the Gauss's high degree of sequence z meeting average and standard deviation, complete the digitized simulation of Gauss's rough surface.

The analogy method of described three-dimensional non-gaussian rough surface, specifically comprises following steps:

1) with the Gauss's high degree of sequence z obtained, utilize the standardized method of gaussian sequence, carry out translation and convergent-divergent in short transverse, obtain standardized normal distribution series Z ₀;

2) with the standardized normal distribution series Z obtained in step 1) ₀with Fourth-order moment statistical nature before high degree of sequence as input, utilize Johnson or Pearson converting system to carry out non-gaussian conversion, obtain non-gaussian sequence Z ₁;

3) step 2 is utilized) the non-gaussian sequence Z that generates ₁, translation and convergent-divergent are carried out to it, obtain and meet average and standard deviation height statistical nature non-gaussian sequence Z, meanwhile, calculate its measure of skewness and kurtosis height statistical nature, judge whether its precision meets the demands; As met accuracy requirement, then complete the digitized simulation of non-gaussian rough surface, otherwise, continue to perform later step;

4) with the non-gaussian sequence Z that step 3) obtains ₁, utilize inversefouriertransform to generate sequence of phase angles Φ, upgrade white noise sequence η (m, n) and Fourier transform A (m, n) thereof;

5) with the transfer function H (m, n) that step 5) acquisition sequence A and claim 1 obtain, by asking the dot-product operation of sequence A and H, and inversefouriertransform being carried out to dot-product operation result, obtaining the high degree of sequence z of new Gaussian surface ₀(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) with Gauss's high degree of sequence z that step 5) generates ₀, utilize the standardized method of gaussian sequence, obtain and meet the new standardized normal distribution sequence Z that average and standard deviation are respectively 0 and 1 ₀;

7) step 2 is repeated) ~ 3), until non-gaussian sequence Z meets the accuracy requirement of measure of skewness and kurtosis.

Described to z ₀short transverse carries out the method for translation and convergent-divergent, specifically comprises following steps:

1) to high degree of sequence z ₀carry out translation, make its average equal given μ;

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

2) to the high degree of sequence z after translation ₀carry out convergent-divergent in short transverse, obtain Gauss's high degree of sequence z that standard deviation is σ.

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

The standardized method of described gaussian sequence, specifically comprises following steps:

1) to high degree of sequence z or z ₀carry out translation, make its average equal 0;

Z ₀=z ₀-mean2 (z ₀) or z=z-mean2 (z)

2) to high degree of sequence z or z after translation ₀carry out convergent-divergent in short transverse, obtain the standard gaussian sequence of standard deviation 1.

Z ₀=σ z ₀/ std2 (z ₀) or Z ₀=σ z/std2 (z)

Described Johnson or the Pearson converting system that utilizes carries out non-gaussian conversion, and concrete steps are as follows:

1) preferential employing Johnson non-gaussian converting system;

Lognormality fit approach S _l: Z ₁=ξ+λ exp [(Z ₀-γ) δ] (ξ < η)

Bounded system fit approach S _b: z=ξ+λ e ^{(η-γ) δ}/ [1+e ^(η-γ)/δ] (ξ < η < ξ+λ)

Boundless system fit approach S _u: z=ξ+λ sinh [(η-γ)/δ]

γ, δ are form parameter, and λ, ξ are respectively scale parameter and location parameter, and γ, δ, λ and ξ determine by average, standard deviation, measure of skewness and kurtosis equal altitudes statistical nature input parameter.

2) if bounded system S _bfit approach cannot be restrained, then adopt Pearson non-gaussian converting system, to improve the precision of non-gaussian converting system:

Utilize probability density function P to meet the following differential equation, obtain probability density function and height distribution function according to average, variance, measure of skewness and kurtosis equal altitudes statistical parameter.

$\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 in formula ₀, c ₁with c ₂determined by average, standard deviation, measure of skewness and kurtosis equal altitudes statistical nature input parameter.

The described inversefouriertransform that utilizes generates sequence of phase angles Φ, upgrades white noise sequence η (m, n) and Fourier transform A (m, n) thereof, specifically comprises following steps:

1) with the non-gaussian sequence Z that claim 3 obtains ₁, by inversefouriertransform, generate new white noise sequence of phase angles Φ (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) with the sequence of phase angles Φ (m, n) that step 1) obtains, phase angle correlated series exp (i Φ) is calculated, and inversefouriertransform is carried out to it, upgrade white noise sequence η (m, n) and Fourier transform A (m, n) thereof;

$\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) the Fourier transform sequence A (m, n) of the white noise sequence η (m, n) in Gauss's rough surface digitized simulation is upgraded.

$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}$

The present invention proposes a kind of digitized simulation method of any 3D rough surface, start with from gaussian filtering, the conversion of high precision non-gaussian, covert loop iteration technology three aspects, angle, overcome prior art bottleneck, improve simulation precision and the simulation precision of 3D micro-roughened surface, ensure that the feasibility of the digitized simulation on any non-gaussian surface of whole measure of skewness-kurtosis plane respective heights statistical nature, accuracy and high efficiency to a certain extent.The contact behaviours that can be faying face provides the pattern model of high precision statistical nature.

Accompanying drawing explanation

Fig. 1 isotropy Gaussian surface and the (β that highly distributes thereof _x=β _y=5, μ=0, σ=1.6, count 128 × 128)

Fig. 2 anisotropic Gaussian surface and the (β that highly distributes thereof _x=5, β _y=50, μ=0, σ=1.6, count 512 × 512)

Fig. 3 isotropy non-gaussian rough surface and the (β that highly distributes thereof _x=β _y=5, μ=0, σ=1.6, S _kz=-1, K _uz=4, count 128 × 128)

X, y orientation average autocorrelation function of Fig. 4 isotropy rough surface and the contrast of theoretical autocorrelation function: (a) x direction autocorrelation function (β _x=5), (b) y direction autocorrelation function (β _x=5)

Fig. 5 anisotropy non-gaussian rough surface and the (β that highly distributes thereof _x=5, β _y=50, μ=0, σ=1.6, Skz=1, K _uz=5 count 512 × 512)

X, y orientation average autocorrelation function of Fig. 6 anisotropic rough surface and the contrast of theoretical autocorrelation function: (a) x direction autocorrelation function (β _x=5), (b) y direction autocorrelation function (β _y=50)

Embodiment

To meet the 3D Gauss of statistical nature or the digitized simulations of non-gaussian rough surface such as certain autocorrelation function, average, standard deviation, measure of skewness and kurtosis, the digitized simulation method of any 3D rough surface of the present invention is described.

1. the digitized simulation of Gauss's rough surface:

1) given autocorrelation function is utilized to generate related function discrete matrix R (m, n)

A) suppose that autocorrelation function is that (σ is the standard deviation of gaussian sequence to f (x, y), β _x, β _ybe respectively the auto-correlation length of X, Y-direction);

$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) the discrete spacing of appointment X, Y-direction is 1 μm, and carries out discrete to autocorrelation function in-m/2≤x≤m/2 with-n/2≤y≤n/2 region, obtains the discrete matrix R (m+1, n+1) of autocorrelation function;

$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) the capable repetition of the m/2+1 in gained R (m+1, n+1) sequence, m/2+2, deletes wherein a line; N-th/2+1 and n/2+2 column weight is multiple, deletes wherein row; Obtain autocorrelation function matrix R (m, n) with simulation surface elevation sequence with identical ranks number.

2) power spectrum density and the transport function of Gaussian surface is asked according to autocorrelation function;

A) Fourier transform is carried out to autocorrelation function discrete matrix, and solve the power spectrum density P of Gauss's rough surface;

$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) be constant C according to the power spectrum density of white noise, assuming that C=1, then transport function is;

$H=\sqrt{P}$

3) generation of white noise sequence and Fourier transform.

A) white noise maker randn (m, n) is utilized to produce white noise sequence η (m, n);

B) the Fourier transform A (m, n) of white noise sequence η is asked:

$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) method of frequency domain dot product (A.*H) negate Fourier transform is utilized to obtain the high degree of sequence of Gaussian surface;

$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) according to assigned altitute average μ, standard deviation and acquisition sequence z ₀actual average mean2 (z ₀) (mean2 (z ₀) ask the function of serial mean) and with the relation of standard deviation std2 (z) (std2 is for asking sequence criteria departure function), to Gauss's high degree of sequence z ₀carry out short transverse translation and convergent-divergent, eliminate Power Spectrum of White Noise density length C and be not equal to 1 impact on high degree of sequence statistical nature, obtain and meet Gauss's high degree of sequence z that average and distribution of standard deviation are μ and σ, complete the digitized simulation of Gauss's rough surface.

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

Utilize above-mentioned given autocorrelation function, the autocorrelation function in x, y direction all gets 5 μm, the x on simulation surface, it is a little 128 × 128 that y gets in direction, the average on surface is 0, standard deviation is 1.6 μm, as shown in Figure 1, its front Fourth-order moment height statistical nature error is as shown in table 1 for the isotropy Gaussian surface pattern that simulation obtains and highly distribution.The autocorrelation function in x, y direction gets 5 μm, 50 μm respectively, and the Gaussian surface pattern that simulation obtains shows anisotropy, has certain texture.As shown in Figure 2, its front Fourth-order moment height statistical nature error is as shown in table 2 for the isotropy Gaussian surface pattern that simulation obtains and highly distribution.

Table 1 isotropy Gauss rough surface (β _x=β _y=5, μ=0, σ=1.6, count 128 × 128) height statistics before Fourth-order moment simulation accuracy

Table 2 anisotropic Gaussian rough surface (β _x=5, β _y=50, μ=0, σ=1.6, count 512 × 512) height statistics before Fourth-order moment simulation accuracy calculate

2. the digitized simulation of non-gaussian rough surface:

6) to the z translation of Gauss's high degree of sequence and the convergent-divergent of above-mentioned simulation, acquisition average and standard deviation are respectively the standardized normal distribution series Z of 0 and 1 ₀;

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

7) adopt Johnson or Pearson non-gaussian converting system to carry out non-gaussian conversion, generate non-gaussian sequence Z ₁;

A) preferential employing Johnson non-gaussian converting system;

Lognormality fit approach S _l: Z ₁=ξ+λ exp [(Z ₀-γ) δ] (ξ < η)

Bounded system fit approach 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;})$

Boundless system fit approach S _u: Z ₁=ξ+λ sinh [(Z ₀-γ)/δ]

γ, δ are form parameter, and λ, ξ are respectively scale parameter and location parameter, and γ, δ, λ and ξ determine by average, standard deviation, measure of skewness and kurtosis equal altitudes statistical nature input parameter.

If b) bounded system S _bfit approach cannot be restrained, then adopt Pearson non-gaussian converting system, to improve the precision of non-gaussian converting system:

Utilize probability density function P to meet the following differential equation, obtain probability density function and height distribution function according to average, variance, measure of skewness and kurtosis equal altitudes statistical parameter.

$\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 in formula ₀, c ₁with c ₂determined by average, standard deviation, measure of skewness and kurtosis equal altitudes statistical nature input parameter.

8) non-gaussian high degree of sequence measure of skewness and kurtosis precision controlling and loop iteration strategy:

A) to non-gaussian high degree of sequence Z ₁carry out translation and convergent-divergent, acquisition average and standard deviation are respectively the non-gaussian high degree of sequence Z of μ and σ;

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

B) the measure of skewness S of non-gaussian sequence Z is calculated _kzwith kurtosis K _uz, and judge that can its precision meet the demands, if meet accuracy requirement, then complete the simulation of non-gaussian rough surface, otherwise continue to perform step below;

$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) non-gaussian sequence Z is utilized ₁quick inversefouriertransform, the sequence of phase angles Φ (m, n) of renewal;

$\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) sequence of phase angles is utilized, calculate phase angle correlated series exp (i Φ), carry out the white noise sequence η in quick inversefouriertransform renewal Gauss rough surface digitized simulation to it, then the Fourier transform A of the white noise sequence upgraded approximates mnexp (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) step 4 is repeated) ~ 8), until non-gaussian sequence Z meets the accuracy requirement of measure of skewness and kurtosis.

Utilize the isotropy Gauss rough surface high degree of sequence z of simulation, the autocorrelation function in its x, y direction all gets 5 μm, and it is a little 128 × 128 that x, y get in direction, the average on surface is 0, and standard deviation is 1.6 μm, obtains measure of skewness S by non-gaussian simulation _kzwith kurtosis K _uzbe respectively the isotropy non-gaussian sequence of-1 and 4, corresponding non-gaussian surface topography and highly distribution are as shown in Figure 3.The iterative loop number of times producing non-gaussian rough surface is 2, and simulated time is less than 0.5 second.In height statistical nature: its front Fourth-order moment height statistical nature error is as shown in table 3, and its error is less than 1.00%; In autocorrelation function statistical nature: as shown in Figure 4, the x direction that the coarse rough surface of Gauss produced is corresponding with non-gaussian rough surface and (the standardization of y orientation average autocorrelation function, maximal value is 1) close to overlapping, and it is substantially identical with corresponding theory autocorrelation function, being slightly offset in bottom, its degree of agreement can being improved by increasing counting of simulation.

Utilize the anisotropic Gaussian rough surface high degree of sequence z of simulation, the autocorrelation function in its x, y direction is respectively 5 μm and 50 μm, it is a little 512 × 512 that x, y get in direction, the average on surface is 0, and standard deviation is 1.6 μm, obtains measure of skewness S by non-gaussian simulation _kzwith kurtosis K _uzbe respectively the anisotropy non-gaussian sequence of 1 and 5, corresponding non-gaussian surface topography and highly distribution are as shown in Figure 5.The iterative loop number of times producing non-gaussian rough surface is 3, and simulated time is less than 0.5 second.In height statistical nature: its front Fourth-order moment height statistical nature error is as shown in table 4, and its error is less than 1.00%, in autocorrelation function statistical nature: as shown in Figure 6, because 512(reflection simulation length of counting is simulated in x direction, originally counted is 128) obviously increase with the ratio of x direction auto-correlation length 5 μm, x orientation average autocorrelation function (the standardization of the coarse rough surface of Gauss that simulation produces and non-gaussian rough surface, maximal value is 1) overlap, almost with corresponding theory autocorrelation function curve co-insides, because 512(reflection simulation length of counting is simulated in y direction, originally counted is 128) slightly reduce with the ratio of y direction auto-correlation length 50 μm, cause when auto-correlation function value is less, the coarse rough surface of Gauss of simulation generation has certain departing to the y orientation average autocorrelation function curve of non-gaussian rough surface and corresponding theoretical autocorrelator trace, this counts due to simulation to reflect simulation length, the simulation length in certain direction and the reduction of corresponding auto-correlation lenth ratio, the similar wave number order in this direction can be caused to tail off, the statistical nature in this direction is caused to depart from.As the statistical precision of autocorrelation function will be improved, can be counted by increase and reduce auto-correlation length.

Table 3 non-gaussian surface (β _x=β _y=5, μ=0, σ=1.6, S _kz=-1, K _uz=4, count 128 × 128) height statistics before Fourth-order moment simulation accuracy calculate

Table 4 non-gaussian surface (β _x=5, β _y=50, μ=0, σ=1.6, S _kz=1, K _uz=5 count 512 × 512) height statistics before Fourth-order moment simulation accuracy calculate

Claims (6)

1. a digitized simulation method for any rough surface of 3D, is characterized in that, comprises following step:

1) when the deflection angle value of any rough surface equals 0 and kurtosis value approximates 3, its surface topography meets Gaussian distribution, otherwise its surface topography meets non-gaussian distribution;

2) if step 1) in rough surface morphology be Gaussian distribution, utilize the analogy method of three-dimensional Gaussian rough surface to obtain rough surface high degree of sequence z, complete the digitized simulation of rough surface;

3) if step 1) in rough surface morphology be non-gaussian distribution, utilize the analogy method of three-dimensional non-gaussian rough surface to obtain rough surface high degree of sequence Z, complete the digitized simulation of rough surface;

Wherein, step 2) in the analogy method of three-dimensional Gaussian rough surface, specifically comprise following steps:

(1) generate discrete autocorrelation function matrix R (m, n) by given autocorrelation function, wherein m, n represent that number is a little got in x, y direction respectively;

(2) Fast Fourier Transform (FFT) is carried out to discrete autocorrelation function matrix R (m, n), obtain the power spectrum density P (m, n) of Gauss's rough surface;

$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) suppose that the power spectrum density constant C of white noise equals 1, utilize the power spectrum density obtained, calculation of transfer function H (m, n);

$H=\sqrt{P}$

(4) utilize white noise maker, obtain white noise sequence η (m, n), utilize inversefouriertransform to generate sequence of phase angles Φ;

$\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) with the sequence of phase angles Φ that step (4) obtains, calculate phase angle correlated series exp (i Φ), and inversefouriertransform is carried out to it, upgrade white noise sequence η (m, and Fourier transform A (m, n) 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) with the transfer function H (m that step (3) obtains, n), step 5) the white noise Fourier transform sequence A (m that obtains, n), ask the dot-product operation of sequence H and sequence A, and inversefouriertransform is carried out to dot-product operation result, obtain the high degree of sequence z of Gaussian surface ₀(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) with Gauss's high degree of sequence z that step (6) generates ₀, short transverse translation and convergent-divergent are carried out to it, obtain the Gauss's high degree of sequence z meeting average and standard deviation, complete the digitized simulation of Gauss's rough surface.

2. the digitized simulation method of any rough surface of a kind of 3D according to claim 1, described to z ₀short transverse carries out the method for translation and convergent-divergent, specifically comprises following steps:

1) to high degree of sequence z ₀carry out translation, make its average equal given μ;

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

Wherein mean2 (z ₀) for asking sequence z ₀the function of average, μ is height average;

2) to the high degree of sequence z after translation ₀carry out convergent-divergent in short transverse, obtain Gauss's height sequence that standard deviation is σ

Row z,

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

Wherein std2 (z ₀) for asking sequence z ₀standard deviation function.

3. the digitized simulation method of any rough surface of a kind of 3D according to claim 1, the analogy method of described three-dimensional non-gaussian rough surface, specifically comprises following steps:

1) with the Gauss's high degree of sequence z obtained in the analogy method of three-dimensional Gaussian rough surface, utilize the standardized method of gaussian sequence, obtain standardized normal distribution series Z ₀;

2) with step 1) the middle standardized normal distribution series Z obtained ₀with Fourth-order moment statistical nature before high degree of sequence as input, utilize Johnson or Pearson converting system to carry out non-gaussian conversion, obtain non-gaussian sequence Z ₁;

3) step 2 is utilized) the non-gaussian sequence Z that generates ₁, translation and convergent-divergent are carried out to it, obtain and meet average and standard deviation height statistical nature non-gaussian sequence Z, meanwhile, calculate its measure of skewness and kurtosis height statistical nature, judge whether its precision meets the demands; As met accuracy requirement, then complete the digitized simulation of non-gaussian rough surface, otherwise, continue to perform later step;

4) with step 3) the non-gaussian sequence Z that obtains ₁, utilize inversefouriertransform to generate sequence of phase angles Φ, upgrade white noise sequence η (m, n) and Fourier transform A (m, n) thereof;

5) with step 5) obtain the transfer function H (m, n) that sequence A and claim 1 obtain, by asking the dot-product operation of sequence A and H, and inversefouriertransform being carried out to dot-product operation result, obtaining the high degree of sequence z of new Gaussian surface ₀(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) with step 5) Gauss's high degree of sequence z of generating ₀, utilize the standardized method of gaussian sequence, obtain and meet the new standardized normal distribution sequence Z that average and standard deviation are respectively 0 and 1 ₀;

7) step 2 is repeated) ~ 3), until non-gaussian sequence Z meets the accuracy requirement of measure of skewness and kurtosis.

4. the digitized simulation method of any rough surface of a kind of 3D according to claim 3, the standardized method of described gaussian sequence, specifically comprises following steps:

1) to high degree of sequence z ₀carry out translation, make its average equal 0;

z ₀＝z ₀-mean2(z ₀)

2) to the high degree of sequence z after translation ₀carry out convergent-divergent in short transverse, obtain the standard gaussian sequence of standard deviation 1, Z ₀=z ₀/ std2 (z ₀).

5. the digitized simulation method of any rough surface of a kind of 3D according to claim 3, described Johnson or the Pearson converting system that utilizes carries out non-gaussian conversion, and concrete steps are as follows:

1) preferential employing Johnson non-gaussian converting system;

Lognormality fit approach S _l: Z ₁=ξ+λ exp [(Z ₀-γ)/δ] (ξ < η)

Bounded system fit approach S _b: z=ξ+λ e ^(η-γ)/δ/ [1+e ^(η-γ)/δ] (ξ < η < ξ+λ)

Boundless system fit approach S _u: z=ξ+λ sinh [(η-γ)/δ]

γ, δ are form parameter, and λ, ξ are respectively scale parameter and location parameter, and γ, δ, λ and ξ determine by average, standard deviation, measure of skewness and kurtosis equal altitudes statistical nature input parameter;

2) if bounded system S _bfit approach cannot be restrained, then adopt Pearson non-gaussian converting system, to improve the precision of non-gaussian converting system:

Utilize probability density function P to meet the following differential equation, obtain probability density function and height distribution function according to average, variance, measure of skewness and kurtosis equal altitudes statistical parameter:

$\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 in formula ₀, c ₁with c ₂determined by average, standard deviation, measure of skewness and kurtosis equal altitudes statistical nature input parameter.

6. the digitized simulation method of any rough surface of a kind of 3D according to claim 3, the described inversefouriertransform that utilizes generates sequence of phase angles Φ, upgrades white noise sequence η (m, n) and Fourier transform A (m thereof, n), specifically following steps are comprised:

1) with the non-gaussian sequence Z obtained in the analogy method of three-dimensional non-gaussian rough surface ₁, by inversefouriertransform, generate new white noise sequence of phase angles Φ (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) with step 1) the sequence of phase angles Φ (m, n) that obtains, calculate phase angle correlated series exp (i Φ), and inversefouriertransform is carried out to it, upgrade white noise sequence η (m, n) and Fourier transform A (m, n) thereof;

$\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) the Fourier transform sequence A (m, n) of the white noise sequence η (m, n) in Gauss's rough surface digitized simulation is upgraded;

$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}.$