论文部分内容阅读
摘要：对岩石三维模型重构的模拟退火算法进行并行性分析，针对其特殊性提出并行算法。在本文中，介绍了在岩石三维模型重构中使用模拟退火算法的情况，其存在的性能问题。探讨其改为并行算法的可能性，找到其与经典模拟退火应用并行性的不同的特殊之处，提出其可并行的条件及改进后的并行算法，并验证。
关键词：模拟退火、并行、三维重构
【分类号】：TM76；TM744
引言
岩石结构中含有大量不同尺寸、随机分布的孔隙，这些孔隙直接影响着岩石的物理力学性质。如果能够定量描述岩石的微细观结构并找到这些微细观结构与岩石宏观力学行为之间的关系，将为科学、准确地描述工程岩体和工程结构的响应，以及预测、控制水利水电和土木工程中的灾害控制提供有力的依据[1]。岩石重构模型是在CT扫描实验的结果上，提取岩石孔隙结构特征，构建了表征岩石孔隙分布与相关性质的两点概率函数和线性路径函数，并以此为控制函数，运用模拟退火算法重构了岩石的三维孔隙结构。
模拟退火算法SA（simulated annealing algorithm）是适合解大规模组合优化问题，特别是NP完全问题的通用有效近似算法[4]。它源于对固体退火过程的模拟，固体退火是先将固体加热至熔化，再徐徐冷却。加温时固体内部粒子能量升高，变为无序状态，冷却时，粒子内能减少，趋于有序；系统在每一个温度下都会有一个趋于平衡态的过程，最后凝固成规整晶体的热力学过程[4]。
将模拟退火算法应用到重构模型上，是以真实的孔隙分布数据作为输入资料，建立了与真实孔隙结构具有相同特征的数字孔隙模型，所建模型能较为真实反映的岩石类材料的内部结构。但该方法的应用受到收敛速度慢的困扰，尤其是当模型尺寸增大，所需的计算时间大幅增加。因此需要找到办法提高模拟退火算法在岩石重构计算中的效率，并行计算正是我们期待的方法。
问题描述
岩石重构的过程是先用CT扫描实验的结果，计算岩石的孔隙率，并构建表征岩石孔隙分布与相关性质的两点概率函数和线性路径函数。根据孔隙率构建一个初始是随机分布的三维立方体模型，在这个立方体中每个个坐标位置的值是0或者1，用来表示空隙或岩石实体。这个初始的模型就是初始状态S，而两点概率函数和线性路径函数就是用于计算能量的控制函数F（S）[2][3]。模拟退火算法中的每一次实验，就是尝试交换这个三维立方体中的一对点（一个0点和一个1点），然后计算能量变化，并判断是否接受为新的状态。
本问题可归结为优化问题，即寻找最佳状态S，使下式最小
其中：
P（r）为参考模型的两点概率函数，P’（r）为当前状态的两点概率函数
L（r）为参考模型的线性路径函数，L’（r）为当前状态的线性路径函数
应用SA算法解决该问题，经典的模拟退火算法是串行的，其算法描述如下：
1 初始化：初始的足够高的温度T，初始状态S，每个温度下的实验次数L（即Markov链长度）。这相当于将固体加热至足够高。
2 循环执行L次实验，每个实验执行以下四个任务。这相当于在一个温度下系统趋于热平衡的过程。
2.1 对S进行“扰动”，计算能量，产生新解S’
2.2 计算S’和S的能量差ΔE
2.3 判断是否接受：如果ΔE<0，接受S’为新解，否则以概率exp（-ΔE/T）接受S’
2.4 接受或拒绝新解
3 根据温度或其他条件判断是否结束，结束则输出当前S，否则降温并继续执行第二步。
并行计算
在将模拟退火算法并行的实践中主要有以下两种策略：
1. 并行操作[4]：把每个任务中的多个操作做并行处理。这个策略受到具体问题的限制，可并行完成的操作数每个问题不同。处理机个数增加到一定数量后，算法性能不再提升。
2. 并行实验[4]：在同一个温度下，n个处理机同时进行实验，产生n个新解S’1， S’2 … S’n ，从中找到能量降低的m个解S’’1， S’’2 … S’’m，然后按某个准则接受其中之一作为新解。这实际上已经与原有模拟退火算法不同：模拟退火算法中每次产生的新解S’是对前一个解S的扰动产生的，然后判断是否接受；而并行实验中的所有S’都是对同一个S扰动产生的，除了要用Metropolis接受准则判断外，还要有一个准则用来选择新解。这个策略可以充分发挥并行计算的优势。像是货郎担问题（Travelling Salesman Problem， TSP）、0-1背包问题（Zero-One Knapsack Problem， ZKP）、调度问题（Scheduling Problem， SCP）等所谓的NP完全问题都可以采用这个策略。
在岩石重构的过程中，生成新的状态并计算其能量的任务难以拆分为多个可以并行的操作，并且并行操作也难以发挥更多处理机的性能。所以我们希望用并行实验的方法来提升模拟退火算法计算岩石重构的效率。
在将并行实验策略应用到岩石重构的退火算法上时，每次并行实验在当前解S的基础上产生n个新解S’1， S’2 … S’n，从中找到能量降低的m个解S’’1， S’’2 … S’’m，根据并行实验的策略，这时应该从m个解中选择一个作为新的解。但不论是哪一个解，都只是在S的基础上交换了两个点，从整体上看这样并行与串行计算的效率几乎相同，甚至由于需要考虑并行线程之间的同步，效率反而更低。
如果能将能量降低的m个解所交换的m对点直接应用在S上而形成新的解，则在原有串行算法中交换一对点的时间内，并行算法理论上可以交换m对点，这显然是我们期待的结果。问题的关键是，虽然已知在S上交換m对点中的任何一对都会使能量下降，但是同时交换全部的m对点能量是否会下降，即m对点之间是否会有干扰。如果没有干扰，则整体能量一定会下降。 所以需要分析能量差ΔE与交换的点之间的关系：
设S为当前解，A点为S中的一个0点，B点为S中的一个1点，A的坐标为（Ax， Ay， Az），B的坐标为（Bx， By， Bz），S’为交换A、B两点后的解 ，ΔE为S和S’的能量差，即
ΔE = F（S’） – F（S）
通过对ΔE的方程进行分析，发现如下情况，ΔE只与A、B两点所在的x方向的两条线（A点一条线、B点一条线）、y方向的两条线，z方向的两条线，一共是六条线上的所有的点的值相关。也就是说，前述的m对点，如果他们的坐标都不相同，则相互之间对于能量没有干扰。
据此，我们提出岩石重构模型的并行模拟退火算法如下：
1 初始化：初始的足够高的温度T，初始状态S，每个温度下的并行实验次数L，每次并行实验的次数N（线程数）。
2 循环执行L次并行实验，每次并行实验同时有N个单独实验。设单独实验的编号为x∈（1，N）。
2.1 选取N对互不干扰的点，分别交给每个单独的实验（主线程做）
2.2 第x个实验做以下工作（每线程单独做）：
2.2.1 在S上交换一对点，并计算能量，产生临时解S’x
2.2.2 计算S’和S的能量差ΔEx。
2.2.3 判断是否接受为临时解：如果ΔEx<0，接受S’x为临时解，否则以概率exp（-ΔEx/T）接受S’x为临时解。
2.3 汇总临时解（主线程做）：设有m个被接受的临时解，即有m对可交换点，交换S中的这m对点，得到新解S’
3. 根据温度或其他条件判断是否结束，结束则输出当前S，否则降温并继续执行第二步。
仿真结果
实验使用的硬件条件如下：AMD Opteron Processor 6378 2.4G 4处理器 512G内存
岩石重构实验中我们采用如下参数：立方体为50x50x50像素，共计125000个点，降温50次，每温度尝试交换100000对点。
岩石重构模拟退火算法中，不论是串行还是并行，其初始状态都是类似的，下面为初始状态z轴方向的截图：
Fig 1 串行计算初始状态
Fig 2 并行6线程计算初始状态
重构计算完成后的结果的z轴方向截图：
Fig 3 串行计算最终结果
Fig 4 并行6线程计算最终结果
从结果图中可以看出，串行和并行计算后的结果满足统计的数据，是可以接受的结果。
下面是串行和并行多线程程序的运行时间：
Fig 5计算时间与线程个数关系图
从上面的结果看，使用并行算法后，随着并行线程数增加，运行时间在减少，在6线程时达到极值，之后随着线程增加运行时间反而增多，这主要是由于我们实验的岩石模型选的太小，用于维护更多线程同步的时间，已经超过了线程所运行的任务时间。在选择更大的岩石模型后，最短运行时间的极值将会在更多线程情况下出现。总之实验结果证明，在保证运行结果正确的前提下，本文所描述的并行算法使得岩石重构的运行时间大幅缩短。
小结
本文讨论了岩石重构模型使用模拟退火算法在并行计算时遇到的问题，对此分析，提出了一个解决方案，并实验验证。岩石重构问题与传统的NP完全问题，例如货郎担问题等在并行处理时，是不同的，不能直接从并行的解中选择其一。而是需要对这些并行的解进行整合叠加，之后产生新的解，才能有效提高运行效率。整合的过程需要对能量函数进行分析，只有满足特定的条件，才能使用本文所描述的并行处理方法。类似岩石重构这样的应用或许是一类典型的问题，我们将继续对这类问题进行深入的研究。
參考文献：
1.鞠杨，杨永明，宋振铎，徐文静. 岩石孔隙结构的统计模型. 中国科学，2008， 38（7）： 1026～1041
2.王金波，孔隙砂岩的三维重构模型研究，硕士学位论文，中国矿业大学（北京），2010
3.郑江韬，岩石孔隙结构三维分形重构，硕士学位论文，中国矿业大学（北京），2012
4.康立山，谢云，尤矢勇，罗祖华，非数值并行算法---模拟退火算法，科学出版社，1997
作者简介：
钱旭，男，1962年10月，中国矿业大学（北京），博士，教授，博士生导师。研究方向：数据库、信息融合技术、软件工程理论与技术、计算机支持的协同工作（CSCW）技术、知识工程。地址：北京市海淀区学院路丁11号 中国矿业大学（北京） 科技楼707A 黄耀辉收
关键词：模拟退火、并行、三维重构
【分类号】：TM76；TM744
引言
岩石结构中含有大量不同尺寸、随机分布的孔隙，这些孔隙直接影响着岩石的物理力学性质。如果能够定量描述岩石的微细观结构并找到这些微细观结构与岩石宏观力学行为之间的关系，将为科学、准确地描述工程岩体和工程结构的响应，以及预测、控制水利水电和土木工程中的灾害控制提供有力的依据[1]。岩石重构模型是在CT扫描实验的结果上，提取岩石孔隙结构特征，构建了表征岩石孔隙分布与相关性质的两点概率函数和线性路径函数，并以此为控制函数，运用模拟退火算法重构了岩石的三维孔隙结构。
模拟退火算法SA（simulated annealing algorithm）是适合解大规模组合优化问题，特别是NP完全问题的通用有效近似算法[4]。它源于对固体退火过程的模拟，固体退火是先将固体加热至熔化，再徐徐冷却。加温时固体内部粒子能量升高，变为无序状态，冷却时，粒子内能减少，趋于有序；系统在每一个温度下都会有一个趋于平衡态的过程，最后凝固成规整晶体的热力学过程[4]。
将模拟退火算法应用到重构模型上，是以真实的孔隙分布数据作为输入资料，建立了与真实孔隙结构具有相同特征的数字孔隙模型，所建模型能较为真实反映的岩石类材料的内部结构。但该方法的应用受到收敛速度慢的困扰，尤其是当模型尺寸增大，所需的计算时间大幅增加。因此需要找到办法提高模拟退火算法在岩石重构计算中的效率，并行计算正是我们期待的方法。
问题描述
岩石重构的过程是先用CT扫描实验的结果，计算岩石的孔隙率，并构建表征岩石孔隙分布与相关性质的两点概率函数和线性路径函数。根据孔隙率构建一个初始是随机分布的三维立方体模型，在这个立方体中每个个坐标位置的值是0或者1，用来表示空隙或岩石实体。这个初始的模型就是初始状态S，而两点概率函数和线性路径函数就是用于计算能量的控制函数F（S）[2][3]。模拟退火算法中的每一次实验，就是尝试交换这个三维立方体中的一对点（一个0点和一个1点），然后计算能量变化，并判断是否接受为新的状态。
本问题可归结为优化问题，即寻找最佳状态S，使下式最小
其中：
P（r）为参考模型的两点概率函数，P’（r）为当前状态的两点概率函数
L（r）为参考模型的线性路径函数，L’（r）为当前状态的线性路径函数
应用SA算法解决该问题，经典的模拟退火算法是串行的，其算法描述如下：
1 初始化：初始的足够高的温度T，初始状态S，每个温度下的实验次数L（即Markov链长度）。这相当于将固体加热至足够高。
2 循环执行L次实验，每个实验执行以下四个任务。这相当于在一个温度下系统趋于热平衡的过程。
2.1 对S进行“扰动”，计算能量，产生新解S’
2.2 计算S’和S的能量差ΔE
2.3 判断是否接受：如果ΔE<0，接受S’为新解，否则以概率exp（-ΔE/T）接受S’
2.4 接受或拒绝新解
3 根据温度或其他条件判断是否结束，结束则输出当前S，否则降温并继续执行第二步。
并行计算
在将模拟退火算法并行的实践中主要有以下两种策略：
1. 并行操作[4]：把每个任务中的多个操作做并行处理。这个策略受到具体问题的限制，可并行完成的操作数每个问题不同。处理机个数增加到一定数量后，算法性能不再提升。
2. 并行实验[4]：在同一个温度下，n个处理机同时进行实验，产生n个新解S’1， S’2 … S’n ，从中找到能量降低的m个解S’’1， S’’2 … S’’m，然后按某个准则接受其中之一作为新解。这实际上已经与原有模拟退火算法不同：模拟退火算法中每次产生的新解S’是对前一个解S的扰动产生的，然后判断是否接受；而并行实验中的所有S’都是对同一个S扰动产生的，除了要用Metropolis接受准则判断外，还要有一个准则用来选择新解。这个策略可以充分发挥并行计算的优势。像是货郎担问题（Travelling Salesman Problem， TSP）、0-1背包问题（Zero-One Knapsack Problem， ZKP）、调度问题（Scheduling Problem， SCP）等所谓的NP完全问题都可以采用这个策略。
在岩石重构的过程中，生成新的状态并计算其能量的任务难以拆分为多个可以并行的操作，并且并行操作也难以发挥更多处理机的性能。所以我们希望用并行实验的方法来提升模拟退火算法计算岩石重构的效率。
在将并行实验策略应用到岩石重构的退火算法上时，每次并行实验在当前解S的基础上产生n个新解S’1， S’2 … S’n，从中找到能量降低的m个解S’’1， S’’2 … S’’m，根据并行实验的策略，这时应该从m个解中选择一个作为新的解。但不论是哪一个解，都只是在S的基础上交换了两个点，从整体上看这样并行与串行计算的效率几乎相同，甚至由于需要考虑并行线程之间的同步，效率反而更低。
如果能将能量降低的m个解所交换的m对点直接应用在S上而形成新的解，则在原有串行算法中交换一对点的时间内，并行算法理论上可以交换m对点，这显然是我们期待的结果。问题的关键是，虽然已知在S上交換m对点中的任何一对都会使能量下降，但是同时交换全部的m对点能量是否会下降，即m对点之间是否会有干扰。如果没有干扰，则整体能量一定会下降。 所以需要分析能量差ΔE与交换的点之间的关系：
设S为当前解，A点为S中的一个0点，B点为S中的一个1点，A的坐标为（Ax， Ay， Az），B的坐标为（Bx， By， Bz），S’为交换A、B两点后的解 ，ΔE为S和S’的能量差，即
ΔE = F（S’） – F（S）
通过对ΔE的方程进行分析，发现如下情况，ΔE只与A、B两点所在的x方向的两条线（A点一条线、B点一条线）、y方向的两条线，z方向的两条线，一共是六条线上的所有的点的值相关。也就是说，前述的m对点，如果他们的坐标都不相同，则相互之间对于能量没有干扰。
据此，我们提出岩石重构模型的并行模拟退火算法如下：
1 初始化：初始的足够高的温度T，初始状态S，每个温度下的并行实验次数L，每次并行实验的次数N（线程数）。
2 循环执行L次并行实验，每次并行实验同时有N个单独实验。设单独实验的编号为x∈（1，N）。
2.1 选取N对互不干扰的点，分别交给每个单独的实验（主线程做）
2.2 第x个实验做以下工作（每线程单独做）：
2.2.1 在S上交换一对点，并计算能量，产生临时解S’x
2.2.2 计算S’和S的能量差ΔEx。
2.2.3 判断是否接受为临时解：如果ΔEx<0，接受S’x为临时解，否则以概率exp（-ΔEx/T）接受S’x为临时解。
2.3 汇总临时解（主线程做）：设有m个被接受的临时解，即有m对可交换点，交换S中的这m对点，得到新解S’
3. 根据温度或其他条件判断是否结束，结束则输出当前S，否则降温并继续执行第二步。
仿真结果
实验使用的硬件条件如下：AMD Opteron Processor 6378 2.4G 4处理器 512G内存
岩石重构实验中我们采用如下参数：立方体为50x50x50像素，共计125000个点，降温50次，每温度尝试交换100000对点。
岩石重构模拟退火算法中，不论是串行还是并行，其初始状态都是类似的，下面为初始状态z轴方向的截图：
Fig 1 串行计算初始状态
Fig 2 并行6线程计算初始状态
重构计算完成后的结果的z轴方向截图：
Fig 3 串行计算最终结果
Fig 4 并行6线程计算最终结果
从结果图中可以看出，串行和并行计算后的结果满足统计的数据，是可以接受的结果。
下面是串行和并行多线程程序的运行时间：
Fig 5计算时间与线程个数关系图
从上面的结果看，使用并行算法后，随着并行线程数增加，运行时间在减少，在6线程时达到极值，之后随着线程增加运行时间反而增多，这主要是由于我们实验的岩石模型选的太小，用于维护更多线程同步的时间，已经超过了线程所运行的任务时间。在选择更大的岩石模型后，最短运行时间的极值将会在更多线程情况下出现。总之实验结果证明，在保证运行结果正确的前提下，本文所描述的并行算法使得岩石重构的运行时间大幅缩短。
小结
本文讨论了岩石重构模型使用模拟退火算法在并行计算时遇到的问题，对此分析，提出了一个解决方案，并实验验证。岩石重构问题与传统的NP完全问题，例如货郎担问题等在并行处理时，是不同的，不能直接从并行的解中选择其一。而是需要对这些并行的解进行整合叠加，之后产生新的解，才能有效提高运行效率。整合的过程需要对能量函数进行分析，只有满足特定的条件，才能使用本文所描述的并行处理方法。类似岩石重构这样的应用或许是一类典型的问题，我们将继续对这类问题进行深入的研究。
參考文献：
1.鞠杨，杨永明，宋振铎，徐文静. 岩石孔隙结构的统计模型. 中国科学，2008， 38（7）： 1026～1041
2.王金波，孔隙砂岩的三维重构模型研究，硕士学位论文，中国矿业大学（北京），2010
3.郑江韬，岩石孔隙结构三维分形重构，硕士学位论文，中国矿业大学（北京），2012
4.康立山，谢云，尤矢勇，罗祖华，非数值并行算法---模拟退火算法，科学出版社，1997
作者简介：
钱旭，男，1962年10月，中国矿业大学（北京），博士，教授，博士生导师。研究方向：数据库、信息融合技术、软件工程理论与技术、计算机支持的协同工作（CSCW）技术、知识工程。地址：北京市海淀区学院路丁11号 中国矿业大学（北京） 科技楼707A 黄耀辉收