反演算法（精选10篇）
反演算法 篇1
1 数据准备
1. 影像准备 :1通过美国地质勘探局网站下载研究区域灞河流域的遥感影像图,下载遥感图为Landsat8影像。2经定位分析得知,屯溪流域位于灞河流域位于陕西省蓝田县灞源乡华岔村西部。查阅该地区TM遥感影像行列号,确定屯溪流域Landsat 8遥感影像所在的行列带号。3按影像行列带号下载影像及头文件。
2. 模板准备 : 准备灞河流域矢量模板,用于对影像进行裁剪获得灞河流域影像图,矢量模板文件为 :灞河裁剪模板 .evf
2 温度反演
Landsa8影像的10波段是热红外波段,接受的是与地表温度高低相对应的地物的热红外辐射,因此,可以通过其来反演出地表温度 , 在温度反演之前需要进行图像拼接与裁剪。
2.1图像拼接和裁剪 :三张原始影像的第10波段影像,进行拼接得到拼接后影像后,打开灞河EVF文件,进行掩膜处理,最后把掩膜的文件放在拼接图像上进行裁剪,得到裁剪结果。
2.2各类算法温度反演
2.2.1 B算法
这种方法简单,易于操作,对影像之外的外来参数依赖性很小。反演效果也较好,是种简单实用的方法。
2.2.2单通道算法
2.2.3单窗算法
2.2.4辐射传输方程法
大气平均作用温度估计误差对地表温度演算精度有较大影响 , 但这种影响直接取决于大气透射率的大小。大气透射率估计误差对数据的地表温度演算精度影响很大,这一误差随地表温度增加而迅速增大 , 温度越高 , 误差越大。
2. 地表比辐射率计算 :水体像元的比辐射率赋值为0.995,自然表面和城镇像元的比辐射率估算则分别根据下式进行计算 :
b1为NDVI值 ;b2为植被覆盖度值,可以通过计算得到地表比辐射率数据。
在NASA官网 (http://atmcorr.gsfc.nasa.gov/) 中输入成影时间以及中心经纬度,可以查询上面公式中所需的参数。利用ENVI主菜单 ->Basic Tools->Band Math,在公式输入栏中输入 :
4. 反演地表温度 :在获取温度为TS的黑体在热红外波段的辐射亮度后,根据普朗克公式的反函数,求得地表真实温度TS :
b1 :公式中的温度为T的黑体在热红外波段的辐射亮度值。由此可以得到真实的地表温度值(摄氏度)。
3 四种反演方法比较
见图1
B算法反演过程较为简便,易操作 ;对影像数据之外的外来参数依赖性较小 ;考虑了地表比辐射率因素的影响 ;单通道算法分类讨论地表比辐射率的影响,对结果反演较好 ;单窗算法既考虑了地表比辐射率的影响,也考虑了大气辐射的影响 ;反演过程所需要的大气参数比传统的辐射传导方程法要少得多,仅需要近地表气温和大气水分含量这两个参数 ; 辐射传输方程法也较为简便,易操作 ;分类讨论了地表比辐射率的影响。
4 论述
单窗算法反演和大气校正法的地表温度最高有54度,误差很大,与其他两者相比,都高出8度左右 ;
1) 总体比较 ;B算法和辐射传输方程法反演出来的温度差异在3度左右 ;两者与单窗算法反演出来的温度相差5-7度左右 ;总体反演温度趋势呈 :单通道算法
2) 各类别比较
见图2
3) 综述 ; 综上可看出,单通道算法反演的温度普遍偏低,但与其他算法最大不超过1度,单窗算法反演温度普遍偏高,比其他算法高出8度左右 ;B算法居中,辐射传输方程法在高温地区反演效果较差,与其他两者相差较多。
摘要：随着科学技术的发展,遥感技术已经成为了对地观测的重要手段,主要的反演地表温度的方法有单窗算法、单通道算法、B算子、辐射传输方程等。本文采用landsat8数据,对灞河流域进行研究分析并比较四种方法。
反演算法 篇2
基于混合优化算法的水声环境参数反演
在水声环境参数反演中,优化算法对反演结果有重要的影响.提出用遗传算法和Levenberg-Marquardt算法组成的`混合优化算法进行浅海环境参数匹配场反演,利用仿真数据对算法的有效性进行了评估,并与遗传算法、模拟退火算法进行了比较,最后采用实验数据验证了方法的有效性.
作 者：杨坤德 朱明永 张鹏 YANG Kun-de ZHU Ming-yong ZHANG Peng 作者单位：西北工业大学航海学院,陕西,西安,710072 刊 名：电声技术 ISTIC英文刊名：AUDIO ENGINEERING 年，卷(期)：2009 33(2) 分类号：O422.2 关键词：匹配场反演 遗传算法 模拟退火算法 混合优化算法关系映射反演方法综述 篇3
【关键字】关系映射反演方法，RMI
一、引言
关系（relationship），映射（mapping）反演（inversion））方法是我国著名数学教育家徐利治先生提出来的。这个方法提出伊始就得到了广泛的重视和应用。关系映射反演方法是化归原则在数学中具体体现，这是一种更强的结构模式。运用关系映射反演方法需要充分地联想和类比，把一事物翻译为另一事物，得到另一事物的解答，最后再把这个解答翻译为对第一个事物的解答，从而迂回地解决问题。这样往往能进入“山穷水复疑无路，柳暗花明又一村”的思维境界。
二、关系映射反演方法
根据徐利治先生的观点，数学中的关系映射反演方法可一般的表述如下：给定一个含有目标原象的关系结构，如果能找到一个可定映映射，将映入或映满，则可从通过一定的数学方法把目标映像确定出来，进而，通过反演又可以把确定出来，这样，原来的问题就得到了解决.整个过程包括五个步骤：关系——映射——定映——反演——得解。[1]
这五个步骤的具体解释是：首先弄清问题原像关系结构和原像未知目标的具体内容；选择适当有效的映射；接着确定未知元素的映像；然后根据被确定了的映像目标通过反演而确定原像目标；最后得到原问题的解答。[5]
三、RMI原则在数学中的运用实例
（一）运用RMI方法的经典实例
例1：哥尼斯堡七巧问题。
哥尼斯堡七巧问题又叫做欧拉七桥问题。解决七桥问题，欧拉运用的就是RMI方法。令表示七桥问题中桥与岛及陆地之间的关系结构，为一次能否走过七座桥的问题。欧拉采用这样的映射：把桥对应为几何线，把联结地点对应为几何点。
原来的问题便对应为能否一笔画出上述平面图的问题。换句话说，便是关于上述点线图的一笔画问题。解此问题需要采用简单的逻辑推理，过程如下：
凡是一笔画中间出现的交点处，曲线一进一出总是通过偶数条，故均可称为“偶点”，只有作为起点和终点的两个点有可能成为“奇点”（即通过的曲线为奇数条）。所以凡是多于两个奇点的平面图都是不可能一笔画出来的。现今图形结构中四个交点都是奇点，因此它是不可能一笔画出的。这就是说，问题答案的是“一笔画是不可能的”。由此对应的反演回去（利用），便可知道原来问题的答案是：不可能不重复地一次通过这七座桥。
（二）笛卡尔发明解析几何。
当一个几何问题不便从几何本身的角度进行求解时，首先把这个几何问题解析表示为一个代数问题，其次在代数领域里求解这个问题得到代数解答，最后把这个代数解答再几何解释为几何结论，这个结论就是原几何问题的解答。
（三）RMI在初等数学中的运用实例
1.函数法。在中学数学范围内， 把一些待处理问题映射为初等函数。在初等函数的映像关系结构中对问题函数进行函数处理，得到函数结论，再利用反函数或函数性质进行反演，使原问题获解的映射反演方法就是函数法。
例1：为何值时，不等式恰好有一解。
分析：如果单纯用不等式的解法去做，会显得很繁琐。我们可以函数的图像，映射为具体的几何图像考虑。
解：令，则可知该函数图像是一条开口向上的抛物线。当与时，该函数图像时平行于横轴的两条直线。
当抛物线顶点在直线的下方时，不等式有无数解。
当抛物线顶点在直线的上方时，不等式有无解。
当抛物线顶点在直线的，不等式恰有一解。
再求其頂点坐标
故当时，原不等式恰有一解。
2.换元法。我们在研究某些复杂问题时， 通过引人一个或几个新变量来代替原式中的某些量，从而把原式用新变量表示，并求得相应的结果，这种解决问题的方法叫做换元法。
3.参数法。当问题中含有多个变量时，引进一个或几个变量，通过中间变量，把间题转化为参数问题， 进而再消去参数，使问题得到解决，这种处理问题的方法叫做参数法。例如求轨迹问题，用此方法比较简单。
四、RMI原则的教与学问题
一般在中学或者大学数学教材中，很少把RMI原则总结出来，所以很多人学了很多年数学或者教了很多年数学，都未必意识到他们所接触到的许多数学题材已经包含着RMI的方法与内容。但事实上，无论是初等数学还是高等数学中，都有不同水准的RMI方法或原则被隐含在其中，不过只有经过分析观察，才能把它抽象出来，并且对它包含的各个具体步骤给以确切的表述和讨论。所以，作为数学教师，要想教会学生们掌握好RMI方法或原则，首要的一步就是要采取“关系结构”的观点去考察数学问题、分析数学教材，并能从其中把联结原象与映像的映射关系揭示出来。当然，这是要依靠教师进行数学教学研究之后，才能很好地去完成这个任务的。
徐利治先生写到：“对于数学科学的学习者或者准备从事数学研究的人来说，应该以培养寻求映射的能力为目标。为了培养和提高理解原象系统（或现实原型）的能力，除了学好数学本科各分支之外，还需要学习自然科学，工程科学等有关分支领域的知识，要有较宽广的科技知识修养。事实上，能否很好地理解原象系统结构或某些应用科技中的现实原型，是决定能否正确地运用RMI方法的首要一步。所以要成为能真正解决问题的数学工作者，这方面的理解和洞察能力的锻炼以及与之有关的知识修养都是必不可少的。”[1]
参考文献：
1.徐利治，郑毓信，《关系映射反演方法》，江苏教育出版社，1988.
2.华东师范大学数学系，《数学分析》（第三版），高等教育出版社.
3.徐利治，数学方法论十二讲，大连理工大学出版社，2007（11）.
4.高兴佑，向长福，关系映射反演方法例谈，高校理科研究
基于遗传算法的土壤结构反演模型 篇4
1 电阻率的测量
测量土壤的电阻率有两种办法, 分别是三极法和四极法[1]。其中三极法是常常应用计算垂直接地棒的电阻公式来反推土壤电阻率的方法, 常用的公式为
式中:a为半径, m;l为长度, m;R为垂直接地极的接地电阻, Ω。上述两个公式分别是根据垂直接地棒的两种计算公式推导而来。
这种方法测量土壤电阻率有一个极大的弊端:在打垂直接地棒时, 土壤的松紧影响它和电棒之间是否紧致, 导致了电棒和土壤之间的接地电阻R的变化。接下来介绍四极法, 这个方法可以消除三级法进行测量时会产生的误差。
F.Wenner在1915年时发明一种新的更加正确的测量土壤电阻率的方法, 他是用布置在同一水平上的4个电极来测量的, 简称为四极法。具体布置方法如图1所示。
在图1中, 4个电极中2、3电极代表电压极, 1、4电极代表电流极, h0代表电极的埋深, a代表极间距。
均匀土壤的计算公式如下[2]:
当埋深h0为零时, 采用式 (1) , 其结果是土壤电阻率的真实值, V23代表电位差, 它是在点电流源I、-I对电压极2、3电极发起的。当埋深不等于零时, 其计算公式如下:
在测得I和V23后, 由式 (2) 即可求出土壤电阻率为
虽然引进的电流线和电压线会造成一定的误差, 但是这样非常有利于实际中操作。因此, 本文为了有利于工程使用和减小误差, 还是把电压极和电流极打在土壤不深处, 并且取a≥10h0, 可以粗略看做h0等于零, 式 (3) 可由式 (1) 代替。
在分层土壤中按照这种方法求出的只是把一个分层土壤看做了一个均匀土壤, 而不再是土壤电阻率的真值, 测量极间距和土壤结构不同对应的视在电阻率也不同, 它与电流的流动范围有关:当电棒之间的距离不大时, 电流极流出的电流绝大部分都是在表层土壤中流动, 得出就是表层土壤电阻率。当电棒之间的距离越来越大, 越来越多的电流流向了深层土壤, 反映出的数据是反映了深层土壤参数, 根据电棒之间距离的变化得到的数据是反映了不同深度土壤的参数。因此, 在实际应用中, 会测量多组电棒距离, 也就是常说的测量深度。
本文测得的一系列数据是土壤的视在电阻率, 对这些数据要用数学方法、磁场理论得到正确的土壤结构。最后根据所推导公式编制计算机程序。
2 视在电阻率的推导
用四极法不同测量深度得到的数据是视在电阻率, 分层土壤中视在电阻率满足下式[3]:
式中, A1 (λ) 和B1 (λ) 与分层土壤参数有关, 因此函数中含有每层土壤的参数。当深度变为零时, 式 (4) 变为[4]
用上述的复镜像法将式 (3) 、式 (4) 进行展开, 这样得到一个比较简单的公式[5]:
式中:ci、di分别为复镜像源的位置和大小;h0为测量电极的埋深;a为极间距。
3 建立视在电阻率的目标函数
为了找到一组正确的土壤结构参数, 实际中要测量m个视在电阻率的数值, 最后总会有一组数据满足下式:
式中:ρci为测量深度ai时用复像法计算的计算值;ρmi为测量深度ai时的测量值。x'为列向量, 其表达式如下:
4 目标函数进行最优化
遗传算法在20世纪80年代中期得到了蓬勃发展[6]。遗传算法的灵感来自于生物进化论, 在自然选择的模式下不断生成和不断的检验。遗传算法是在编码群体进化的基础上, 根据适应度, 使遗传机制和生物进化选择得以实现。在实现过程中, 一一迭代个体位串, 通过随机选择把位串中重要的基因提出来重新组合, 使新一代的群体个体比上一代要优化, 就这样不断的群体优化, 逐渐靠近最优化, 最终实现所要解决的问题[7]。适应度函数的设计、参数编码、遗传操作的设计、控制参数的设计和遗传参数的设定、参数编码称为遗传算法五要素。它是一个最适应环境的个体, 在一代一代不断进化中产生, 得到问题最优解, 但没有显示完全随机搜索的特点, 在整个进化过程中遗传算法的遗传操作是随机的。
遗传算法在天然气管道控制优化应用中, Goldberg提出了遗传算法的工作流程和结构形式, 通常称为标准遗传算法 (standard GA-SGA) 或规范遗传算法 (canonical GA-GA) 。人们通常把问题的特征和领域知识用到GA的应用过程中, 各种特定的GA便由此形成, 这些GA具有全局搜索以及求解不同种类的优化问题的能力[8]。标准的遗传算法的基本流程图以及结构如图2所示。
由图2可知, 遗传算法是要经过多次迭代才能得出结果, 其运行过程的基本内容如下:
1) 选择合适的算法编码策略, 用结构空间来替代参数集合。
2) 定义算法的适应度函数f (x) 。
3) 确定遗传参数如变异的概率Pm、交叉的概率Pc等;确定遗传算法的策略, 通过选择算法参数如群体的变异、交叉、选择方法、群体的大小n来确定算法策略。
4) 采用随机初始化的方法生成群体P。
5) 计算群体串位解码后的个体适应值f (x) 。
6) 根据已选择的策略, 在群体中用将交叉算子、变异、选择等方法进行繁殖产生下一代。
7) 如果群体性能以及已完成迭代次数不满足指标, 返回步骤6) , 就可修改遗传算法策略后再返回步骤6) 。
选取随机变异算子、单点配对交叉算子、赌盘选择算子为遗传算法策略的3个基本算子[9], 用浮点数编码方法作为编码方式[10]基本过程, 如图3所示。
5 计算实例并验证
算例一
对于双层土壤, 表1为一组实际工程中用四极法测的在不同测量深度下的视在电阻率数据。表2为根据本文的反演模型计算出的水平分层土壤的模型, 并与CDEGS的结果进行了对比。算法的遗传参数选取:群体的大小为200, 最大代数为500, 交叉概率为0.7, 变异概率为5%。
算例二
新疆国投哈密电厂一期工程在拟建场地共布置土壤电阻率测试点25个, 测试点均布置在工程钻孔处, 且测试点编号与钻孔编号一致, 主要选择在构架、主厂房、烟囱等主要建筑物地段。
依据《电力工程物探技术规程》 (DL/T5159-2002) , 本次土壤电阻率测定采用了四极电测深法, 测定 (解释) 最大深度为30 m。为了保证解释深度满足规程要求, 四极电测深法各极距的选定:测量极距MN/2为0.5 m;供大极距AB/2为1.5 m、3 m、5 m、7 m、9 m、13 m、20 m、45 m。
哈密电厂实测视在电阻率值1如表3所示。
哈密电厂实测视在电阻率值2如表4所示。
本次土壤电阻率测定使用重庆地质仪器厂生产的DDC-8型电子自动补偿仪。对现场使用的仪器设备, 按规程质量体系文件的要求, 使用前进行了认真的检查, 各项指标均满足测试要求。
根据表3和表4的数据, 用本文理论编制的计算程序计算得到的结果, 如表5所示。在设计该项目时取上层土壤电阻率5125Ω·m, 上层土壤厚度2.5 m, 下层土壤电阻率41Ω·m。
6 结论
1) 土壤分层模型在以往的接地参数计算中近似的看作均匀土壤会导致很大的误差, 因为大接地网, 电流会向深层土壤中流动。
2) 反演模型计算出的水平分层土壤的模型, 与CDEGS的测量结果相符, 验证了本文所建模型的正确性。
参考文献
[1]鲁志伟, 文习山.三角法测量双层土壤接地电阻的误差分析[J].高电压技术, 2002, 28 (7) :15-16.LU Zhiwei, WEN Xishan.Error analysis of grounding resistance of double-layer soil based on triangulation measurement[J].High Voltage Engineering, 2002, 28 (7) :15-16.
[2]鞠勇, 张波, 崔翔, 等.双层土壤中接地网接地电阻补偿测量法的修正曲线[J].电网技术, 2003, 27 (1) :44-46.JU Yong, ZHANG Bo, CUI Xiang, et al.The correction curve of grounding resistance compensation measurement in two-layer soil[J].Power System Technology, 2003, 27 (1) :44-46.
[3]解广润.电力系统接地技术[M].北京:中国电力出版社, 1999:185-205.XIE Guangrun.Power system grounding technology[M].Beijing:Chinese Power Press, 1999:185-205.
[4]TAKEHIKO T.calculation of earth resistivity for a deep-driven rod in a multi-layer earth structure[J], IEEE Transaction on Power Delivery, 1991, 6 (2) :608-614.
[5]张波, 崔翔.复镜像法中的的一种自适应采样方法[J].华北电力大学学报, 2002, 29 (4) :1-4.ZHANG Bo, CUI Xiang.An adaptive method of complex image sampling method[J].Journal of North China Electric Power University, 2002, 29 (4) :1-4.
[6]马永杰, 云文霞.遗传算法研究进展[J].计算机应用研究, 2012, 29 (4) :1210-1206.MA Yongjie, YUN Wen Xia.The progress of study on Genetic Algorithm[J].Application Research of Computers, 2012, 29 (4) :1210-1206.
[7]赵亮, 雎刚, 吕剑虹.一种改进的遗传多目标优化算法及其应用[J].中国电机工程学报, 2008, 28 (2) :96-102.ZHAO Liang, JU Gang, LJianhong.An improved multi-objective optimization with genetic algorithm and its application[J].Proceedings of the CSEE, 2008, 28 (2) :96-102.
[8]杨金明, 吴捷, 钟丹虹.多目标优化问题中一种改进的遗传算法[J].华南理工大学学报 (自然科学版) , 2001, 29 (12) :65-67.YANG Jinming, WU Jie, ZHONG Danhong.Multi-objective optimization of an improved genetic algorithm[J].Journal of South China University of Technology (Natural science edition) , 2001, 29 (12) :65-67.
[9]崔明义.一种基于景观特征的浮点数编码遗传算法的研究[J].计算机应用研究, 2007 (8) :48-51.CUI Mingyi.A study of the application of computer float coded genetic algorithm based on the landscape characteristics[J].Application Research of Computer, 2007 (8) :48-51.
反演算法 篇5
基于MPI并行实数编码混合遗传算法的波阻抗反演
地球物理反演的`局部线性方法易使解陷入局部极值,并严重依赖初始模型,而传统的遗传算法在优化应用中存在局部搜索能力弱、早熟收敛等问题,为此,提出了一种解决地球物理反演问题的并行实数编码混合遗传算法(MRCGA).该方法采用拟网格法初始种群、综合交叉策略和线性算子,实现了并行实码混合遗传算法.理论模型试算证明了该算法反演地震波阻抗的有效性.
作 者：白俊雨 赵俊省 潘凌飞 郭嵩魏 Bai Junyu Zhao Junsheng Pan Lingfei Guo Songwei 作者单位：成都理工大学地球探测与信息技术教育部重点实验室,四川,成都610059 刊 名：勘探地球物理进展 英文刊名：PROGRESS IN EXPLORATION GEOPHYSICS 年，卷(期)：2009 32(4) 分类号：P631.4 关键词：MPI 遗传算法 实数编码 初始种群 综合交叉策略 线性算子 波阻抗反演反演算法 篇6
通过重力异常反演地球内部横向密度的不均匀性是了解地球内部的精细结构的方法之一，而反演方法很多，线性反演方法主要有最小二乘法、脊回归法等。非线性有遗传算法、退火算法等，杨长福(2004年)将以单元密度和网格厚度为模型参数，进行重力异常的多层密度及其界面的同时反演。
目前反演多层密度界面的方法主要有两种：一是把多层密度界面的重力场，通过分离场的方法逐层分离，逐层反演;二是直接反演，如王一新的正则化反演法和陈胜早的N层变密度反演法。这些方法都很复杂，且要求密度分布已知，使用起来很不方便。蚁群算法(Ant Colony Algorithm)是人们对蚁群集体行为的研究成果，而提出的一种基于种群的模拟进化算法，属于随机搜索算法。该算法思想由意大利学者M.Dorigo等人首先提出，充分利用了蚁群搜索食物的过程与著名的旅行商问题(TSP)之间的相似性，通过人工模拟蚂蚁搜索食物的过程来求解TSP。杨宁(2006年)用蚁群算法解算DGPS整周模糊度，彭珊鸽(2007年)用蚁群算法来优化点状注记的配置问题，刘宁(2008年)将蚁群算法应用反演断层参数，本文将以单元密度为模型参数，利用蚁群算法进行重力异常的密度反演。
1 正演模型及计算方法
在实际研究地质体所在的空间时，一般将区域进行剖分，如图1所示，z向分为多个水平层，每层按照对x，y向分成若干个规则长方体，每个长方体大小形状相同，长方体各边分别平行于三个坐标轴，对于每个长方体可以认为许多个点源组成的，在单一的长方体中我们认为密度是均匀的，用积分的方法可以得到长方体产生的重力异常。
现介绍长方体外任意一点的重力异常的计算，取图1中第k个长方体并置于图2所建立坐标系中，设q(ε，η，τ)为长方体内任一点坐标，其长方体的空间取值范围为：ε∈[a1，a2]，η∈[b1，b2]，τ∈[h1，h2]，ρk为长方体的剩余密度，根据引力位公式可以计算出长方体外任一点坐标p(x，y，z)的引力位。
其中：，G为万有引力常数，M为长方体的质量。
通过对上式求导，进行积分后的p(x，y，z)的重力异常
由重力的可加性，则所有长方体在p(x，y，z)点产生的总重力异常可表示为
2 蚁群算法
蚁群算法的思想最先用于求解旅行商问题(TSP)，下面就求解平面上n个城市的TSP问题(1，2，3…表示城市的序号)为例来说明该算法思想。假设m只蚂蚁放入n个随机选择的城市中;dij为城市i、j之间的距离;蚂蚁在运动过程中，根据各条路径上的信息素选择下一个它没有访问的城市，同时完成一步(或者完成一个循环)对信息素进行更新，选择下一个城市的主要依据有：τij(t)为t时刻连接城市i和j的路径(i，j)上的残留信息量，初始时刻各路径上信息量相等，设τij(0)=C(C为常数);ηij表示城市i转移到城市j的期望程度，可根据某种启发式算法具体确定，在TSP问题中一般取ηij=1/dij。蚂蚁k(k=1，2，…，m)根据各条路径上的信息量决定转移方向，t时刻蚂蚁k从城市i向城市j转移的概率Pkij(t)计算如下：
其中，allowed={0，1，…，n-1}-tabuk表示蚂蚁k下一步允许选择的城市。α-残留信息量的相对重要程度;β-期望值的相对重要程度;与自然蚁群系统不同之处在于人工蚁群系统具有一定的记忆力，tabuk(k=1，2，…，m)用于记录蚂蚁k所走过的城市，集合tabuk随着进化过程进行动态调整。人工蚁群保留了自然蚁群信息素挥发特点，随着时间的推移，以前留下的信息逐渐消逝，参数ρ(0≤ρ<1)表示信息素的持久性，1-ρ则表示信息素的衰减度。在每只蚂蚁完成对所有n个城市的访问后(即一次循环结束)，各路经商的信息素量根据公式(5)(6)进行调整：
其中，Δτkij表示第k只蚂蚁在本次循环中留的路径(i，j)上的信息素量，Δτij表示本次循环中路径(i，j)上的信息素的增量。
其中，Q是一个常数，表示蚂蚁所留的信息素量，Lk表示第k只蚂蚁在本次循环中所走的路径长度。在初始时刻Δτij=0(i，j=0，1，…，n-1)，ηij表示由城市i转到j的期望程度，可根据启发算法具体确定。M.Dorigo曾给出三种不同的实现方法，由于ant-quantity system、ant-density system利用的是局部信息，而ant-cycle system利用的则是整体信息，因此在求解TSP时ant-cycle system性能较好，因此通常采用它为基本模型。
3 利用蚁群算法反演密度
3．1目标函数
从本文前面论述(3)式可以看出，在给定长方体参数的情况下，根据已知重力异常确定剩余密度参数时，也就是在每一个独立观测方程含n个未知参数ρ(ρ1，ρ2，ρ3，…，ρn)，其解具有多值性，为寻求最优解，利用最小二范数准则建立目标函数为：
式中，N为观测值个数，gj(ρ)为第j个观测点上由假设密度参数ρ计算所得的重力异常，gj(T)为理论重力异常值，Φ(ρ)为目标函数。
3．2反演密度
由于最初的蚁群算法是源于离散型网络是的路径问题，因此在一般函数优化问题中，必须对原有模型许多细节进行修正。
设一组密度参数ρ(ρ1，ρ2，ρ3，…，ρn)有n个变量，其中ρi的取值范围为[Ri-1，Ri+1]，若把每个密度参数的取值范围分成S份，那么密度参数的组合情况就有Sn种。我们就要在Sn种参数组合中找到一种最佳参数组合(最短路径)ρ，使目标函数Φ(ρ)的值最小。在TSP问题的目的是找一条最短的闭合路径，我们可以发现在进行信息素更新时，信息素由路径的长度确定。而在密度参数反演问题中，蚂蚁所选择的路径长度，我们利用目标函数值大小来确定，定义蚂蚁k在某次循环中所选择路径的长度Lk为Lk=Φwk(P)，其中Φwk(P)表示该次循环中，蚂蚁k选择的路径w的目标函数值。
4 试验数据模拟
在如图3中坐标系是图2坐标系的水平投影，9个长方体的长、宽、高分别取1km，对于各个长方体的密度分布如图3所示，单位为：10-3 kg/m3。在本次试验中各观测点为长方体上表面矩形的中心位置，由公式(3)可以计算观测点的重力异常值。在本例中取所有参数在[0.4，1.5]范围之内，将每个密度参数取值范围分成10000份，然后采用随机函数产生一组参数，利用(8)式求出目标函数，在由(4)、(5)、(6)、(7)式建立蚁群算法的递推关系并且进行密度参数的更新。
5 结论
(1)蚁群算法经过改进能够用于已知重力异常的密度反演。
(2)密度反演结果与与理论数值上有一定的差异;这种差异同样与给定的搜索空间和空间划分的精细程度有关。
参考文献
[1]柯小平, 王勇, 许厚泽.用遗传算法反演地壳的变密度模型[J].武汉大学学报信息科学版2004, 29 (11) :981-984.
[2]刘宁, 张永志.位错模式的蚁群算法反演断层参数[J].大地测量与地球动力学, 2008, 28 (6) :1-6.
[3]杨长福.用脊回法反演重力异常的多层密度及其界面[J].西北地震学, 2004, 26 (4) :293-297.
[4]ColoniA, Dorigo M and ManiezzoV.Ant system:Optimization by a colony of cooPeratingagent[J].IEEE Trans onSystems, Man and Cyberneties-PartB:Cybemeties.1996, 26 (1) :29-41.
[5]侯遵泽, 杨文采, 刘家琦.中国大陆地壳密度差异多尺度反演[J].地球物理学报, 1998, 41 (5) :642-651.
[6]王家林, 王一新, 林桂康.利用场变换的分离场法反演多层密度界面[J].石油物探, 1986, 25 (2) :69-80.
[7]王一新, 王家林, 张曙明.研究多层密度界面的正则化非线性反演方法[J].石油物探, 1987, 26 (1) :78-90.
[8]陈胜早.用多层变密度模型直接反演法研究下扬子地区的岩石圈构造[J].地球物理学报, 1989, 32 (1) :20-35.
[9]杨宁, 田蔚风, 金志华, 蚁群算法在DGPS动态整周模糊度解算中的应用[J].航天控制, 2006, 24 (4) :4-7.
[10]彭珊鸽, 宋鹰, 吴凡.基于蚁群算法的点状注记智能化配置[J].测绘科学, 2007, 32 (5) :80-81.
[11]Dorigo M, Gambardella, Maria L.Ant colonies for the traveling salesman problem[J].Biosys-tems, 1997, 43 (2) :73-81.
反演算法 篇7
遥感探测技术可以获得大范围、时间连续的地表温度,能够准确反映地表温度的时空分布[5]。特别是热红外地表温度反演技术发展迅速,其产品已广泛用于气象、农业、生态等领域[6]。然而,现有的成熟算法只适用于能见度高的晴天条件,无法用于灰霾天气条件下地表温度。
霾是一种大量极细微的干尘粒等均匀地浮游在空中,使水平能见度小于10.0 km的空气普遍混浊的天气现象[7]。研究表明,霾中的物质成分除了细尘以外,还包括硫酸与硫酸盐、硝酸与硝酸盐、碳氢化合物、黑碳等粒子,因此其发生的频数与人类活动所造成的气溶胶污染密切相关[8,9]。在我国,霾天气主要分布在中东部和南部地区,西部和东北部地区相对较少,而且在同一区域内,大中城市的霾天气较乡村明显偏多,“浊岛”现象明显。霾天气的季节分布基本上是冬多夏少,这是因为我国大部分地区冬季受大陆高压控制,大气层结相对稳定,雨日雨量较少,易于低层大气中气溶胶粒子富集,并形成霾天气。研究表明灰霾影响地表温度反演精度,反演的均方根误差可达2.2 K[10]。近年来灰霾天气频发,为使卫星遥感技术能持续精确监测地表温度,需开展灰霾天气条件下地表温度遥感反演机理研究。
因此,现以FY-3A VIRR和MERSI传感器数据为主要数据源,考虑传感器热红外通道光谱特点,结合气溶胶光学厚度,改进已有的反演算法,研究适用于灰霾天气影响下风云卫星资料的高精度地表温度反演算法。
1 研究数据及模型
1.1 研究数据
用于地表温度反演的数据为FY-3A MERSI和VIRR L1级数据。风云三号卫星是我国第二代极轨气象卫星,于2008年5月27日从太原发射升空。其上同时搭载了可见光红外扫描辐射计(visible and infrared radiometer,VIRR)和中分辨率光谱成像仪(medium resolution spectral imager)MERSI。
VIRR具有10个探测通道,在星下点为1 km。VIRR第3、4、5通道为三个红外大气窗区窗口,可分别用来获取白天和夜晚的地表信息。VIRR热红外通道对下垫面辐射变化敏感,通道4和5(波长范围分别为10.3~11.3μm和11.5~12.5μm)相当于AVHRR的分裂窗通道4和5,主要用于地表特征参数的提取。
MERSI可以实现宽视场(2 000 km)和高分辨率(250 m)的全球云层和地表的监测,具备从可见光到热红外20个通道的探测能力,光谱范围为0.41~12.5μm,可探测大气、陆地、海洋的可见光反射辐射和热红外发射的辐射[11]。MERSI第5通道为热红外通道,其中心波长为11.25μm,空间分辨率为250 m,定标精度1 K,可以得到地表热辐射信息,并定量反演成地表温度。
1.2 研究模型———MODTRAN
MODTRAN是在LOWTRAN的基础上发展起来的,可用于计算大气透过率或指定大气路径下的辐射值。其有效光谱频率为0到50 000 cm-1(或大于0.2μm的波长范围),最大光谱分辨率为2 cm-1(LOWTRAN光谱分辨率为20 cm-1)。MODTRAN用吸收系数、谱线密度参数和平均线宽三个与温度相关的参数,计算分子透过率,将光谱范围分割成1cm-1,该步长下可联合Correlated-K和二向反射分布函数集使多次散射项更为精确。
MODTRAN中包含10种气溶胶模式,选择模式自带的春夏、秋冬两种气溶胶剖面和改变气溶胶视距,模拟出各种天气条件下的气溶胶透过率信息。对于不同地表,可通过改变MODTRAN输入参数来实现。在模式中输入大气模式、地表温度、地表比辐射率、大气总水汽含量、观测角度、气溶胶模式、水平气象视距以及传感器的通道响应函数等信息,就可以模拟出卫星观测数据,输出卫星接收到的辐射能。其输入数据包括:
(1)11组地表温度,中纬度冬季取值范围为[267.2 K,292.2 K],取样间隔为2.5 K。
(2)比辐射率参考蔡国印[12,13]等有关地表比辐射率的研究成果,取各类地物在FY-3 VIRR4、5和MERSI第5通道的平均比辐射率作为输入(如表1所示)。
(3)9组比辐射率扰动,取值范围为[-0.016,0.016],取样间隔为0.004。
(4)VIRR第4、5通道及MERSI第5通道的光谱响应函数。
(5)9组大气水汽含量,取值范围为[0.1,2.5],采样间隔为0.3 g/cm2。
(6)12组卫星观测天顶角,取值范围[0,55],采样间隔5°。
(7)25组气象水平视距,取值范围[1,25],采样间隔1 km。
(8)4种气溶胶模式,分别为非气溶胶模式、乡村气溶胶模式、城市气溶胶模式和沙尘气溶胶模式。
2 FY-3A资料分裂窗模型构建
2.1 改进分裂窗算法
改进分裂窗算法是Becker和Li基于Wan和Doizer等人[14]利用辐射传输方程LOWTRAN程序进行地表辐射亮温计算,发展出的分裂窗算法上经过订正得来的。Becker和Li(1995)针对NOAA传感器两热红外通道数据提出了半经验分裂窗模型,并且得到广泛运用。
根据Becker局地分裂窗算法,地表温度Ts(K)可以表达为两个热红外通道(i=1和i=2)亮温的函数:
式中,T1和T2是两个热红外通道1、2的亮温;ε1、ε2是两个热红外通道1、2的地表比辐射率;A0、P、M是局地分裂窗系数,A0是与大气状况有关的系数,P和M可以表达为两个通道的波段比辐射率ε1和ε2的函数;ε是两个热红外通道的平均比辐射;Δε是两个热红外通道的比辐射率之差,α、β、γ'、α'、β'是常数,可以通过最佳拟合获得。
早期的分裂窗算法用于海表温度反演,精度可以达到1 K以内。为了使地表温度反演的结果更加精确和合理,Becker等在原有的分裂窗算法模型中又加入了水汽和观测角度项。根据Becker and Li提出的分裂窗算法,地表温度Ts(K)可以表达为两个热红外通道i(i=4,5)亮温的函数式中,T4和T5分别为两个热红外通道4、5的亮温(K)。
式中,W(g/cm2)为大气水汽含量;θ为卫星观测角度;ε4、ε5分别为两个热红外通道4、5的地表比辐射率;A0、P、M都是分裂窗系数,A0是与大气状况有关的系数,P可以表达为卫星观测天顶角θ、水汽含量W以及地表比辐射率ε4、ε5的函数;M则可表示成ε4、ε5以及水汽W的函数。a1~a13为常数,可以通过基于最小二乘法的多元一次线性拟合得到。
2.2 气溶胶对模拟结果的影响
2.2.1 气溶胶对亮温的影响
为探究气溶胶对亮温的影响,统计了不同气溶胶模式和能见度条件下,组合通道之间的亮温差值ΔT随地表温度升高的变化趋势,结果如图1所示。图中气溶胶模式包括清洁大气、乡村气溶胶、城市气溶胶以及沙尘气溶胶,通道组合包括VIRR4、5通道组合(图中表示为V4&V5)以及VIRR4、MERSI5(图中表示为V4&M5)通道组合。基于模拟结果,发现相比于V4&V5,V4&M5的ΔT随地表温度变化不显著,其ΔT在任何水平气象视距(VIS)和地表温度条件下均保持2 K左右的差值。
而VIRR4、5通道之间的ΔT随地表温度和水平气象视距(VIS)的变化比较复杂,每种气溶胶模式展现了亮温差值变化的不同特征。在清洁大气中,当VIS=1 km时,两通道ΔT基本可以忽略不计,并且对地表温度敏感性较弱。当VIS=5 km时,ΔT达到0.5 K左右,但仍然不受地表温度的影响。当大气条件改善到VIS=15 km,ΔT开始对地表温度敏感,并且随地表温度的升高而升高。在乡村气溶胶模式下,ΔT基本大于0,当VIS低于10 km时,ΔT随地表温度变化不明显,基本上保持在0.5 K以下。当空气质量良好,VIS=15 km,ΔT变化巨大,在地表温度为280 K时取得极大值,ΔT=8 K。整体而言,VIS=15 K时,ΔT都在5 K以上,然而随地表温度变化没有明显的规律。在城市气溶胶类型,ΔT随温度升高而升高。其中,当VIS=10 km时,ΔT随温度升高经历了由负转正的反转关系。沙尘气溶胶条件下的ΔT关系与上述各类都有所差别,其ΔT在267.5 K时取得极大值,其中VIS=15 K时两通道亮温差异性最强,达到了7 K以上。在其他地表温度区间内,VIS=15 km的ΔT都要大于其他气溶胶光学厚度。
综合所有通道组合和气溶胶类型,总体规律可总结为VIRR4通道亮温要高于VIRR5和MERSI5通道亮温,并且在同一气溶胶模式下,气溶胶含量越低,两通道之间的ΔT越大。
2.2.2 气溶胶对分裂窗算法的影响
气溶胶主要通过吸收地表长波辐射的方式影响卫星接收地表辐射能量。不同类型气溶胶粒子大小和分布浓度都有显著差异,所以气溶胶对分裂窗算法精度影响主要体现在气溶胶类型和气溶胶光学厚度这两方面。首先探究气溶胶类型对算法模型的影响,能见度都设为模式缺省值。分别以VIRR4/VIRR5和VIRR4/MERSI5为组合,构建分裂窗算法模型。分别在无气溶胶、乡村气溶胶、城市气溶胶和沙尘气溶胶这四种模型下进行拟合,与不考虑气溶胶类型所拟合的结果对比,最后的拟合结果如图2所示,不同能见度条件下模拟地表温度的均方根误差(RMSE)。蓝色表示将VIRR4、5通道作为分裂窗的组合通道,红色表示VIRR4、MERSI5组合通道,每种气溶胶类型下拟合结果的RMSE显示在图像横坐标下方。从图中可以明显看出,同样的大气条件下同一传感器两个热红外通道的反演精度要优于不同传感器的通道之间的组合,VIRR4、5通道组合平均的RMSE相比于VIRR4、MERSI5要低0.5 K。当处于不包含气溶胶的理想大气时,反演精度最高,两种通道组合拟合的结果显示RMSE均低于0.11 K。在考虑到气溶胶的存在时,拟合的精度有所降低。而我们不去考虑气溶胶的状况,将所有气象条件下的地表温度反演全用一个公式去表达,得到的RMSE是最大的。所以气溶胶对地表温度反演的影响应该被考虑到分裂窗算法中。
将各种气溶胶模式下模拟的地表温度结果与输入的地表温度进行对比,两者的误差区间和分布范围由图3显示。图中每行表示气溶胶模式,列代表选择的通道组合类型,(a)为VIRR4和VIRR5通道组合,(b)为VIRR4和MERSI5通道组合。这里,横坐标表示的是输入输出温度的差值,而纵坐标表示的是这些差值所对应的频数。当每种气溶胶模式分开模拟式,误差基本符合均值为0的正态分布,如果忽略气溶胶的影响,整体分布趋势向右倾斜,即模式输入值要大于模拟结果,且最大差值可达2 K。再将同一气溶胶模式下的不同通道组合对误差的影响进行比对。结果显示:VIRR4、5通道的拟合结果要明显优于VIRR4、MERSI5的组合结果,主要体现在误差离散型小,各气溶胶模式下的误差极值也要小于后者。
在确认不同气溶胶类型带来的反演误差后,进一步分析各类气溶胶条件下分裂窗算法模型反演误差随气溶胶光学厚度的变化趋势。现在分别模拟了VIRR4、5通道和VIRR4、MERSI5通道组合分别在无气溶胶、乡村气溶胶、城市气溶胶和沙尘气溶胶的大气状态下,当水平气象视距从1 km(浓雾天气)变化到25 km(清洁大气),分裂窗算法的拟合系数随水平气象视距的变化,如图4所示。图中显示,当VIS小于7 km,各项系数随VIS变化剧烈,且每种气溶胶模式都有类似的变化规律。这说明当VIS比较低,即气溶胶光学厚度比较大的大气条件下,气溶胶对反演模型的精度有巨大影响。
在不同气溶胶模式中,分别统计地表温度模拟结果在不同VIS条件下与MODTRAN模式预设地表温度之间的RMSE(均方根误差)、BIAS(平均绝对偏差)、R2(相关系数),统计结果如图5所示。图中用红色表示VIRR4、5通道组合结果,黑色表示为VIRR4、MERSI5通道组合结果。结果表明,RMSE、BIAS随VIS增长而减小,相关系数随VIS增大而提高,这种关系在VIS小于7 km的气象条件下尤为明显,在图像中表现为曲率明显高于7 km以上时以上三个精度评价因子随VIS变化的曲率。具体到各通道组合之间的对比,图5(a)、(b)中红线都位于黑线之上,这表示VIRR4、MERSI5拟合结果的RMSE和BIAS都要优于VIRR4、5的组合,这在图5(c)中表现为前者具有更高的相关系数。图5中的各类气溶胶模式用不同的散点形状表示,2、3、5、8分别代表着无气溶胶、乡村气溶胶、城市气溶胶和沙尘气溶胶。如图所示,城市气溶胶模式下拟合结果与模式输入结果差距最大。VIS=1 km时,RMSE可达1.8K,BIAS也在1 K左右。而变化最平缓的是无气溶胶模式,在VIS从1~25 km的过程中,整个RMSE变化小于0.5 K,BIAS变化小于0.4 K,R2都保持在0.99以上。
[(a)列表示VIRR4、5组合下系数a1~a13与水平气象视距关系;(b)列表示VIRR4、MERSI5组合下系数a1~a13与水平气象视距关系;行1:无气溶胶;行2:乡村气溶胶;行3:城市气溶胶;行4:沙尘气溶胶]
3 结果与分析
3.1 包含气溶胶信息分裂窗算法建立
经过上述模拟实验,可以发现在气象水平视距小于7 km时,气象水平视距和分裂窗拟合系数a1、a11、a12、a13有良好的幂函数关系,所以以气象水平视距为自变量,各系数为因变量,以幂函数为模型,采用最大似然法将两者进行拟合,拟合的结果如表2所示。表中coff表示幂函数的指数,R2表示采用利用幂函数建模后水品气象视距和上述四个系数之间的相关系数。结果显示变量之间相关系数都在0.9以上,所以可用水平气象视距表示上述系数。在此关系上建立包含气溶胶信息的新分裂窗算法模型。
改进后的算法模型可表现成如下的形式
式中,W为大气水汽含量;θ为卫星观测角度;ε4、ε5分别为两个热红外通道4、5的地表比辐射率;A0、P、M都是分裂窗系数,A0是与大气状况有关的系数,P可以表达为卫星观测天顶角θ、水汽含量W以及地表比辐射率ε4、ε5的函数;M则可表示成ε4、ε5、VIS及水汽W的函数。a1~a13为常数,可以通过基于最小二乘法的多元一次线性拟合得到。
用新建的模型再次对地表温度进行拟合,拟合的系数和结果显示在表3中。
3.2 改进后模型评价
将原有的算法反演带来的误差与带入气溶胶这一参量之后的算法误差进行对比分析,分析结果如下表4所示。由于水平气象视距在7 km以上对分裂窗算法反演带来的误差影响可以忽略不计,所以表中仅统计了当气象水平视距7 km以内模型改进前后的两种通道组合在四种气溶胶模式下均方根误差,平均偏差。考虑到原先的模型中对1~7 km气象水平视距下的模型分别进行了统计,所以RMSE和BIAS取7个算法模型的平均值带入表格中,作为改进前的统计结果。结果发现,改进后的方程不但连续性增强,而且总体效果要优于分开模拟效果的平均水平。尤其是在城市气溶胶类型下新拟合方程的RMSE和BIAS有明显改善,两种通道组合的平均RMSE=0.415 K,而改进前RMSE=0.847 K。对乡村气溶胶和沙尘气溶胶改善不是很明显是因为原先的方程就有较高的精度。
注:阴影表示系数适用于能见度小于7 km的大气条件下,空白部分表示采用原来的系数。
4 结论
以FY-3 VIRR及MERSI资料为数据源,基于其热红外通道光谱特征,通过MODTRAN大气辐射传输模式,模拟中纬度冬季不同地表条件下的FY-3热红外通道(VIRR4、VIRR5和MERSI5)亮温数据,构建分裂窗地表温度反演模型。模拟时,在原有的模型中加入气溶胶光学厚度这一参量,探究气溶胶对构建算法的影响,并利用气溶胶改进分裂窗算法模型,最终建立包含气溶胶信息的分裂窗算法来反演地表温度。主要结论如下。
(1)基于FY-3A VIRR4、5通道和MERSI5通道亮温数据,在原有的分裂窗算法模型中加入气象水平视距这一参量,并分别拟合了无气溶胶、城市气溶胶、乡村气溶胶和沙尘气溶胶条件下算法的系数,使得改进的模型更加合理。改进后的模型拟合效果相比原模型均有改进。其中对城市气溶胶这一类型改善巨大,RMSE由0.814 K降到0.442 K,BIAS从0.998 K降到0.668 K。
(2)分析了不同气溶胶模式下热红外通道之间亮温差随地表温度升高的变化。结果显示,VIRR4、5通道亮温组合的亮温差值在不同气溶胶模式下受地表温度影响较大,而VIRR4和MERSI 5的组合在不同地表温度和气溶胶类型下均保持2 K左右的亮温差。
摘要：FY-3A上搭载的中分辨率成像光谱仪(MERSI)和扫描辐射计(VIRR)都包含热红外通道,可用于地表温度的反演。结合热红外通道的光谱响应函数,用MODTRAN模拟清洁大气、城市气溶胶、乡村气溶胶和沙尘气溶胶条件下的亮温值,通过分裂窗算法求解不同气溶胶类型和光学厚度下的地表温度。结果表明:气溶胶类型和光学厚度之间的差异会造成热红外通道的亮温差发生变化,这一现象在VIRR的热红外通道表现得尤为显著。通过拟合系数和气溶胶光学厚度之间的关系建立新的反演算法,发现:在分裂窗算法基础上加入气溶胶光学厚度这一参量之后,地表温度反演的均方根误差有所降低,尤其是在气溶胶浓度高时改进效果更佳,其中对城市气溶胶这一类型改善效果明显(RMSE由0.814 K下降到0.442 K,BIAS从0.998K降到0.668 K)。
反演算法 篇8
1 改进遗传算法
反演计算往往需要多次调用正分析计算, 多次的正反算结合, 通常会大大增加计算量及计算时间, 这就对进行反演计算时所采用的算法提出了更高的要求。遗传算法在许多优化问题中已得到了广泛应用, 但与遗传算法理论本身的应用潜力相比, 其推广应用还有很大空间[10]。遗传算法是全局优化自适应概率搜索算法, 它具有智能寻优、鲁棒性等优点, 但也存在早收敛、结果不精确等不足[11,12,13,14]。
一般遗传算法包含3个基本算子:繁殖 (reproduction) 、交叉 (crossover) 和变异 (mutation) [15], 其中, 交叉概率和变异概率是十分重要的计算参数, 它们的选择直接影响到算法的收敛性。当交叉概率变大时, 计算过程中新个体的产生就会加速, 但如果太大, 会较大程度地损坏遗传模式, 使得群体中较高适应度的新个体也不断发生变化。如果交叉概率很小, 计算搜索过程又将会变慢。当变异概率特别小时, 很难产生新一代的个体;当其特别大时, 遗传算法将会变成一般意义上的随机搜索算法。但交叉概率pc和变异概率pm的确定既重要又不易。如将自适应技术融入遗传算法, pc和pm就可以依据适应度自行变化。当个体适应值高于群体平均适应值时, 采用较小的pc和pm, 从而使个体被遗传到下一代的机会更大;当个体的适应值较低, 则采用较大的pc和pm, 这样较低适应值的个体将被舍弃。
交叉概率pc和变异概率pm自适应调整公式:
式中, fmax为群体中最好的适应值;favg为每代群体中的平均适应值;f'为个体适应值;pc1, pc2, pm1, pm2为给定值。
2 泵房结构有限元模型
对淮安三站泵房结构的动力参数反演建立在有限元正分析基础之上, 采用有限元软件ANSYS建立淮安三站泵房结构有限元动力计算模型 (图1、图2) , 为减小反演分析计算工作量, 采用弹簧考虑泵房周围土体作用的弹簧约束模型。上部厂房结构采用Beam4单元进行离散。块基结构包括底板、廊道、进出水流道、辅机层、上下游胸墙、工作桥以及厂房地坪等, 采用Solid45单元离散。
模型的坐标原点位于两机组中间, 废黄河零点的高程, 块基结构下游侧位置, 并采用笛卡尔坐标系统, 符合右手定则:x轴沿横河向, 以自右岸指向左岸为正;y轴沿竖直向, 以竖直向上为正;z轴沿顺河向, 以自上游指向下游为正。
3 泵房结构材料参数及边界条件反演
3.1 反演方法及反演程序流程
在淮安三站的实测资料中, 虽已测得泵房结构竖直向、顺河向、横河向各方向自振频率, 可联合各方向自振频率构造目标函数消除单方向基频目标函数计算精度不高、计算效果差的问题, 但由于AN-SYS软件计算系统频率时, 需要人为辨别系统各方向频率及扭转频率, 采用此反演方法程序编制不易实现。所以, 综合考虑第1阶频率与加速度峰值构造目标函数, 提出双重目标函数反演方法。
双重目标函数反演方法主要思想:第1阶段以实测的泵房结构第1阶频率构造第一重目标函数, 采用改进遗传算法, 反演出满足实测系统频率的多组待定解, 即为第一重目标函数最优解。第2阶段再对得到的多组第一重目标函数解, 以现场测试得到的系统加速度峰值构造第二重目标函数, 最终得出最优解。双重目标函数反演流程如图3所示。
目标函数的选择十分重要, 直接决定了反演计算的计算量以及反演结果的正确性。上述所提出的双重目标函数反演法, 在减小反演计算量的同时也保证了反演结果的唯一性。其中所采用的双重目标函数分别为:
相对于单目标函数反演方法, 即只考虑一种类型的实测数据构造一个目标函数的反演方法, 双重目标函数反演方法使用了2个目标函数, 联合考虑了结构的自振特性及动力响应, 使其具有独特的优点:在同时反演多个参数时, 单目标函数可能因为参数之间的互相影响而导致反演得到的最优解并非真实的最优解, 而双重目标函数则由于考虑了不同类型的实测数据, 将会使计算结果更为可靠。
3.2 反演结果
在采用双重目标函数反演法进行泵房结构动力材料参数及边界条件反演时, 根据实际情况, 在泵房顺河向块基 (高程-5.0~+10.8 m处) 施加弹簧单元, 建立用来模仿周围土体的、弹性约束作用的弹簧约束模型。
反演计算过程中将反复调用有限元动力计算, 计算量较大, 耗费时间较长, 故在群体规模上不宜取得过大, 本文种群大小M取30。反演变量上下限范围也不宜取得过广且应以包含最优解为界。本文泵房结构的动弹模上下限取为2.0×104~4.3×104MPa。弹簧模型中弹簧刚度系数参照《公路桥涵地基与基础设计规范》 (JTG D63—2007) 中桩基土弹簧计算公式K=ab1mz计算。其中, a为土层厚度;b1为桩基计算宽度;m为地基土的比例系数;z为土层中点距地面距离。经计算本文弹簧模型中弹簧刚度系数为6.1×107N/m。弹簧刚度系数取为3.1×107~7.9×107N/m。
经双重目标函数的改进遗传算法反分析计算, 得到反演参数分布 (图4—图6) 及反演计算过程 (图7) 。
从图4可以看出, 初始种群中反演参数的分布比较均匀, 这符合遗传算法随机产生第一代群体的规律, 但是最优个体附近个体较少。从图5可以看出, 进化了19代后, 最优个体附近已经出现了较多较好的个体了。
图7显示出了进化到第20代时得到了相对最优解, 且一直保持到第30代此解都没有变化, 说明本文所提的双重目标函数反演方法具有较高的搜索效率, 从目标函数值可看出该法具有较高的反演精度。经双重目标函数反演的最终的最优解为:
3.3 利用反演参数的泵房结构动力特性计算
利用弹簧模型反演的参数进行正分析, 并将得到的泵房结构自振频率与实测得到的自振频率值进行比较, 具体见表1。
从计算结果与实测值比较来看, 除第五阶和第七阶频率值外, 其他各阶误差均在10%以内, 说明反演得到的计算参数基本可靠, 基本满足计算精度要求, 计算所采用的弹簧模型能较好地反映实际情况。
为了更好地说明双重目标函数反演方法的计算效果, 本文采用仅考虑一阶频率的单目标函数进行了上述反演工作, 利用所得数据进行正分析计算, 将所得结果列于表2。
从表2中可以看出, 采用仅考虑一阶频率的单目标函数反演方法所得结果进行的正分析计算与实测值的比较中, 除第一阶频率的误差为4.92%, 其余各阶误差均超过20%。原因在于, 反演计算中的反演参数, 即动弹模与弹簧刚度系数, 它们均影响系统的刚度。在计算过程中, 较大动弹模值与较小弹簧刚度系数值的组合及较小动弹模值与较大弹簧刚度系数值的组合, 可能会使系统的计算第一阶频率值相等。因此, 在此情况下, 采用单目标函数反演方法所得结果并非一定为最优解。而双重目标函数反演方法由于采用了不同类型的实测数据进行控制, 很好地避免了上述问题。
反分析结果通常为正分析所用。对于很多实际工程, 因实际情况限制, 现场测试点较少, 从而难以反映结构的空间振动形态。因此, 现场测试结合数值计算结果来综合考虑显得尤为重要。
4 结语
反演算法 篇9
微动指来自于地表和底下深处的随机振动源引起的振幅非常小的地面微弱振动。包括人类等人工源引起的短周期震动和海浪及地球内部活动等自然力引起的长周期震动。
微动探测的主要方法有空间自相关法 (SAC) [1]和频率-波数法 (F-K) [2]两种。本文研究的是基于空间自相关法的微动勘探反演算法中遗传算法的应用于优化。
遗传算法 (Genetic Algorithms, 简称GA) [3]是一种基于生物自然选择与遗传机理的随机搜索与优化方法, 是由美国Michigan大学的John Holland教授创建的。它来源于达尔文的进化论、魏茨曼的物种选择学说和孟德尔的群体遗传学说, 其基本思想是模拟自然界遗传机制和生物进化论而形成的一种过程搜索最优解的算法。它本身不仅具有极强的鲁棒性, 还隐含着并行性, 同时还具有全局寻优能力。因此, 从一定意义上说, 遗传算法是一种普遍适用于各种问题的简单而有效的搜索方法。
1 遗传算法优化研究
1.1 遗传算法基本原理
基本遗传算法有五大要素构成:参数的编码方式、初始种群规模的设定、适应度函数的设计方法、各种遗传算子的设计和控制参数的设定, 他们组成了遗传算法的核心内容。[4]
(1) 编码:遗传算法不能直接处理问题空间的参数, 必须把它们转换成遗传空间的由基因按一定结构组成的染色体或个体。这一转换操作就叫做编码。
(2) 个体适应度:在遗传算法中, 搜索进化的过程是根据评估函数值来对每个个体进行优劣评估的, 这个过程不需要其他外部信息作为参考。这里的评估函数值就是适应度。
(3) 算子;遗传操作中算子主要包括:选择算子、交叉算子以及变异算子三种类型。
(4) 控制参数:遗传算法中有其运行的基本参数被称作操作参数, 在研究优化问题之前必须提前设定这些参数。主要的参数有4种:
N:种群规模, 即种群中个体的数量
T:算法停止条件, 即种群进化代数
Pc:交叉概率
Pm:变异概率
遗传算法实质上是一种以群体和遗传操作为基础进行繁衍、监测和评价的迭代算法。在搜索之前, 将问题的解映射为一个解集空间, 解集空间是由一定数量的二进制个体组成的种群, 以这些代表问题假设解的种群开始, 我们要解决的一个问题就是根据问题的目标函数构造一个适应度函数, 以适应度函数为依据, 对种群中的个体进行评估, 从一大堆数据中寻找一个解。图1为遗传算法的流程图。
1.2 遗传算法优化研究
凭借遗传算法具备的非线性全局优化的能力, 它已经被广发的应用到了诸多领域中, 在微动面波反演的研究中也引入了此方法。但是在解决反演问题的过程中也暴露了一些不足之处, 首先标准遗传算法的计算过程中有早熟收敛[5]的现象, 这一现象致使算法全局寻优的优良性能不能完全显示出来, 其次遗传算法的收敛速度相比于其他传统算法比较慢, 效率低。因此, 针对以上暴露的不足, 本文从四个方面入手, 提出了遗传算法的优化策略, 具体的策略实施过程如下:
(1) 自适应函数的优化
在遗传算法中目标函数通常会作为评估函数值即适应性函数, 在遗传进化中, 超级个体可能会在群体中出现, 就是适应值比群体平均值大很多的个体, 如果根据适应值的比例选择时, 这类超级个体就会占有群体中的绝对比例, 这样算法在运算时就会过早的收敛于一个局部最优值, 这就是所说的早熟现象。下面通过引进一种线性变换来使原有适应性函数得到改变, 避免早熟收敛的发生。这里设f为已有的适应值函数, F为改变的适应值函数用线性组合F=af+b表示。并设定
公式 (1) 中k为一个固定的增量系数。
由上式得:
原有适应度经过线性的变换之后他们之间的差距得到了调整, 群体内多样性的个体得到了保证, 而且这个方法简便计算, 容易实现。
(2) 选择策略的优化
利用最优保存和赌轮[6]选择两种方式相结合的策略作为选择的方法。在群体中首先找到适应值最佳和最差的个体, 把最差的个体用最佳个体替换, 然后允许最佳的两个个体不经过交叉编译而进入下一代, 剩下的个体按正常进行排列, 最后再进行赌轮选择, 这样群体中平均适应值不断得到提高, 而且还使最佳个体的适应值不再减小。
(3) 变异策略的优化
在标准遗传算法中, 群体中单一性的产生有一部分原因是近亲繁殖造成的, 因此为了避免这类事件的发生, 可以让父母A, B多次交叉变异, 在产生的不同的多个个体中筛选出最佳个体放入子代中, 对上述步骤反复执行, 直到子个体数量得到满足。
(4) 遗传操作算子的优化
遗传算法有收敛过早的缺点, 这一缺点是由于同位基因数目不同概率不等, 在进化中某些特定基因越来越少, 造成基因缺失, 致使收敛过早。针对这一现象, 本文引入多元算子到遗传算法中来。例如对二元字符串结构进行变异, 通过同或异或的运算保留了特定基因的多样性, 极大地避免了收敛过早的问题。例如:
经过同或异或后两种逻辑状态必定互补。
根据以上优化策略优化后的遗传算法执行顺序如下:
(1) 子群体的个数要首先确定, 对子群体初始化, 子群体每个的规模是N, 进一步对重新分配的子群体确定进化代数;
(2) 根据规则1, 把各子群体的父代适应值求解;
(3) 对各个子群进行选择操作;
(4) 获取各子群体繁殖的下一代, 变异的规则采用二元变异;
(5) 查看终止条件是否满足, 如过满足, 就退出;
(6) 如果进化代数满足了重新分配子群体的数量, 就开始对子群体重新分配;
(7) 回到 (2) 。
1.3 实验测试
遗传算法的随机优化计算的特点, 使其受参数的影响很大, 参数的不同会使算法对问题的计算产生不同的影响, 恰当的参数可以加速算法早熟收敛, 相反不合适的参数会致使遗传算法收敛变慢, 因此使得计算后无法寻得最优解。然而确定选取什么样的参数目前还没有有效的解决办法。本节从两个不同的方面对优化后的遗传算法进行了测试, 验证优化后的效果和其对参数的敏感度。实验所用函数为
遗传算法对比实验
将遗传算法同优化后的遗传算法进行对比。种群规模设为100, 交叉概率0.72, 变异概率0.15。实验结果如下表。
通过以上实验, 可以看出, 优化后的遗传算法增强了自身的寻优能力, 减少了迭代次数, 明显改善了早熟收敛的情况, 另外优化后的算法参数的选择范围扩大, 参数的选择对算法的影响减小了, 使得算法能够更方便的在实际中应用。
微动反演中遗传算法的实际应用与优化。
选择表2中的模型进行反演, 理论数据根据knoppoff算法的速度递增层模型正演计算得到理论频散曲线。
层状介质微动面波速度反演中反演参数是横波速度VS以及厚度d。纵波速度VP通过下面关于横波速度VS与泊松比μ关系的公式计算:
密度ρ根据纵波速度和密度的统计关系得到。因此实际反演中只对各层的横波速度和层厚度参数反演。因此, 输入参数为各层横波速度和层厚度的上下限。算法目标函数计算式为
式中Viobs是第i个频点上观测的微动面波速度, Vithe是每代在确定模型参数后, 计算得到的同一频点的微动面波速度值, N是频点的个数, 每次迭代时的每个个体是F, 就是模型的参数。
反演结果 (表2所示) 就是获得一组模型的参数, 此组模型参数让上式表述的目标数取值达到最小。所以, 这里就是求解目标函数的最小值。从表2中可以看到优化后的遗传算都各参数的反演结果的相对误差均优于标准遗传算法, 各层的相对误差均不超过10%。结合实际, 进行多次的试算, 具体采用以下参数:种群大小为40-60个, 遗传代数为160代;交叉概率为Pc1=0.92, Pc2=0.75, Pc3=0.6, 遗传概率Pm1=0.1, Pm2=0.09, Pm3=0.01。
由表2和图2所示可以看出, 优化的遗传算法反演结果的相对误差均低于标准遗传算法。从图3和4中可知, 两种方法一般在40-60代就得到最佳解, 其中优化的遗传算法反演1的收敛速度最快, 约在第25代即获得最佳解的。优化的算法运算速度比标准算法提升了2个数量级。
2小结
反演算法 篇10
近年来,在许多专家和学者的努力研究下,SNMR技术得到了进一步的完善,其中,在SNMR信号反演方面,一维正反演理论较为成熟,并正在深入研究二维、三维正反演模型[5—8]。SNMR信号反演都是以正演方程An=E为理论依据进行求解的,其中,A为核函数矩阵、E为各激发脉冲矩对应的SNMR信号初始振幅值、n为各层地质体含水量。文献[9]提出了改进的模拟退火算法(the improved simulated annealing algorithm,ISA)反演,提高了现有反演算法的稳定度和收敛速度;文献[10]提出了QT反演算法,利用各个激发脉冲矩对应的全部采样点数据进行反演,充分挖掘了接收信号信息,在一定程度上提高了反演精度;文献[11]采用了积分门技术接收信号,提高各个采样点数据的精度,并进行全衰减反演,是对QT反演的一种改进。探测区域背景电阻率空间分布情况决定核函数矩阵A的特性,而测点的电阻率值往往是通过一些物探方法或者周围区域现有的钻井信息获得,这使得核函数矩阵A必然存在一定误差。但是,当前文献刊登的反演算法都是在假设探测区域背景电阻率空间分布为精确值的条件下进行的,这影响了相关算法的反演精度,尤其是在背景区域电阻率值小于50欧米的条件下,背景电阻率存在较大估计误差时,会导致反演结果无效。在求解A和E都存在误差的矩阵方程时,总体最小二乘法(total least square,TLS)的性能优于最小二乘[12],为了提高SNMR信号反演结果的准确性,本文采用总体最小二乘反演模型;由于核函数矩阵的条件数很大,An=E是一个病态方程求解问题,而正则化方法是求解病态方程的有效手段[13],为了提高解的稳定性,本文构造了正则化-总体最小二乘(regularization-total least square,R-TLS)反演模型,并将该模型转化为受条件约束的非线性优化问题,同时,提出了改进的和声搜索算法(the improved harmony search algorithm,IHS)以求解该约束非线性优化反演问题。为了提高反演精度,往往将含水层的厚度设定为较小值,使得含水层数大于激发脉冲矩数,反演方程An=E成为欠定方程,在这种情况下,上述改进的和声搜索算法仍适用。
1 层状介质中非线性反演
1.1 激发磁场在层状介质中的空间分布
在层状介质模型中,通过Maxwell方程组可以推导出回线源周围点(r,z)处,第i层介质的磁场垂直分量Hiz和径向分量Hir满足表达式[5]:
式中,I0为激发电流强度;为拉莫尔频率,u0为真空磁导率;ρi为第i层电阻率值;J0()·为第一类0阶Bessel函数;J1()·为第一类1阶Bessel函数;ai和bi为与地质构造相关的参数。
1.2 地面核磁共振信号初始振幅模型
在水平层状介质中,核磁共振信号初始振幅离散化模型为[8]:
式(3)中,qi为第i个激发脉冲矩,nj和Δzj分别为第j层的含水量和厚度,K(qi,zj)为圆柱坐标系下离散化后的核函数。
式(4)中,M0为磁化强度;b1⊥=u0H1⊥/I0,H1⊥为激发磁场强度垂直于地磁场的分量;θimnj为微分单元(rm,Φn,zj)处的扳倒角;Δrm和ΔΦn分别为微分单位(rm,Φn)处的径向和横向长度;M和N分别为径向和横向微分单元的数目。
在地面核磁共振探测中,通过发射不同大小的激发脉冲矩qi,以激发不同深度处水中氢离子,将式(3)转化为矩阵方程:
式(5)中,
1.3 基于约束总体最小二乘法的SNMR反演模型
SNMR反演中,核函数精度依赖于各个地电层的电阻率值的选取,而电阻率值往往是通过垂直电法等物探方法探测得到,由于测量误差和各种噪声的干扰,使得电阻率值存在一定误差,进而影响到核函数的精度,实测SNMR信号初始振幅值同样存在误差。为了提高SNMR反演精度,建立了约束总体最小二乘模型,其目标函数为:
式(6)中,A*和E*为不含误差的核函数和SNMR信号初始振幅值;为F范数;δ为一个正常数;L为限制模型平坦度的偏导数矩阵,通常取1次或者2次偏导数。式(6)对应的约束Lagrange方程为:
式(7)中,μ为Lagrange因子。当δ趋于0时,式(7)转化为[13[14]14]:
式(8)中,ξ为约束R-TLS方程的正则化因子。
1.4 IHS算法求解含水量
和声搜索算法(harmony search algorithm,HS)[15]是2001年Geem Z.W.提出的一种模拟乐师即席创作的智能优化算法,具有不依赖初值、随机性和遍历性等优点,非常适合用于求解多维非线性优化问题。为了提高标准HS算法的求解性能,本文对记忆库选取概率(harmony memory considering rate,HMCR)、音调被调概率(pitch adjusting rate,PAR)、音调被调步长(bw)参数,依赖于迭代次数进行动态调整,使其前期具有良好的遍历性、后期求解具有高精度特性。为了求解式(8)目标函数中的含水量值,IHS算法实现步骤如下:
步骤1:算法初始化。随机产生数目为记忆库大小(harmony memory size,HMS)的解向量放入和声记忆库(harmony memory,HM)中,解向量长度与模型中含水层个数相等为MN,即
本文HMS=5,ni,j为0~1之间的随机数;将HM中各个解向量带入式(8),求其适应度值向量Fitness=(F1,F2,…,FHMS),Fitness中最大值为Fmax;最大迭代次数为Nmax=5 000;初始迭代次数为Num=0;反演精度阈值为Threshold=10-8。
步骤2:生成一个新的和声。HM*中每个元素ni*(i=1,2,…,MN)以概率HMCR从HM的中随机选取一个值,以概率1-HMCR选自[0,1]之间的随机数。若ni*取自HM,则以概率PAR对其进行bw步长的扰动,即ni*更新为ni*-bw或ni*+bw,由于含水量取值范围在0至1之间,所以ni*更新公式为
为了使HS算法前期具有良好的遍历性、后期求解具有高精度特性,本文对标准HS算法进行改进,使参数HMCR、PAR和bw根据迭代次数进行动态调整。即
步骤3:计算适应度函数,更新HM。将步骤2中新和声HM*带入公式(8)计算适应度函数值F*,如果F*<Fmax,则用HM*个更新HM中Fmax对应的解向量,并用F*更新Fmax。
步骤4:判断迭代是否停止。如果迭代次数达到Nmax,或者当前最优适应度函数值小于Threshold,则停止迭代,反演结果为Fitness最小值对应的解向量;否则,返回执行步骤2。
2 野外实测数据实验及结果分析
在野外实测数据反演实验中,反演数据出自老挝首都万象盆地(Vientiane basin,Laos)一次地面找水探测结果,选取测区1的Site1测点[3]。该次探测采用地面核磁共振和垂向电测深(vertical electrical sounding,VES)联合反演,并结合钻探取证,以探明地下水文地质构造。用反演得到的含水量值的方均根RMS评价反演结果性能,RMS公式为:
Site1测点区域地磁场强度为43 770 n T,使用1 864 Hz的拉莫尔频率,磁倾角为24°,铺设边长为100的方形线圈进行探测,实测数据信噪比为6.9d B,实测SNMR信号初始振幅值曲线如图1中的红色“-o-”所示。VES反演结果表明,测点位置地下垂深0~4 m区域电阻率值为550Ω·m,垂深4~11 m区域电阻率值为5 000Ω·m,垂深11~70 m区域电阻率值为20Ω·m。钻探结果显示,该区域垂深0~38 m为废弃河道形成的冲积岩,其主要成分是砾石、砂和黏土;垂深38~70 m为岩盐,其含水分布如图2中红色曲线所示。反演中,为了提高反演精度,含水层厚度设定为1 m。
从图1可以看出本文反演结果和实测数据拟合的较好;图2中蓝色曲线为本文算法反演结果,其方均根RMS为3.12%,绿色曲线为法国Samovar v6.2反演软件反演结果[3],其方均根RMS为3.65%,从图中可以看出,本文算法反演结果优于法国Samovar v6.2反演软件,与钻探结果吻合度更高,能更精确的反映地下水文地质构造。
3 结论
通过对层状地电模型中SNMR反演进行研究,发现了现有文献中的反演算法存在不足,即背景区域电阻率的分布存在估计误差,而现有文献中没能提出解决方法,仍采取忽略此估计误差的方法进行反演,降低了反演的准确度,尤其是在背景区域为低电阻率的条件下,由于估计误差的存在,经常会给出错误的反演结果。同时,考虑到反演核函数矩阵的条件数很大,为了在背景区域电阻率分布值序列和测量信号的初始振幅值序列都存在误差的条件下,提高反演精度,建立了约束R-TLS反演模型,通过公式推导将该模型转化为受条件约束的非线性优化问题,并提出了IHS算法以求解该约束非线性优化问题,IHS使得标准HS算法前期具有良好的遍历性、后期求解具有高精度特性。由于约束条件的限制,该IHS算法在反演矩阵方程为欠定方程时仍能精确求解,消除了反演模型中对含水层最大划分层数的限制。在对老挝首都万象盆地的实测数据进行反演验证中,本文算法的反演结果含水量值的方均根为3.12%,文献[3]中法国Samovar v6.2反演软件反演结果含水量值的方均根为3.65%,两种反演结果均与钻探结果接近,但本文算法略显优势。
摘要：地面核磁共振反演能抽象为一个求解矩阵方程An=E的问题,A为与背景电阻率空间分布有关的核函数矩阵,E为测量信号的初始振幅值,n为带求解的含水量分布值,由于A和E都存在误差,为了提高n的求解精度和稳定性,构造了正则化-总体最小二乘模型,并将该模型转化为受条件约束的非线性优化问题,设计了改进的和声搜索算法以求解该问题,在含水层数大于激发脉冲矩数的欠定方程或者病态方程的求解中,该算法仍然适用。野外实测数据反演中,导电层电阻率分布情况来自垂向电测深勘探结果,观测信号的信噪比为6.9 d B,算法的反演结果含水量值的方均根为3.12%,法国Samovar v6.2反演软件反演结果含水量值的方均根为3.65%,两种反演结果均与钻探结果接近,但本文算法略显优势。