CN102092929B - 用于非球面加工的离子束修形加工方法 - Google Patents

Abstract

本发明公开了一种用于非球面加工的离子束修形加工方法，包括以下步骤：先通过检测、实验获取面形误差分布数据R和去除函数P；然后建立离子束修形加工路径，并测算该光学曲面各点的入射角φ；利用计算模型和驻留时间迭代算法计算驻留时间密度函数和材料去除量M，再计算各加工点处离子束的驻留时间；再根据已建立的加工路径和驻留时间，对待加工光学元件的光学表面进行线性三轴联动的数控修形加工；再进行面形误差检测，根据加工前、后的面形误差分布数据和M对P进行校正，使用校正后的去除函数重复前述步骤进行迭代加工，直到满足面形收敛精度要求，加工结束。本发明的加工方法稳定性好，加工精度高，可控性强，对设备要求低。

Description

用于非球面加工的离子束修形加工方法

技术领域

本发明涉及光学元件加工技术领域，尤其涉及一种曲面光学元件的离子束修形加工方法。

背景技术

离子束修形是应用于光学零件确定性加工的一种新技术。离子束修形是在真空环境中，使用离子源发射的离子束轰击光学镜面，利用离子轰击产生的物理溅射现象去除光学零件表面的材料。离子束修形具有纳米量级的加工精度，是高确定性、高稳定性和非接触的加工方式。离子束修形克服了传统方法修形加工过程中的边缘效应、刀具磨损和压力负载等缺点。

如图1中的A图所示，目前的离子束修形加工方法在实现对非球面等光学曲面1的加工时，为了保持去除函数稳定不变，需要使离子束2垂直入射光学曲面1，即离子束2方向保持与光学曲面1的法线方向处重合。由于非球面等光学曲面1上法线方向是空间改变的，因此加工时，离子束2的入射方向也必须是空间改变的，以处处重合光学曲面1的法线方向(如图1所示)。所以，离子束2修形加工非球面等光学曲面1时，离子源3相对工件需要5个运动自由度，即调整位置的X、Y和Z三个线性运动自由度和调整姿态的θ_x和θ_y两个旋转自由度。因此，为实现非球面等光学曲面的修形加工，目前的离子束修形加工中必须使用五轴联动，加工设备必须具备五轴联动的运动系统。而五轴联动的运动系统不仅操作复杂，制备成本高，而且加工效率和加工效果还有待提高。

发明内容

本发明所要解决的技术问题是克服现有技术的不足，提供一种稳定性好、加工精度高、可控性强、对设备要求低的用于非球面加工的离子束修形加工方法。

为解决上述技术问题，本发明提出的技术方案为一种用于非球面加工的离子束修形加工方法，包括以下步骤：

(1)检测面形误差：利用面形检测装置对待加工光学元件的非球面光学曲面进行面形误差检测，获取该光学曲面的面形误差分布数据R，记为r(x，y)；

(2)确定去除函数：应用离子束修形加工过程进行去除函数试验获取去除函数P(或者直接使用已经获取并保留下来的去除函数)，记为p(x，y)；

(3)建立加工路径：采用现有路径规划方法(例如直线式光栅路径)，根据待加工光学元件的形状和尺寸建立离子束修形加工路径；

(4)测算入射角：在垂直于待加工光学元件光轴的平面内建立X-Y直角坐标系，根据其光学曲面的曲面方程z＝h(x，y)，测算该光学曲面各点法线与垂直入射离子束的夹角，即入射角φ；

(5)确定驻留时间密度函数：利用计算模型r(x，y)＝η(x，y)·[p(x，y)*τ(x，y)]和驻留时间迭代算法，计算本次加工控制所用的驻留时间密度函数τ(x，y)和材料去除量M；

本步骤中采用了与现有技术不同的计算模型，现有技术中，由于去除函数p(x，y)保持不变，因此在求解驻留时间密度函数τ(x，y)时，仅仅通过卷积方程r(x，y)＝p(x，y)*τ(x，y)即可进行求解，即在现有的五轴联动系统加工中的材料去除量r(x，y)等于去除函数p(x，y)与驻留时间密度函数τ(x，y)的空间二维卷积；而本发明方法的后续数控修形加工步骤中，由于是采用线性三轴联动系统进行加工，去除函数p(x，y)是空间变化的，因此现有的卷积方程难以继续用于求解，利用现有技术无法实现线性三轴加工方式的驻留时间解算，因此我们采用了计算模型r(x，y)＝η(x，y)·[p(x，y)*τ(x，y)]；该计算模型表明线性三轴加工中的材料去除量r(x，y)等于通常的去除函数p(x，y)与驻留时间密度函数τ(x，y)的空间二维卷积，再乘以相应的材料去除率系数η(x，y)；因此，通过该计算模型可以求解出线性三轴加工中用于控制加工过程的驻留时间密度；

在所述计算模型中，r(x，y)即为步骤(1)中检测获取的面形误差分布数据，p(x，y)即为步骤(2)中获取的去除函数，η(x，y)为材料去除率系数(材料去除率系数是由于离子束倾斜入射导致的材料去除率变化引入的系数)，材料去除率系数η(x，y)的物理意义是指加工过程中在曲面点(x，y)处的材料去除速率相对于平面的材料去除速率之比；

根据Sigmund建立的溅射理论(Peter Sigmund，Theory of Sputtering.I.Sputtering Yield ofAmorphous and Polycrystalline Targets.Phy.Rey.，184(2)：p.383-416.1969)，溅射过程中溅射的强弱可用溅射产额Y₀来描述，溅射产额Y₀定义为平均一个入射离子溅射出的材料表面原子数；溅射产额Y₀同样是入射角φ的函数，且随入射角φ的变化而变化；此处我们引入相对溅射产额Y，相对溅射产额Y定义为Y(φ)＝Y₀(φ)/Y₀(0)，其中Y(φ)为相对溅射产额函数，Y₀(φ)为入射角φ处平均一个入射离子溅射出的材料表面原子数，Y₀(0)为入射角0°处平均一个入射离子溅射出的材料表面原子数；根据该定义，Y(0)＝1，即入射角为零时(垂直入射)，相对溅射产额为1；本步骤中，所述的材料去除率系数η(x，y)仅与入射角φ有关，即η(x，y)＝f[θ(x，y)]，其值等于相对溅射产额大小，即η(x，y)＝Y(φ)；根据Sigmund理论，相对溅射产额函数Y(φ)随着入射角φ的增大而增大，最小值Y(0)＝1，据此材料去除率系数η(x，y)总是大于或等于1，即η(x，y)≥1；可见，线性三轴加工中的材料去除量较大，采用线性三轴加工方式的材料去除效率更高，计算所得的驻留时间也将较小，因此采用线性三轴加工方式可以明显缩短离子束修形加工时间，提高修形加工的效率；

本步骤中材料去除率系数η(x，y)的计算步骤为：①计算光学曲面上加工点(x，y)处的入射角φ(x，y)；②计算相对溅射产额随入射角度变化的相对溅射产额函数Y(φ)；③根据计算所得的入射角φ(x，y)和相对溅射产额函数Y(φ)确定出光学曲面上加工点(x，y)处的材料去除率系数η(x，y)；其中，步骤②中相对溅射产额函数Y(φ)的计算可以通过实验的方法确定，也可通过现有的理论方法计算；

(6)确定加工路径上各加工点的驻留时间：根据步骤(5)中获取的驻留时间密度函数τ(x，y)，利用现有方法计算加工路径上各加工点处离子束的驻留时间；

(7)数控修形加工：根据步骤(3)中已建立的加工路径和步骤(6)计算所得的驻留时间，对待加工光学元件的光学表面进行数控修形加工；数控修形加工过程中，离子束沿平行于待加工光学元件光轴方向入射光学曲面，通过离子束发生装置的线性三轴联动系统对所述光学表面进行数控修形加工；

(8)检测加工后面形误差：利用面形检测装置对数控修形加工后的光学曲面进行面形误差检测，获取加工后的面形误差分布数据E，记为r′(x，y)；

(9)校正去除函数：根据加工前、后的面形误差分布数据和材料去除量M对所述去除函数P进行校正，校正公式为

P′＝γP

γ＝min||R-E-γM||

其中，P′为校正后的去除函数，γ为校正系数；

(10)迭代加工：使用校正后的去除函数P′替代上述步骤(2)中确定的去除函数P，并重复上述步骤(3)～(9)对所述光学元件进行迭代加工，直到所述加工后的面形误差分布数据E满足面形收敛精度要求，加工结束。

上述本发明的用于非球面加工的离子束修形加工方法中，所述步骤(5)中的驻留时间迭代算法为改进后的驻留时间迭代算法，具体的操作步骤为：

(a)输入求解条件：向所述的计算模型r(x，y)＝η(x，y)·[p(x，y)*τ(x，y)]中先输入求解条件，包括步骤(1)中检测获取的面形误差分布数据r(x，y)(即加工中期望的材料去除量)、步骤(2)中获取的去除函数p(x，y)以及材料去除率系数η(x，y)；

(b)设置迭代法计算的初始条件：

τ₀(x，y)＝(1/B)·r(x，y)/η(x，y)，其中B为加工过程的体积去除效率(即去除函数强度)，B＝∫∫p(x，y)dxdy；

仿真残留面形误差初始值e₀(x，y)＝r(x，y)-η(x，y)·[p(x，y)*τ0(x，y)]，

迭代次数k的初始值为0；

(c)计算第k步的驻留时间密度校正量

${\mathrm{\Δ}}_{\mathrm{\τ}}^k(x,y)=(1/B)\·e_k(x,y)/\mathrm{\η}(x,y),$

e_k(x，y)为第k步后的仿真残留面形误差；

(d)计算第k+1步的驻留时间密度τ_k+1(x，y)：

ξ为迭代松弛因子；

(e)驻留时间非负性检查：检查驻留时间密度τ_k+1(x，y)，如有小于零的部分，则置零；

(f)计算第k+1步的材料去除量m_k+1(x，y)：m_k+1(x，y)＝η(x，y)·[p(x，y)*τ_k+1(x，y)]；

(g)计算第k+1步的残留面形误差e_k+1(x，y)：e_k+1(x，y)＝r(x，y)-m_k+1(x，y)；

(h)判断第k+1步的驻留时间密度τ_k+1(x，y)和残留面形误差e_k+1(x，y)是否满足要求(面形收敛精度要求)，如满足要求，则输出计算结果，停止程序；如果不满足要求，则置k＝k+1，并跳到(c)步继续计算。

上述本发明方法中提到的去除函数p(x，y)、面形误差分布数据r(x，y)、驻留时间密度函数τ(x，y)、入射角φ(x，y)和材料去除率系数η(x，y)等函数，既可以为连续函数，也可为矩阵，如果是矩阵，则相应的函数计算公式则采用矩阵计算公式即可。

本发明的上述技术方案是基于以下技术原理和思路：离子束修形加工是利用离子束轰击光学表面时发生的物理溅射现象去除光学表面的材料，由于溅射现象不仅在离子束垂直入射时发生，在倾斜入射时也同样发生，而且发生程度更强烈，即倾斜入射同样能去除表面材料，而且去除速率更大；据此，离子束修形加工过程中不必使离子束垂直入射光学曲面，即使倾斜入射也同样可以达到加工效果。因此，上述本发明的加工方法提出了一种仅使用线性三轴联动系统即可对非球面等光学曲面进行修行加工的离子束修形加工方法，使用本发明的方法进行加工时，离子束方向不必处处与光学曲面的法线方向重合，可以保持离子束方向在空间上不变(即保持离子束平行于待加工光学元件的光轴即可，参见图1中的B图)。由于保持离子束方向的空间不变，因此离子束入射到光学曲面上时在各加工点处的入射角度是变化的，相应的，去除函数也是空间变化的，因此，本发明的加工工艺需要考虑入射角变化引起的去除速率变化，即在驻留时间解算时需要对去除函数的变化进行相应的计算和调整(参见上述本发明的技术方案)。

与现有技术相比，本发明的优点在于：

1.本发明采用线性三轴联动系统即能加工以往需要五轴联动系统才能加工的非球面等光学曲面，这便大大降低了非球面等光学曲面的离子束加工设备的成本；

2.本发明仅采用三轴联动系统进行加工，相比五轴联动系统，具有较高的动态特性和定位精度，可潜在地提高最终的加工精度；

3.由于倾斜入射可提高材料去除速率，因此本发明的方法加工还可明显缩短加工时间，提高加工效率；

4.由于本发明不是直接对面形误差函数进行补偿，而是对驻留时间的解算方法作了进一步的改进和优化，因此可以通过计算材料去除量来对去除函数进行修正，并用修正后的去除函数进行迭代加工，这样能够进一步提高离子束修形加工的加工精度和加工效率。

综上，本发明的加工方法具有广阔的应用前景。

附图说明

图1为离子束修形加工方式的对比示意图；其中，A图为传统五轴联动加工方式的加工示意图，B图为本发明的三轴联动加工方式的加工示意图；图中图例的含义为：

1、光学曲面；2、离子束；3、离子源。

图2为本发明实施例中检测获取的待加工非球面镜的面形误差分布图。

图3为本发明实施例中离子束修形工艺过程的去除函数分布图。

图4为本发明实施例中计算所得的非球面镜上的入射角度分布图。

图5为本发明实施例中计算所得的离子束修形工艺过程的相对溅射产额曲线图。

图6为本发明实施例中计算所得的非球面镜上的材料去除率系数分布图。

图7为本发明实施例中计算所得的驻留时间密度分布图。

图8为本发明实施例中计算所得的材料去除量分布图。

图9为本发明实施例中的非球面镜经过一次实际修形加工之后测得的面形误差分布图。

图10为本发明实施例中的非球面镜经过四次实际修形加工之后测得的面形误差分布图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明。

实施例

本实施例的离子束抛光工艺过程是在一台普通的离子束抛光设备上进行，其工艺参数为：工作气体为氩气，工作真空0.8×10^-2Pa，离子能量800eV，束电流25mA。待加工光学元件为一块尺寸为250mm×90mm的离轴非球面镜，材料为普通微晶玻璃；其镜面的离轴量为16mm，顶点曲率半径为285mm，偏心率为0.93。该镜面经过了传统的光学研抛，但是面形误差仍然较大，不能满足要求。以下采用本发明的离子束修形加工方法对该非球面镜进行离子束修形加工，具体的操作步骤如下：

1.检测面形误差：利用激光干涉仪检测上述待加工光学元件镜面的面形误差，检测获取的面形误差分布数据r(x，y)用矩阵R表示，其矩阵分布如图2所示，由图2可见，面形误差的PV值是3.3λ(λ＝632.8nm)，RMS值是0.49λ；

2.确定去除函数：应用上述的离子束抛光设备及其修形工艺过程进行去除函数试验，获取的去除函数p(x，y)用矩阵P表示，其矩阵分布如图3所示，由图3可见，该去除函数的体积去除速效率B＝31μm·mm²·min^-1；

3.建立加工路径：采用直线式光栅路径进行扫描加工，采用的路径行距和间隔均为1mm；

4.测算入射角：在垂直于待加工光学元件光轴的平面内建立X-Y直角坐标系，根据待加工光学元件镜面的曲面方程测算该镜面各点法线与垂直入射离子束的夹角，即入射角φ；测算结果用与上述面形误差矩阵R相同大小的矩阵Ф表示，矩阵分布如图4所示；

5.确定驻留时间密度函数：利用计算模型r(x，y)＝η(x，y)·[p(x，y)*τ(x，y)]和驻留时间迭代算法，计算本次加工控制所用的驻留时间密度函数τ(x，y)和材料去除量M；

5.1计算材料去除率系数η(x，y)，具体步骤为：先根据上述获取的入射角φ的矩阵Φ计算相对溅射产额随入射角度变化的相对溅射产额函数Y(φ)；根据Sigmund理论(Drueding T W，Fawcett S C，Wilson S R等.Ion beam figuring of small optical components.Optical Engineering，34(12)：p.3565-3571.1995)，相对溅射产额函数Y(φ)的计算公式为：

$Y\left(\mathrm{\φ}\right)=\frac{Y_0\left(\mathrm{\φ}\right)}{Y_0\left(0\right)}=\frac{a}{\mathrm{\μ}}{B}_1^{-1/2}\left(\mathrm{\φ}\right)\exp \left(\frac{A^2\left(\mathrm{\φ}\right)}{2B_1\left(\mathrm{\φ}\right)}\right)$

其中，

a为离子的入射深度，σ和μ分别为入射离子在平行和垂直于离子入射方向上的能量分散参数；为了先确定参数a、σ和μ的值，我们利用仿真软件SRIM2003(www.srim.org)进行仿真计算，我们仿真了5000粒能量为800eV的Ar⁺离子入射SiO₂材料的能量分散情况，仿真结果为a＝3.7nm，σ＝1.8nm，μ＝1.3nm；将仿真后的a、σ和μ值代入上述溅射产额函数Y(φ)的计算公式，即可计算出相对溅射产额函数Y(φ)，计算结果如图5所示；根据测算所得的入射角矩阵Φ(参见图4)和相对溅射产额函数Y(φ)(参见图5)，利用公式η(x，y)＝Y(φ)即可计算出该镜面上各加工点处的材料去除率系数η(x，y)，计算结果同样以矩阵H表示，其矩阵分布如图6所示；

5.2确定驻留时间密度函数：根据获取的去除函数矩阵P(参见图3)和检测获取的待加工光学元件镜面的面形误差矩阵R(参见图2)，利用本发明改进的驻留时间迭代计算方法计算加工控制所用的驻留时间密度分布，驻留时间密度函数τ(x，y)可用矩阵T表示；具体操作步骤为：

5.2.1输入求解条件：输入检测获得的面形误差矩阵R、去除函数矩阵P和上面计算所得的材料去除率系数H作为求解条件；

5.2.2计算并设置迭代法计算的初始条件：

驻留时间密度函数初始值T₀＝(1/B)·R/H，

仿真残留面形误差矩阵初始值E₀＝R-H·(P*T₀)，

迭代次数k的初始值为0；

5.2.3计算第k步的驻留时间密度校正量

E_k为第k步后的仿真残留面形误差矩阵；

5.2.4计算第k+1步的驻留时间密度T_k+1： $T_{k+1}=T_k+\mathrm{\ξ}\·{\mathrm{\Δ}}_{\mathrm{\τ}}^k(x,y);$

5.2.5驻留时间非负性检查：检查驻留时间密度T_k+1，如有小于零的部分，则置零；

5.2.6计算第k+1步的材料去除量M_k+1：M_k+1＝H·(P*T_k+1)；

5.2.7计算第k+1步的残留面形误差E_k+1：E_k+1＝R-H·(P*T_k+1)；

5.2.8判断第k+1步的驻留时间密度T_k+1和残留面形误差E_k+1是否满足要求，如满足要求，则输出计算结果：驻留时间密度矩阵T＝T_k+1；仿真残留面形误差矩阵E′＝E_k+1，材料去除量M＝M_k+1，并停止程序；如果不满足要求，则置k＝k+1，并跳到第5.2.3步继续计算；

利用上述步骤最后计算所得的驻留时间密度函数可用矩阵T表示，矩阵分布如图7所示，驻留时间密度的单位为min/mm²；计算所得的材料去除量可用矩阵M表示，矩阵分布如图8所示；

6.计算加工路径上各点的驻留时间：根据解算得到的驻留时间密度矩阵T(参见图7)，可以计算出加工路径上各加工点处离子束的驻留时间，计算结果以矩阵Г表示，驻留时间的单位为min.，由于加工路径上加工点的行距和间隔都是1mm，所以，驻留时间矩阵Γ和驻留时间密度矩阵T(图7)的数值是一样的，仅是物理意义和单位不同；

7.数控修形加工：根据已建立的加工路径和计算所得的驻留时间矩阵Γ，利用现有技术方法生成数控代码，再以数控方式对待加工光学元件表面进行数控修形加工；数控修形加工过程中，离子束沿平行于待加工光学元件光轴方向入射光学曲面，通过离子束发生装置的线性三轴联动系统对所述光学表面进行数控修形加工；

8.检测加工后面形误差：利用激光干涉仪检测加工后光学镜面的面形误差，检测所得的面形误差用矩阵E表示，检测结果如图9所示，面形误差的RMS值减小到0.41λ；

9.校正去除函数：根据实际加工结果校正去除函数P：根据公式γ＝min||R-E-γM||，可以计算出系数γ＝1.05；利用公式P′＝γP对去除函数P进行校正；

10.迭代加工：使用校正后的去除函数P’代替原来计算中的去除函数P，重复上述步骤(3)～(9)对光学元件再进行迭代加工三次。

经过以上四次实际修形加工后，该离轴非球面镜的面形误差逐步降低到了0.075λ，检测结果如图10所示，该面形误差已经满足指标要求，即可停止对该工件的加工。经过上述方法加工后，面形误差已大大减小，达到加工精度要求。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不限于上述实施例。

Claims (2)

1.一种用于非球面加工的离子束修形加工方法，包括以下步骤：

(1)检测面形误差：利用面形检测装置对待加工光学元件的非球面光学曲面进行面形误差检测，获取该光学曲面的面形误差分布数据R，记为r(x，y)；

(2)确定去除函数：应用离子束修形加工过程进行去除函数试验获取去除函数P，记为p(x，y)；

(3)建立加工路径：采用现有路径规划方法，根据待加工光学元件的形状和尺寸建立离子束修形加工路径；

(4)测算入射角：在垂直于待加工光学元件光轴的平面内建立X-Y直角坐标系，根据其光学曲面的曲面方程z＝h(x，y)，测算该光学曲面各点法线与垂直入射离子束的夹角，即入射角φ；

(5)确定驻留时间密度函数：利用计算模型r(x，y)＝η(x，y)·[p(x，y)*τ(x，y)]和驻留时间迭代算法，计算本次加工控制所用的驻留时间密度函数τ(x，y)和材料去除量M；

在所述计算模型中，r(x，y)即为步骤(1)中检测获取的面形误差分布数据，p(x，y)即为步骤(2)中获取的去除函数，η(x，y)为材料去除率系数，η(x，y)＝Y(φ)，Y(φ)＝Y₀(φ)/Y₀(0)，其中Y(φ)为相对溅射产额函数，Y₀(φ)为入射角φ处平均一个入射离子溅射出的材料表面原子数，Y₀(0)为入射角0°处平均一个入射离子溅射出的材料表面原子数；

(6)确定加工路径上各加工点的驻留时间：根据步骤(5)中获取的驻留时间密度函数τ(x，y)，利用现有方法计算加工路径上各加工点处离子束的驻留时间；

(7)数控修形加工：根据步骤(3)中已建立的加工路径和步骤(6)计算所得的驻留时间，对待加工光学元件的光学表面进行数控修形加工；数控修形加工过程中，离子束沿平行于待加工光学元件光轴方向入射光学曲面，通过离子束发生装置的线性三轴联动系统对所述光学表面进行数控修形加工；

(8)检测加工后面形误差：利用面形检测装置对数控修形加工后的光学曲面进行面形误差检测，获取加工后的面形误差分布数据E，记为r′(x，y)；

(9)校正去除函数：根据加工前、后的面形误差分布数据和材料去除量M对所述去除函数P进行校正，校正公式为

P′＝γP

γ＝min||R-E-γM||

其中，P′为校正后的去除函数，γ为校正系数；

(10)迭代加工：使用校正后的去除函数P′替代上述步骤(2)中确定的去除函数P，并重复上述步骤(3)～(9)对所述光字元件进行迭代加工，直到所述加工后的面形误差分布数据E满足面形收敛精度要求，加工结束。

2.根据权利要求1所述的用于非球面加工的离子束修形加工方法，其特征在于，所述步骤(5)中的驻留时间迭代算法的具体操作步骤为：

(a)输入求解条件：向所述的计算模型r(x，y)＝η(x，y)·[p(x，y)*τ(x，y)]中先输入求解条件，求解条件包括所述步骤(1)中检测获取的面形误差分布数据r(x，y)、所述步骤(2)中获取的去除函数p(x，y)以及材料去除率系数η(x，y)；

(b)设置迭代法计算的初始条件为：τ₀(x，y)＝(1/B)·r(x，y)/η(x，y)，其中B为加工过程的体积去除效率，且B＝∫∫p(x，y)dxdy；

仿真残留面形误差初始值e₀(x，y)＝r(x，y)-η(x，y)·[p(x，y)*τ0(x，y)]，迭代次数k的初始值为0；

(c)计算第k步的驻留时间密度校正量

e_k(x，y)为第k步后的仿真残留面形误差；

(d)计算第k+1步的驻留时间密度τ_k+1(x，y)：

ξ为迭代松弛因子；

(e)驻留时间非负性检查：检查驻留时间密度τ_k+1(x，y)，如有小于零的部分，则置零；

(f)计算第k+1步的材料去除量m_k+1(x，y)：m_k+1(x，y)＝η(x，y)·[p(x，y)*τ_k+1(x，y)]；

(g)计算第k+1步的残留面形误差e_k+1(x，y)：e_k+1(x，y)＝r(x，y)-m_k+1(x，y)；

(h)判断第k+1步的驻留时间密度τ_k+1(x，y)和残留面形误差e_k+1(x，y)是否满足要求，如满足要求，则输出计算结果，停止程序；如果不满足要求，则置k＝k+1，并跳到所述(c)步继续计算。

Images