摘 要：采用基于稀疏矩阵的大规模非负最小二乘法，对大口径、微浮雕结构光学元件加工中的驻留时间进行了分析与求解，并对该算法开展了正则化研究。仿真结果表明：与传统非负最小二乘法相比，基于稀疏矩阵的大规模非负最小二乘法精度高、效率快。采用该算法仿真加工平均振幅为１．１７７ ６倍波长的大径、微浮雕结构光学元件，误差面形均方根收敛至０．０６７倍波长。
驻留时间算法是计算机控制光学表面成形（ＣＣＯＳ）的关键技术之一。国内外常见驻留时间求解方法有脉冲迭代法、傅里叶变换法、数值迭代法和矩阵求解法等。脉冲迭代方法存在计算发散的问题［１］。傅里叶变换法利用时域卷积等同频域乘积的原理，将时域卷积运算变换到频域作乘法运算，计算效率高，但需多次调整参数来实现驻留时间解的非负［２］。数值迭代法先根据初始面形误差给驻留时间函数一个初始值，然后根据一定的公式反复迭代，直至满足收敛条件，计算效率较低，某些情况下会发生振荡而无法收敛［３］。矩阵方程法是２维卷积方程离散简化为线性矩阵方程求解的方法，通过求解矩阵方程得到驻留时间函数。光学加工中，大口径、微浮雕结构的光学元件（如透镜阵列、相位板等）将带来表面轨迹点数量的增加。特别是微浮雕结构的光学元件，面形结构起伏变化较大，所对应的去除量分布均匀性差；并且要较一般折反射元件加工去除量偏大，存在大量的小空间周期结构，由此对应的面形采样点与轨迹点密度偏大。以上两种情况，将最终导致形成的驻留时间求解问题规模大、病态性强。在保证求解结果收敛性、稳定性的前提下，如何高效率计算大规模数据处理是ＣＣＯＳ技术加工大口径、微浮雕结构的光学元件亟待解决的关键技术之一。为此，本文基于矩阵方程法的思想，针对ＣＣＯＳ应用中非回转对称去除函数形成的驻留时间矩阵方程特征，研究了一种基于稀疏矩阵的大规模驻留时间算法，并采用正则化方法降低方程的病态性，最后运用该算法进行了大口径微浮雕结构面形的仿真加工。
１ 大规模驻留时间矩阵方程模型
以工件上某一角点为参考点，建立工件坐标系（Ｘ，Ｙ）。同理，建立非回转对称去除函数坐标系（Ｘ′，Ｙ′），如图１所示。
抛光时，工件表面任意一点（ｘ，ｙ）的去除量，相当于抛光模在工件表面作进给，其中心位于工件表面各点（ξ，η）对点（ｘ，ｙ）产生去除量的叠加［４－５］。当抛光模沿着走刀路径（图２）和驻留时间函数扫描完光学元件表面所有的轨迹点后，光学元件表面某点的实际材料去除量为
ｉ＝１式中：Ｎｔ是抛光模驻留光学元件表面的轨迹点总数；Ｒ（ｘｋ－ξｉ，ｙｋ－ηｉ）为抛光模驻留在（ξｉ，ηｉ）处时，抛光模在点（ｘｋ，ｙｋ）处的材料去除量；Ｔ（ξｉ，ηｉ）为抛光模在轨迹点（ξｉ，ηｉ）处的驻留时间。引入符号表示
则式（１）可写为矩阵形式
式中：Ｎｒ表示面形采样点的个数。式（２）左边为期望的材料去除量函数，即将驻留时间反卷积求解模型变为驻留时间矩阵方程求解模型［５］。
抛光过程中抛光模驻留在某一轨迹点抛光时，光学表面的材料去除区域限定在一定范围内，很多轨迹点上的去除量为零，映射到式（２）中，Ｒ矩阵中存在大量的零值。实际加工中获得矩阵Ｒ中非零值所占带宽较窄，只有几条对角线上存在非零数值，因此Ｒ是一个稀疏矩阵。当元件的尺寸越大，抛光模的轨迹点数Ｎｔ和光学元件表面形采样点Ｎｒ越多时，矩阵Ｒ中非零值所占带宽越窄，稀疏性越明显。以往驻留时间算法主要侧重线性程组解的精度，而忽略线性方程本身的特征，导致计算效率较低，精度提高受到限制。本文采用了一种新方法———基于稀疏矩阵的大规模非负最小二乘法用于驻留时间的求解。
２ 基于稀疏矩阵的大规模非负最小二乘法
算法以非负最小二乘为基础，以优化计算残差的２范数为优化目标，其标准形式为
将其进行变换为
由于ＲＴＲ是一个半正定矩阵，式（４）是一个非线性的凸规划问题。根据Ｋａｒｕｓｈ－Ｋｕｈｎ－Ｔｕｃｋｅｒ条件可以转换为单调线性互补问题，即
若Ｒ是列满秩矩阵，ＲＴＲ是正定矩阵，则式（５）为严格单调线性互补问题，而实际获得的Ｒ满足以上两个条件。针对此类问题，算法中采用块转轴运算法则进行求解，基本思想如下：
令Ｆ和Ｇ是｛１，…，ｎ｝的子集，并且满足Ｆ∩Ｇ＝｛１，…，ｎ｝，Ｆ∩Ｇ＝；将矩阵Ｒ分解为Ｒ＝［ＲＦ，ＲＧ］，这里ＲＦ∈Ｒｍ×｜Ｆ｜，ＲＧ∈Ｒｍ×｜Ｇ｜，｜Ｆ｜，｜Ｇ｜分别对应的ＲＦ与ＲＧ的列数。假设（ＴＦ，ｙＧ）是式（５）一对互补基础解，令ＴＦ是式（６）的无约束最小二乘解，则ｙＧ由式（７）计算。
若ＴＦ≥０，ｙＧ≥０，Ｔ是式（４）的最优解；若存在任意Ｔｉ＜０，ｉ∈Ｆ或ｙｉ＜０，ｉ∈Ｇ，Ｔ是式（４）的不可行解。对于不可行解，交换Ｆ和Ｇ中负值元素，由式（６），（７）更新ｙＧ，ＴＦ；然后反复循环，直到得到最优解［７－８］。
３ 正则化方法
由于基于稀疏矩阵的大规模非负最小二乘法中，采用了Ｃｈｏｌｅｓｋｙ分解法计算式（５）中矩阵乘法，计算速度特别快的同时，要求矩阵方程的条件数小于１０８，否则求解失败。针对这种情况，文中引入了Ｔｉｋｈｏｎｏｖ正则化方法来降低方程的病态性，它的基本形式为［９］
式中：α为正则化参数，一般为正数。当Ｌ为单位矩阵时，称其为标准Ｔｉｋｈｏｎｏｖ正则化。文中根据标准Ｔｉｋ－ｈｏｎｏｖ正则化，引入正则化参数，将矩阵方程转换为
将原来的矩阵Ｒ扩展为（Ｎｒ＋Ｎｔ）×Ｎｔ矩阵，同时矩阵Ｈ扩展为（Ｎｒ＋Ｎｔ）×１。改进后的式（９）彻底变成了一个求解亏秩矩阵最小二乘问题，降低了式（２）的病态性，同时方程的解要满足使式（１０）最小［６］
由此可见，正则化参数α降低驻留时间解的能量与线性方程组的病态性，但同时相当于减小了计算残差的权值，相当于在面形计算残差和驻留时间解的能量幅度之间加入一个调节权值。较大的α可以降低驻留时间解的能量与线性方程组的病态性，但同时相当于减小了面形计算残差的权值，将引起残留面形误差值变大。正则化参数α的取值至关重要。
文中采用先验法的思想，通过对不同微浮雕结构面形进行数值计算发现：当１０－６≤α≤０．０５时，误差面形均方根（ＲＭＳ）值和峰谷（ＰＶ）值与无正则化的结果相比，偏差在１％以下，这与实际加工过程中产生的误差相比，可以忽略。当然文中给出的正则化参数的范围，并不能完全解决遇到的问题。因此，该范围对于其他面形只具有参考意义，实际加工中还要根据具体情况进行选择。
４ 驻留时间算法比较
图３为典型３２０ｍｍ×３２０ｍｍ的微浮雕面形，输入ＷＹＫＯ－Ｖｉｓｉｏｎ３２分析软件中显示，ＲＭＳ值为０．２５１λ，２０点ＰＶ值为１．９０４λ（其中λ＝６３２．８ｎｍ）。在图３所示的微浮雕面形中心处取１２０ｍｍ×１２０ｍｍ区域，保持其他抛光参数与条件不变，分别采用传统非负最小二乘法、基于稀疏矩阵的大规模非负最小二乘法分别计算驻留时间，得到的两种算法仿真加工后的残留误差分布如图４所示，其中传统的非负最小二乘法采用ＭＡＴＬＡＢ中ｌｓｑｎｏｎｎｅｇ函数来实现。
由图４可见，两种算法误差分布基本一致，ＲＭＳ及ＰＶ数值表明，文中研究的算法精度比传统非负最小二乘法高。计算过程中，二者所用时间差异很大，传统非负最小二乘法计算时间６ｍｉｎ，而基于稀疏矩阵的大规模非负最小二乘法仅需０．１ｍｉｎ。
为了进一步比较计算效率，在保持其他条件不变的情况下，改变矩阵Ｒ规模，同时采用两种算法进行计算，所用时间比较如图５所示。图中横坐标表示矩阵的规模（行数），纵坐标表示算法计算所用时间。图中曲线显示：在矩阵规模小于２ ０００时，两种算法的效率相差较小；而当矩阵规模大于２ ０００时，传统非负最小二乘法所用时间呈现近似底数大于１的指数函数增大，基于稀疏矩阵的大规模非负最小二乘法所用时间远小于传统非负最小二乘法，增长趋势很缓慢。
对于大规模矩阵计算时，传统非负最小二乘法在算法实现过程中需要大量的内存空间。基于稀疏矩阵的大规模非负最小二乘法则采用顺序压缩法［１０］，即采用了３个数组存储一个稀疏矩阵。由于实际计算中矩阵Ｒ中的非零元素占总元素数量的比例不到５％，而且随着矩阵规模的增大，所占的比例越来越小。因此，基于稀疏矩阵的大规模非负最小二乘法对矩阵Ｒ的存储方法，更适合用于大规模的计算。
５ 仿真结果与分析
在图３所示３２０ｍｍ×３２０ｍｍ微浮雕面形中，最大振幅为２．２λ，平均振幅为１．１７７ ６λ。采用本文研究的算法对该面形进行仿真加工，取α＝０．０１，矩阵规模９１ ０２２×９１ ０２２，加工结果（图６）输入ＷＹＫＯ－Ｖｉｓｉｏｎ３２分析软件中显示，得到的ＲＭＳ值为０．２２９λ，２０点ＰＶ值为１．６７５λ。从图６（ａ）中可以看出，仿真加工面形总体结构与目标面形相一致。图６（ｂ）中，面形中心轮廓线在ｙ方向吻合度较高；在ｘ方向大体趋势与目标相一致，但面形轮廓中引入了一些小空间周期的结构，能量集中在中高频段，这主要是由抛光斑与轨迹参数的耦合程度引起。将仿真加工面形（图６（ａ））与目标面形（图３（ａ））进行相减，可得如图７（ａ）所示的误差面形。误差面形的ＲＭＳ值为０．０６７λ，２０点ＰＶ值为０．８１５λ，表明加工取得了比较理想的效果，从而证明了本文驻留时间算法的有效性与合理性。
６ 结 论
本文根据微浮雕结构面形的特点，分析了去除函数非回转对称时所形成的驻留时间求解方程的特征———大规模、稀疏、病态性强。采用了基于稀疏矩阵的大规模非负最小二乘法，并对之开展了正则化研究，解决了驻留时间求解的病态性问题，大幅度提高了计算效率、降低了内存需求。采用基于稀疏矩阵的大规模非负最小二乘仿真加工３２０ｍｍ×３２０ｍｍ微浮雕结构的光学元件，仿真结果表明，误差面形ＲＭＳ收敛至０．０６７λ，从而证明了该算法的合理性。但是文中采用正则化降低矩阵方程病态性的同时，也抑制了驻留时间解向量收敛性的提高，从而降低了加工面形的精度。建立正则化与驻留时间解向量收敛性的平衡机制，研究二者之间的自适应算法，是今后驻留时间求解算法的一项重要工作。
参考文献：
［１］ Ｄｕｍａｓ Ｐ Ｒ，Ｈａｌｌｏｃｋ Ｒ Ｗ，Ｐｉｓａｒｓｋｉ Ａ．Ａｐｐｌｉｃａｔｉｏｎｓ ａｎｄ ｂｅｎｅｆｉｔｓ ｏｆ“ｐｅｒｆｅｃｔｌｙ ｂａｄ”ｏｐｔｉｃａｌ ｓｕｒｆａｃｅｓ［Ｃ］／／Ｐｒｏｃ ｏｆ ＳＰＩＥ．２００８，７１０２：１－１１．
［２］ 周林，戴一帆，解旭辉，等．离子束加工去除模型和驻留时间算法［Ｊ］．纳米技术与精密工程，２００７，５（２）：１０７－１１２．（Ｚｈｏｕ Ｌｉｎ，Ｄａｉ Ｙｉｆａｎ，ＸｉｅＸｕｈｕｉ，ｅｔ ａｌ．Ｍｏｄｅｌ ａｎｄ ｍｅｔｈｏｄ ｔｏ ｄｗｅｌｌ ｔｉｍｅ ｉｎ ｉｏｎ ｂｅａｍ ｆｉｇｕｒｉｎｇ．Ｎａｎｏｔｅｃｈｎｏｌｏｇｙ ａｎｄ Ｐｒｅｃｉｓｉｏｎ Ｅｎｇｉｎｅｅｒｉｎｇ，２００７，５（２）：１０７－１１２）［３］ 李全胜，成晔，蔡复之，等．计算机控制光学表面成形驻留时间算法研究［Ｊ］．光学技术，１９９９，３（５）：６－６０．（Ｌｉ Ｑｕａｎｓｈｅｎｇ，Ｃｈｅｎｇ Ｙｅ，ＣａｉＦｕｚｈｉ，ｅｔ ａｌ．Ｄｗｅｌｌ ｔｉｍｅ ａｌｇｏｒｉｔｈｍ ｉｎ ｃｏｍｐｕｔｅｒ ｃｏｎｔｒｏｌｌｅｄ ｏｐｔｉｃａｌ ｓｕｒｆａｃｉｎｇ．Ｏｐｔｉｃａｌ Ｔｅｃｈｎｉｑｕｅ，１９９９，３（５）：６－６０）
［４］ 孙希威．磁流变抛光机床数控系统关键技术研究［Ｄ］．哈尔滨：哈尔滨工业大学，２００６：７３－８４．（Ｓｕｎ Ｘｉｗｅｉ．Ｒｅｓｅａｒｃｈ ｏｎ ｋｅｙ ＣＮＣ ｔｅｃｈｎｉｑｕｅ ｏｆ
ｔｈｅ ｍａｇｎｅｔｏｒｈｅｏｌｏｇｉｃａｌ ｆｉｎｉｓｈｉｎｇ ｍａｃｈｉｎｅ ｔｏｏｌ．Ｈａｒｂｉｎ：Ｈａｒｂｉｎ Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｔｅｃｈｎｏｌｏｇｙ，２００６：７３－８４）
［５］ 石峰，戴一凡，彭小强，等．基于矩阵运算的光学零件磁流变加工的驻留时间算法［Ｊ］．国防科技大学学报，２００９，３１（２）：１０３－１０６．（Ｓｈｉ Ｆｅｎｇ，Ｄａｉ Ｙｉｆａｎ，Ｐｅｎｇ Ｘｉａｏｑｉａｎｇ，ｅｔ ａｌ．Ｄｗｅｌｌ ｔｉｍｅ ａｌｇｏｒｉｔｈｍ ｂａｓｅｄ ｏｎ ｖｅｃｔｏｒ ｆｏｒ ＭＲＦ ｐｒｏｃｅｓｓ ｏｆ ｏｐｔｉｃｓ．Ｊｏｕｒｎａｌ ｏｆ Ｎａｔｉｏｎａｌ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｄｅ－ｆｅｎｓｅ Ｔｅｃｈｎｏｌｏｇｙ，２００９，３１（２）：１０３－１０６）
［６］ 武建芬，卢振武，张红鑫，等．光学元件离子束加工驻留时间优化求解［Ｊ］．光学学报，２０１０，３０（１）：１９２－１９７．（Ｗｕ Ｊｉａｎｆｅｎ，Ｌｕ Ｚｈｅｎｗｕ，Ｚｈａｎｇ
Ｈｏｎｇｘｉｎ，ｅｔ ａｌ．Ｏｐｔｉｍｉｚｅｄ ｄｗｅｌｌ ｔｉｍｅ ｓｏｌｕｔｉｏｎ ｆｏｒ ｏｐｔｉｃｓ ｉｎ ｉｏｎ ｂｅａｍ ｆｉｇｕｒｉｎｇ．Ａｃｔａ Ｏｐｔｉｃａ Ｓｉｎｉｃａ，２０１０，３０（１）：１９２－１９７）
［７］ Ｐｏｒｔｕｇａｌ Ｌ Ｆ，Ｊｕｄｉｃｅ Ｊ Ｊ，Ｖｉｃｅｎｔｅ Ｌ Ｎ．Ａ ｃｏｍｐａｒｉｓｏｎ ｏｆ ｂｌｏｃｋ ｐｉｖｏｔｉｎｇ ａｎｄ ｉｎｔｅｒｉｏｒ－ｐｏｉｎｔ ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｌｉｎｅａｒ ｌｅａｓｔ ｓｑｕａｒｅｓ ｐｒｏｂｌｅｍｓ ｗｉｔｈｎｏｎｎｅｇａｔｉｖｅ ｖａｒｉａｂｌｅｓ［Ｊ］．Ｍａｔｈ Ｃｏｍｐ，１９９４，６３（２０８）：６２５－６４３．
［８］ Ｃａｎｔａｒｅｌｌａ Ｊ，Ｐｉａｔｅｋ Ｍ．Ａ ｓｏｌｖｅｒ ｆｏｒ ｌａｒｇｅ ｓｐａｒｓｅ ｌｅａｓｔ ｓｑｕａｒｅｓ ｐｒｏｂｌｅｍｓ ｗｉｔｈ ｎｏｎ－ｎｅｇａｔｉｖｅ ｖａｒｉａｂｌｅｓ［Ｊ／ＯＬ］．Ｃｏｍｐｕｔｅｒ Ｓｃｉｅｎｃｅ，２００４．ｈｔ－ｔｐ：／／ａｒｘｉｖ．ｏｒｇ／ａｂｓ／ｃｓ／０４０８０２９．
［９］ 邓伟杰，郑立功，史亚莉，等．基于线性代数和正则化方法的驻留时间算法［Ｊ］．光学·精密工程，２００７，１５（７）：１００９－１０１５．（Ｄｅｎｇ Ｗｅｉｊｉｅ，Ｚｈｅｎｇ Ｌｉｇｏｎｇ，Ｓｈｉ Ｙａｌｉ，ｅｔ ａｌ．Ｄｗｅｌｌ ｔｉｍｅ ａｌｇｏｒｉｔｈｍ ｂａｓｅｄ ｏｎ ｍａｔｒｉｘ ａｌｇｅｂｒａ ａｎｄ ｒｅｇｕｌａｒｉｚａｔｉｏｎ ｍｅｔｈｏｄ．Ｏｐｔｉｃｓ ａｎｄ Ｐｒｅｃｉｓｉｏｎ Ｅｎｇｉｎｅｅｒｉｎｇ，
２００７，１５（７）：１００９－１０１５）
［１０］ 管纪文，刘大有．数据结构［Ｍ］．北京：高等教育出版社，１９８５：３２－４３．（Ｇｕａｎ Ｊｉｗｅｎ，Ｌｉｕ Ｄａｙｏｕ．Ｄａｔａ ｓｔｒｕｃｔｕｒｅ．Ｂｅｉｊｉｎｇ：Ｈｉｇｈｅｒ ＥｄｕｃａｔｉｏｎＰｒｅｓｓ，１９８５：３２－４３）
本文作者：罗丽丽， 何建国， 王亚军， 张云飞， 黄 文， 吉 方
