CN102496173A - 一种用于采样点分离的数字偏移方法 - Google Patents

Abstract

本发明提供了一种计算机图像渲染方法，该方法基于(t，s)序列，使用数字偏移实现采样点的分离，与基于Halton序列，使用Cranley-Patterson旋转实现采样点分离的方法相比，本发明能够显著提高图像的渲染质量。

Description

一种用于采样点分离的数字偏移方法

技术领域

本发明涉及一种计算机图像渲染方法，特别是一种利用拟随机数(quasi-randomnumbers)确定采样点位置的渲染方法。

背景技术

在计算机图像渲染领域，为模拟各种光影效果，例如，景深、运动模糊、全局光照等，需要对具有如下形式的积分进行求值：

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}---\left(1\right)$

其中x为一s维变量，I^s为s维单位立方体，函数f(·)由所要模拟的效果及输入的场景数据决定。公式(1)一般无法精确求解，在实际的渲染过程中总是使用数值计算的方法求其近似值：

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}\≈\frac{1}{N}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(x_i\right)---\left(2\right)$

确定公式(2)中采样点位置x_i的方法可分为两类，一类是使用随机数的方法，另一类是使用拟随机数的方法。拟随机数是一种具有良好结构的分布模式，在公式(2)的积分估计中能够提供比随机数更高的效率。

采样点的分离(splitting)是另一种提高公式(2)效率的方法，该方法将公式(2)中的变量x分解为两个或更多的部分x＝(x⁽¹⁾，x⁽²⁾，...，x⁽ⁿ⁾)，其中x^(p)的维度为s_p，1≤p≤n，在对高维部分求值的过程中保持低维部分不变：

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}=\underset{I^{s_1}}{\∫}\underset{I^{s_2}}{\∫}...\underset{I^{s_n}}{\∫}f(x^{\left(1\right)},x^{\left(2\right)},...,x^{\left(n\right)}){\mathrm{dx}}^{\left(1\right)}{\mathrm{dx}}^{\left(2\right)}...{\mathrm{dx}}^{\left(n\right)}$

$\≈\frac{1}{N_1}\underset{i_1=0}{\overset{N_1-1}{\mathrm{\Σ}}}\frac{1}{N_2}\underset{i_2=0}{\overset{N_2-1}{\mathrm{\Σ}}}...\frac{1}{N_n}\underset{i_n=0}{\overset{N_n-1}{\mathrm{\Σ}}}f(x_{i_1}^{\left(1\right)},x_{i_2}^{\left(2\right)},...x_{i_n}^{\left(n\right)})---\left(3\right)$

其中N₁，N₂，...，N_n为采样点在不同层次上的分离次数。

为同时获得使用拟随机数和采样点分离的好处，美国专利US Patent No.7,358,971公开日期2007年9月27日提供了一种基于Halton序列，使用Cranley-Patterson旋转(Cranley-Patterson rotations)实现采样点分离的方法。

上述Halton序列为拟随机数序列的一种，(t，s)序列为另一种拟随机数序列且能提供优于Halton序列的质量。生成b进制(t，s)序列中第i个点的第j维分量的小数部分数字矢量的方法如下：

$d_i^{\left(j\right)}=C^{\left(j\right)}\left(\begin{array}{c}{b}_0\left(i\right)\\ b_1\left(i\right)\\ \·\\ \·\\ \·\\ b_{M-1}\left(i\right)\end{array}\right)---\left(4\right)$

其中

为(t，s)序列的第j维生成矩阵(generator matrices)，其系数

为有限域F_b中元素。b_r(i)表示整数i在b进制表示下的第r位数字，即

M为所选择的在b进制表示下的最大精度。已知小数部分数字矢量

求其对应值

的方法为：

$x_i^{\left(j\right)}=(b^{-1},b^{-2},...,b^{-M})d_i^{\left(j\right)}---\left(5\right)$

(t，s)序列的例子包括：Sobol序列、Faure序列、Niederreiter序列和Niederreiter-Xing序列等。不同的(t，s)序列的区别在于其生成矩阵的构造方法不同。

任何拟随机数都可直接用于公式(2)，但要用于公式(3)则须作适当的变换。Cranley-Patterson旋转和数字偏移(digital shifts)均为满足所需要求的变换方式，但Cranley-Patterson旋转会破坏拟随机数的结构，而数字偏移则不会。

b进制数字偏移的计算方法为：设有两组在b进制表示下的数字矢量，分别为α＝(α_M-1，α_M-2，...，α₀)^T和β＝(β_M-1，β_M-2，...，β₀)^T，用

表示数字偏移运算，将用α对β作数字偏移得到的结果记为

则γ＝(γ_M-1，γ_M-2，...，γ₀)^T，其中γ_r＝(α_r+β_r)mod b，0≤r≤M-1。除直接进行数字偏移外，有时还需先将左操作数的数字位向右移动E位后再作运算，其中E为一非负整数，将此类数字偏移运算表示为

若

则当M-E≤r≤M-1时γ_r＝β_r，当0≤r＜M-E时γ_r＝(α_r+E+β_r)mod b。

发明内容

本发明提供了一种计算机图像渲染方法，该方法基于(t，s)序列，使用数字偏移实现采样点的分离。

实现本发明技术目的的方案，包括以下步骤：

(A)对于渲染图像时需要求值的积分

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}$

使用数值计算的方法求其近似值，包括以下步骤：

将变量x分解为两个或更多的部分x＝(x⁽¹⁾，x⁽²⁾，...，x⁽ⁿ⁾)，其中x^(p)的维度为s_p，1≤p≤n，在对高维部分求值的过程中保持低维部分不变

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}=\underset{I^{s_1}}{\∫}\underset{I^{s_2}}{\∫}...\underset{I^{s_n}}{\∫}f(x^{\left(1\right)},x^{\left(2\right)},...,x^{\left(n\right)}){\mathrm{dx}}^{\left(1\right)}{\mathrm{dx}}^{\left(2\right)}...{\mathrm{dx}}^{\left(n\right)}$

$\≈\frac{1}{N_1}\underset{i_1=0}{\overset{N_1-1}{\mathrm{\Σ}}}\frac{1}{N_2}\underset{i_2=0}{\overset{N_2-1}{\mathrm{\Σ}}}...\frac{1}{N_n}\underset{i_n=0}{\overset{N_n-1}{\mathrm{\Σ}}}f(x_{i_1}^{\left(1\right)},x_{i_2}^{\left(2\right)},...x_{i_n}^{\left(n\right)});$

(B)确定步骤(A)中所述变量

的值的方法，其中0≤i_p≤N_p-1，1≤p≤n，包括以下步骤：

$x_{i_1,1}^{\left(1\right)}=\frac{{\tilde{i}}_1}{N_1};$

$x_{i_1,j}^{\left(1\right)}=h\left(d_{{\tilde{i}}_1}^{(j-1)}\right),$ 2≤j≤s₁；

${\tilde{i}}_1=i_1;$

$x_{i_p,1}^{\left(p\right)}=\frac{b^{m_p}}{N_p}h(d_{{\tilde{i}}_1}^{(s_1+...{+s}_{p-1})}\⊕d_{{\tilde{i}}_2}^{(s_2+...+s_{p-1})}\⊕...\⊕d_{{\tilde{i}}_{p-1}}^{\left(s_{p-1}\right)}{\⊕}_{E_1\left(m_p\right)}b\left({\tilde{i}}_p{b}^{M-m_p}\right)),$ 2≤p≤n；

$x_{i_p,j}^{\left(p\right)}=h(d_{{\tilde{i}}_1}^{(s_1+...+s_{p-1}+j-1)}\⊕d_{{\tilde{i}}_2}^{(s_2+...+s_{p-1}+j-1)}\⊕...\⊕d_{{\tilde{i}}_{p-1}}^{(s_{p-1}+j-1)}{\⊕}_{E_j\left(m_p\right)}{d}_{{\tilde{i}}_p}^{(j-1)}),$

2≤j≤s_p，2≤p≤n；

${\tilde{i}}_p=b^{M}h(d_{{\tilde{i}}_1}^{(s_1+...+s_{p-1})}\⊕d_{{\tilde{i}}_2}^{(s_2+...+s_{p-1})}\⊕...\⊕d_{{\tilde{i}}_{p-1}}^{\left(s_{p-1}\right)}{\⊕}_{E_1\left(m_p\right)+(M-m_p)}b\left(i_p\right)),$ 2≤p≤n；

其中，

表示

的第j维分量的值；

$d_i^{\left(j\right)}:C^{\left(j\right)}{(b_0\left(i\right),b_1\left(i\right),...,b_{M-1}\left(i\right))}^T,$ 其中

为一b进制(t，s)序列的第j维生成矩阵，b_r(i)表示整数i在b进制表示下的第r位数字，即

M为所选择的在b进制表示下的最大精度；

b(i)表示整数i在b进制表示下的各位数字所组成的数字矢量，即b(i)：＝(b_M-1(i)，b_M-2(i)，...，b₀(i))^T；

函数h(·)将b进制数字矢量转化为对应的小数值，对一b进制数字矢量d＝(d_M-1，d_M-2，...，d₀)^T，h(d)：＝(b^-1，b^-2，...，b^-M)(d_M-1，d_M-2，...，d₀)^T；

符号

表示数字偏移运算，符号表示先将左操作数的数字位向右移动E位后再作数字偏移运算；

m_p为满足

的最小非负整数；

E_j(m_p)为满足

为线性无关矢量组的非负整数；

(C)若其中s_max为所选择的最大维度，则将(x⁽¹⁾，x⁽²⁾，...，x⁽ⁿ⁾)分组为

其中p₀＝0，p_q+1＝n，

$\underset{p=p_v+1}{\overset{p_v+1}{\mathrm{\Σ}}}{s}_p\≤s_{\max },$ 0≤v≤q；

对于变量组使用步骤(B)所述方法确定其值；

对于变量组

其中1≤v≤q，确定其值的方法包括以下步骤：

$x_{i_{p_v+1},1}^{(p_v+1)}=\frac{b^{m_{p_v+1}}}{N_{p_v+1}}h\left(d_{i^{\left(v\right)}}^{\left(1\right)}{\⊕}_{E_1\left(m_{p_v+1}\right)}b\left({\tilde{i}}_{p_v+1}{b}^{M-m_{p_v+1}}\right)\right),$ 1≤v≤q；

$x_{i_{p_v+1},j}^{(p_v+1)}=h\left(d_{i^{\left(v\right)}}^{\left(j\right)}{\⊕}_{E_j\left(m_{p_v+1}\right)}{d}_{{\tilde{i}}_{p_v+1}}^{(j-1)}\right),$ 1≤v≤q， $2\≤j\≤s_{p_v+1};$

${\tilde{i}}_{p_v+1}=b^{M}h\left(d_{i^{\left(v\right)}}^{\left(1\right)}{\⊕}_{E_1\left(m_{p_v+1}\right)+(M-m_{p_v+1})}b\left(i_{p_v+1}\right)\right),$ 1≤v≤q；

$x_{i_{p_v+\mathrm{\Δ}},1}^{(p_v+\mathrm{\Δ})}=\frac{b^{m_{p_v+\mathrm{\Δ}}}}{N_{p_v+\mathrm{\Δ}}}h(d_{i^{\left(v\right)}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+1)}\⊕d_{{\tilde{i}}_{p_v+1}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1})}\⊕d_{{\tilde{i}}_{p_v+2}}^{(s_{p_v+2}+...+s_{p_v+\mathrm{\Δ}-1})}$

$\⊕...\⊕d_{{\tilde{i}}_{p_v+\mathrm{\Δ}-1}}^{\left(s_{p_v+\mathrm{\Δ}-1}\right)}{\⊕}_{E_1\left(m_{p_v+\mathrm{\Δ}}\right)}b\left({\tilde{i}}_{p_v+\mathrm{\Δ}}{b}^{M-m_{p_v+\mathrm{\Δ}}}\right)),$

1≤v≤q，2≤Δ≤p_v+1-p_v；

$x_{i_{p_v+\mathrm{\Δ}},j}^{(p_v+\mathrm{\Δ})}=h(d_{i^{\left(v\right)}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+j)}\⊕d_{{\tilde{i}}_{p_v+1}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+j-1)}\⊕d_{{\tilde{i}}_{p_v+2}}^{(s_{p_v+2}+...+s_{p_v+\mathrm{\Δ}-1}+j-1)}$

$\⊕...\⊕d_{{\tilde{i}}_{p_v+\mathrm{\Δ}-1}}^{(s_{p_v+\mathrm{\Δ}-1}+j-1)}{\⊕}_{E_j\left(m_{p_v+\mathrm{\Δ}}\right)}{d}_{{\tilde{i}}_{p_v+\mathrm{\Δ}}}^{(j-1)}),$

1≤v≤q，2≤Δ≤p_v+1-p_v， $2\≤j\≤s_{p_v+\mathrm{\Δ}};$

${\tilde{i}}_{p_v+\mathrm{\Δ}}=b^{M}h(d_{i^{\left(v\right)}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+1)}\⊕d_{{\tilde{i}}_{p_v+1}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1})}\⊕d_{{\tilde{i}}_{p_v+2}}^{(s_{p_v+2}+...+s_{p_v+\mathrm{\Δ}-1})}$

$\⊕...\⊕d_{{\tilde{i}}_{p_v+\mathrm{\Δ}-1}}^{\left(s_{p_v+\mathrm{\Δ}-1}\right)}{\⊕}_{E_1\left(m_{p_v+\mathrm{\Δ}}\right)+(M-m_{p_v\mathrm{\Δ}})}b\left(i_{p_v+\mathrm{\Δ}}\right)),$

1≤v≤q，2≤Δ≤p_v+1-p_v；

$i^{\left(v\right)}={\mathrm{\π}}_v({\tilde{i}}_1{b}^{m_2+...+m_{p_v}}+{\tilde{i}}_2{b}^{m_3+...+m_{p_v}}+...+{\tilde{i}}_{p_v-1}{b}^{m_{p_v}}+{\tilde{i}}_{p_v}),$ 1≤v≤q；

其中π_v(·)为一置换操作，对一整数

其中

$b_{M-1}^{\left(v\right)}\left(x\right)=(b_{M-1}\left(x\right)+{\mathrm{\ξ}}_v),\mathrm{mod}b$

$b_{M-2}^{\left(v\right)}\left(x\right)=(b_{M-2}\left(x\right)+{\mathrm{\ξ}}_{v,b_{M-1}\left(x\right)}),\mathrm{mod}b$

$b_{M-3}^{\left(v\right)}\left(x\right)=(b_{M-3}\left(x\right)+{\mathrm{\ξ}}_{v,b_{M-1}\left(x\right),b_{M-2}\left(x\right)}),\mathrm{mod}b$

$b_0^{\left(v\right)}\left(x\right)\stackrel{\stackrel{\stackrel{\·}{\·}}{\·}}{=}(b_0\left(x\right)+{\mathrm{\ξ}}_{v,b_{M-1}\left(x\right),b_{M-2}\left(x\right),...,b_1\left(x\right)}),\mathrm{mod}b$

其中ξ_v，ξ_v，0，ξ_v，1，...，ξ_v，b-1，ξ_v，0，0，...，ξ_{v，b-1，b-1}，...，ξ_{v，b-1，...，b-1}，1≤v≤q，为相互独立的随机数。

本发明具有以下优点：

由于本发明基于(t，s)序列，使用数字偏移实现采样点的分离，与基于Halton序列，使用Cranley-Patterson旋转实现采样点分离的方法相比，本发明能够显著提高图像的渲染质量。

具体实施方式

以下结合实施例对本发明进行详细说明。

表1列出了计算2进制Niederreiter序列的C++代码示例：

表1

在表1所示代码中，最大精度设为32位。设定一个最大精度的原因是，在数值计算中，既不可能也无必要进行无限精度的计算。

在表1所示代码中，最大维度设为12维。设定最大维度的作用在于，一方面，计算(t，s)序列需要用到的生成矩阵的数据需事先存储在计算机可读的介质中，而现实中的任何设备都不可能存储无限多的数据；另一方面，拟随机数的质量虽优于随机数，但随着维度的增大，拟随机数的质量会逐渐下降，所以在实际中也无必要使用无限维的拟随机数。

对不同的(t，s)序列，最大维度的选取要根据其质量参数(quality parameters)t和设定的最大精度来决定，其原则是，所选择的最大维度的质量参数不应超过设定的最大精度。

表2列出了2进制Sobol序列、Niederreiter序列和Niederreiter-Xing序列在不同维度中的质量参数。该参数越小表示质量越高，最小为0。

表2

维度	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	20
																	Sobol	0	0	1	3	5	8	11	15	19	23	27	31	35	40	45	71
Niederreiter	0	0	1	3	5	8	11	14	18	22	26	30	34	38	43	68
																	Niederreiter-Xing	0	0	1	1	2	3	4	5	6	8	9	10	11	13	15	21

从上表可以看出，对于2进制Sobol序列和Niederreiter序列，如果最大精度设为32位，则最大维度应选为12维，因为12维的Sobol序列和Niederreiter序列的质量参数分别为31和30，未超过32，而13维的Sobol序列和Niederreiter序列的质量参数分别为35和34，超过了32。

为充分发挥拟随机数的效率，在实现采样点分离的数字偏移运算中有时需先将左操作数的数字位向右移动若干位后再作数字偏移运算。将左操作数的数字位向右移动的位数记为E_j(m)，其中1≤j≤s，m为满足b^m≥N的最小非负整数，s为采样的维度，N为采样点的数量，b为(t，s)序列的基数，则E_j(m)须满足

1≤j≤s，0≤k＜E_j(m)}为线性无关矢量组，其中为(t，s)序列的第j维生成矩阵的系数。

对任一(t，s)序列，当m＞t时，任何满足

的非负整数E_j(m)均符合上述要求。尽管对E_j(m)的选择并不唯一，但为了使各个维度的采样质量尽可能相同，一种选择E_j(m)的方式是：若m≤t，则E_j(m)＝0，1≤j≤s；若m＞t，设m-t＝qb+r，其中0≤q，0≤r＜b，当1≤j≤r时E_j(m)＝q+1，当r＜j≤s时E_j(m)＝q。

需要注意的是，上述E_j(m)的选择方式可用于任何(t，s)序列，但对于Sobol序列和Niederreiter序列这种质量参数在低维部分小于高维部分的(t，s)序列，应根据不同维度的质量参数来选择E_j(m)以提高效率。表3列出了在对2进制Niederreiter序列使用数字偏移实现采样点分离的过程中，计算左操作数的数字位向右移动位数的C++代码示例：

表3

使用(t，s)序列需要设定一个最大维度，但渲染图像时所需采样点的维度可能超过该设定值。对于这种情况，需要对(t，s)序列的序号执行一次置换操作，置换后的序号可用于原(t，s)序列生成新的分量。表4列出了实现该置换操作的C++代码示例：

表4

表5列出了在对2进制Niederreiter序列使用数字偏移实现采样点分离的过程中，计算所需偏移量的C++代码示例，其中用到了在表1和表4中所示代码：

表5

表6列出了基于2进制Niederreiter序列，使用数字偏移实现采样点的分离的C++代码示例，其中用到了在表1、表3和表5中所示代码：

表6

表7列出了在渲染过程中使用表6所示代码的C++代码示例：

表7

在表6所示代码中，对采样点在不同层次上的分离均使用了相同的(t，s)序列，即2进制Niederreiter序列，除使用相同的(t，s)序列外也可使用基数相同但类型不同的(t，s)序列，例如，2进制Niederreiter-Xing序列。Niederreiter-Xing序列的质量虽优于Sobol序列和Niederreiter序列，但在使用上存在一个限制，即必须事先确定总的采样维度，因而不适宜作为图像渲染中一种普遍适用的方法，不过如果已知当前采样点的分离不会有更高维的采样，则可用Niederreiter-Xing序列替换Niederreiter序列以获得更好的渲染效果。

一种(t，s)序列的扩展称之为随机化(t，s)序列(randomized(t，s)sequences)，其作用在于能估计积分近似值的误差大小。虽然渲染图像时不必去估计积分近似值的误差大小，但这种随机化方法可用于在动画的渲染中缓解图像中的瑕疵，即对动画渲染中不同的图像帧使用不同的随机数对采样点进行随机化，以淡化可能存在的瑕疵图案。表8列出了在表5所示代码的基础上增加了随机化步骤的C++代码示例：

表8

Claims (4)

1.一种计算机图像渲染方法，其特征为，包括以下步骤：

(A)对于渲染图像时需要求值的积分

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}$

使用数值计算的方法求其近似值，包括以下步骤：

将变量x分解为两个或更多的部分x＝(x⁽¹⁾，x⁽²⁾，...，x⁽ⁿ⁾)，其中x^(p)的维度为s_p，1≤p≤n，在对高维部分求值的过程中保持低维部分不变

$\underset{I^s}{\∫}f\left(x\right)\mathrm{dx}=\underset{I^{s_1}}{\∫}\underset{I^{s_2}}{\∫}...\underset{I^{s_n}}{\∫}f(x^{\left(1\right)},x^{\left(2\right)},...,x^{\left(n\right)}){\mathrm{dx}}^{\left(1\right)}{\mathrm{dx}}^{\left(2\right)}...{\mathrm{dx}}^{\left(n\right)}$

$\≈\frac{1}{N_1}\underset{i_1=0}{\overset{N_1-1}{\mathrm{\Σ}}}\frac{1}{N_2}\underset{i_2=0}{\overset{N_2-1}{\mathrm{\Σ}}}...\frac{1}{N_n}\underset{i_n=0}{\overset{N_n-1}{\mathrm{\Σ}}}f(x_{i_1}^{\left(1\right)},x_{i_2}^{\left(2\right)},...x_{i_n}^{\left(n\right)});$

(B)确定步骤(A)中所述变量

的值的方法，其中0≤i_p≤N_p-1，1≤p≤n，包括以下步骤：

$x_{i_1,1}^{\left(1\right)}=\frac{{\tilde{i}}_1}{N_1};$

$x_{i_1,j}^{\left(1\right)}=h\left(d_{{\tilde{i}}_1}^{(j-1)}\right),$ 2≤j≤s₁；

${\tilde{i}}_1=i_1;$

$x_{i_p,1}^{\left(p\right)}=\frac{b^{m_p}}{N_p}h(d_{{\tilde{i}}_1}^{(s_1+...{+s}_{p-1})}\⊕d_{{\tilde{i}}_2}^{(s_2+...+s_{p-1})}\⊕...\⊕d_{{\tilde{i}}_{p-1}}^{\left(s_{p-1}\right)}{\⊕}_{E_1\left(m_p\right)}b\left({\tilde{i}}_p{b}^{M-m_p}\right)),$ 2≤p≤n；

$x_{i_p,j}^{\left(p\right)}=h(d_{{\tilde{i}}_1}^{(s_1+...{+s}_{p-1}+j-1)}\⊕d_{{\tilde{i}}_2}^{(s_2+...+s_{p-1}+j-1)}\⊕...\⊕d_{{\tilde{i}}_{p-1}}^{(s_{p-1}+j-1)}{\⊕}_{E_1\left(m_p\right)}{d}_{{\tilde{i}}_p}^{(j-1)}),$

2≤j≤s_p，2≤p≤n；

${\tilde{i}}_p=b^{M}h(d_{{\tilde{i}}_1}^{(s_1+...+s_{p-1})}\⊕d_{{\tilde{i}}_2}^{(s_2+...+s_{p-1})}\⊕...\⊕d_{{\tilde{i}}_{p-1}}^{\left(s_{p-1}\right)}{\⊕}_{E_1\left(m_p\right)+(M-m_p)}b\left(i_p\right)),$ 2≤p≤n；

其中，

表示的第j维分量的值；

$d_i^{\left(j\right)}:=C^{\left(j\right)}{(b_0\left(i\right),b_1\left(i\right),...,b_{M-1}\left(i\right))}^T,$ 其中

为一b进制(t，s)序列的第j维生成矩阵，b_r(i)表示整数i在b进制表示下的第r位数字，即

M为所选择的在b进制表示下的最大精度；

b(i)表示整数i在b进制表示下的各位数字所组成的数字矢量，即b(i)：＝(b_M-1(i)，b_M-2(i)，...，b₀(i))^T；

函数h(·)将b进制数字矢量转化为对应的小数值，对一b进制数字矢量d＝(d_M-1，d_M-2，...，d₀)^T，h(d)：＝(b^-1，b^-2，...，b^-M)(d_M-1，d_M-2，...，d₀)^T；

符号

表示数字偏移运算，符号

表示先将左操作数的数字位向右移动E位后再作数字偏移运算；

m_p为满足

的最小非负整数；

E_j(m_p)为满足

为线性无关矢量组的非负整数；

(C)若

其中s_max为所选择的最大维度，则将(x⁽¹⁾，x⁽²⁾，...，x⁽ⁿ⁾)分组为 $(x_{i_{p_0+1}}^{(p_0+1)},...,x_{i_{p_1}}^{\left(p_1\right)},x_{i_{p_1+1}}^{(p_1+1)},...,x_{i_{p_2}}^{\left(p_2\right)},...,x_{i_{p_q+1}}^{(p_q+1)},...,x_{i_{p_q+1}}^{\left(p_{q+1}\right)}),$ 其中p₀＝0，p_q+1＝n，

$\underset{p=p_v+1}{\overset{p_v+1}{\mathrm{\Σ}}}{s}_p\≤s_{\max },$ 0≤v≤q；

对于变量组

使用步骤(B)所述方法确定其值；

对于变量组

其中1≤v≤q，确定其值的方法包括以下步骤：

$x_{i_{p_v+1},1}^{(p_v+1)}=\frac{b^{m_{p_v+1}}}{N_{p_v+1}}h\left(d_{i^{\left(v\right)}}^{\left(1\right)}{\⊕}_{E_1\left(m_{p_v+1}\right)}b\left({\tilde{i}}_{p_v+1}{b}^{M-m_{p_v+1}}\right)\right),$ 1≤v≤q；

$x_{i_{p_v+1},j}^{(p_v+1)}=h\left(d_{i^{\left(v\right)}}^{\left(j\right)}{\⊕}_{E_j\left(m_{p_v+1}\right)}{d}_{{\tilde{i}}_{p_v+1}}^{(j-1)}\right),$ 1≤v≤q， $2\≤j\≤s_{p_v+1};$

${\tilde{i}}_{p_v+1}=b^{M}h\left(d_{i^{\left(v\right)}}^{\left(1\right)}{\⊕}_{E_1\left(m_{p_v+1}\right)+(M-m_{p_v+1})}b\left(i_{p_v+1}\right)\right),$ 1≤v≤q；

$x_{i_{p_v+\mathrm{\Δ}},1}^{(p_v+\mathrm{\Δ})}=\frac{b^{m_{p_v+\mathrm{\Δ}}}}{N_{p_v+\mathrm{\Δ}}}h(d_{i^{\left(v\right)}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+1)}\⊕d_{{\tilde{i}}_{p_v+1}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1})}\⊕d_{{\tilde{i}}_{p_v+2}}^{(s_{p_v+2}+...+s_{p_v+\mathrm{\Δ}-1})}$

$\⊕...\⊕d_{{\tilde{i}}_{p_v+\mathrm{\Δ}-1}}^{\left(s_{p_v+\mathrm{\Δ}-1}\right)}{\⊕}_{E_1\left(m_{p_v+\mathrm{\Δ}}\right)}b\left({\tilde{i}}_{p_v+\mathrm{\Δ}}{b}^{M-m_{p_v+\mathrm{\Δ}}}\right)),$

1≤v≤q，2≤Δ≤p_v+1-p_v；

$x_{i_{p_v+\mathrm{\Δ}},j}^{(p_v+\mathrm{\Δ})}=h(d_{i^{\left(v\right)}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+j)}\⊕d_{{\tilde{i}}_{p_v+1}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+j-1)}\⊕d_{{\tilde{i}}_{p_v+2}}^{(s_{p_v+2}+...+s_{p_v+\mathrm{\Δ}-1}+j-1)}$

$\⊕...\⊕d_{{\tilde{i}}_{p_v+\mathrm{\Δ}-1}}^{(s_{p_v+\mathrm{\Δ}-1}+j-1)}{\⊕}_{E_j\left(m_{p_v+\mathrm{\Δ}}\right)}{d}_{{\tilde{i}}_{p_v+\mathrm{\Δ}}}^{(j-1)}),$

1≤v≤q，2≤Δ≤p_v+1-p_v，

${\tilde{i}}_{p_v+\mathrm{\Δ}}=b^{M}h(d_{i^{\left(v\right)}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1}+1)}\⊕d_{{\tilde{i}}_{p_v+1}}^{(s_{p_v+1}+...+s_{p_v+\mathrm{\Δ}-1})}\⊕d_{{\tilde{i}}_{p_v+2}}^{(s_{p_v+2}+...+s_{p_v+\mathrm{\Δ}-1})}$

$\⊕...\⊕d_{{\tilde{i}}_{p_v+\mathrm{\Δ}-1}}^{\left(s_{p_v+\mathrm{\Δ}-1}\right)}{\⊕}_{E_1\left(m_{p_v+\mathrm{\Δ}}\right)+(M-m_{p_v\mathrm{\Δ}})}b\left(i_{p_v+\mathrm{\Δ}}\right)),$

1≤v≤q，2≤Δ≤p_v+1-p_v；

$i^{\left(v\right)}={\mathrm{\π}}_v({\tilde{i}}_1{b}^{m_2+...+m_{p_v}}+{\tilde{i}}_2{b}^{m_3+...+m_{p_v}}+...+{\tilde{i}}_{p_v-1}{b}^{m_{p_v}}+{\tilde{i}}_{p_v}),$ 1≤v≤q；

其中π_v(·)为一置换操作，对一整数

其中

$b_{M-1}^{\left(v\right)}\left(x\right)=(b_{M-1}\left(x\right)+{\mathrm{\ξ}}_v),\mathrm{mod}b$

$b_{M-2}^{\left(v\right)}\left(x\right)=(b_{M-2}\left(x\right)+{\mathrm{\ξ}}_{v,b_{M-1}\left(x\right)}),\mathrm{mod}b$

$b_{M-3}^{\left(v\right)}\left(x\right)=(b_{M-3}\left(x\right)+{\mathrm{\ξ}}_{v,b_{M-1}\left(x\right),b_{M-2}\left(x\right)}),\mathrm{mod}b$

$b_0^{\left(v\right)}\left(x\right)\stackrel{\stackrel{\stackrel{\·}{\·}}{\·}}{=}(b_0\left(x\right)+{\mathrm{\ξ}}_{v,b_{M-1}\left(x\right),b_{M-2}\left(x\right),...,b_1\left(x\right)}),\mathrm{mod}b$

其中ξ_v，ξ_v，0，ξ_v，1，...，ξ_v，b-1，ξ_v，0，0，...，ξ_{v，b-1，b-1}，...，ξ_{v，b-1，...，b-1}，1≤v≤q，为相互独立的随机数。

2.如权利要求1.所述方法，还包括，对采样点在不同层次上的分离使用相同的(t，s)序列或基数相同但类型不同的(t，s)序列。

3.如权利要求1.所述方法，其中，所述(t，s)序列包括Sobol序列、Faure序列、Niederreiter序列和Niederreiter-Xing序列。

4.如权利要求1.所述方法，其中，所述(t，s)序列还包括随机化的(t，s)序列。