时间域激电法25维正演算法的研究

时间域激电法25维正演算法的研究

0时间域激电法2.5维正则化反演的效率2.5维雍正算法相对成熟。时间间隔激发现象法数据的反演通常采用最小衰减二乘法，只选择最小调整误差，反演结果的分辨率较高。有些采用了正式的反演方案，但当实测数据量大、正反演网格剖面精细且反演波数选择较多时，伪反演法的计算效率较低，这将影响算法的实用性。因此，如何提高2.5维雍正反演的计算效率是一个需要进一步研究的问题。对于2.5维地电断面的正反演计算问题,需要对波数进行循环计算,当波数较多时计算量很大(罗延钟和孟永良,1986).对于这种大计算量问题,可分为两种方案节省时间,其一是在满足计算精度的同时尽可能减少波数个数来缩小计算量,从而实现快速正反演.此外也可通过提升计算速度来达到减少计算时间的目的.本文利用MPI并行方法进行处理,在满足精度的同时,加快计算速度,达到了减少计算时间的目的.在正则化反演中,由于每次拉格朗日乘子搜索都需要求解新模型、正演、计算拟合差等工作,因此当搜索次数较多时会减慢反演速度,本文对常规搜索方案修改优化,使得反演速度得以提升.本文对时间域激电法2.5维正反演算法的诸多细节进行了深入研究,实现了基于MPI并行计算的时间域激发极化法2.5维正则化反演.经理论模型算例测试,反演效果良好,达到了提高正反演计算效率的目的.12.5维冲突法1.1域的电位u时间域激发极化法2.5维正演是基于“等效电阻率”的概念(傅良魁,1991),在波数域采用有限单元法计算傅氏电位V,再通过反余弦变换得到空间域的电位U(罗延钟和张桂青,1987;麻昌英等,2014;刘云等,2014).视极化率数据根据两次电阻率法2.5维正演结果为式中η电阻率法正演采用二次场算法(鲁杏等,2014).将复杂的电阻率模型分解成背景电阻率模型和剩余电阻率模型.背景电阻率模型通常取均匀半空间或层状介质,其电位值可通过解析式或数字滤波法求得.复杂电阻率模型对应总电位V,背景电阻率模型对应一次电位V可以整理出二次电位所满足的方程(4)为当一次电位V1.2电阻率法正演精度.垂直接触带模型有解析解,因此可利用垂直接触带模型来检验有限元正演精度.垂直接触带左侧电阻率为10Ω·m,右侧电阻率为100Ω·m.图1为三极剖面法视电阻率有限元结果与解析解对比图.三极装置中AM=100m,MN=10m.可以看出有限元正演结果与解析解吻合较好,表明电阻率法正演精度是可靠的.激发极化法正演是基于“等效电阻率”概念通过两次电阻率法正演得到的,因而激发极化法正演精度也是可靠的.2极化率模型的反演对于时间域激发极化法数据的反演,本文采用分步反演的思路,先反演视电阻率数据获得电阻率模型,固定电阻率模型再反演视极化率数据,最后获得极化率模型.正则化反演是对模型参数的分布加入某种约束的反演方法.使用较多的方案是对模型的圆滑程度做约束,获得模型的圆滑解(Constableetal,1987;deGroot-Hedlin,1990;何梅兴等,2011).本文采用的正则化反演目标函数定义为式中U为目标函数,μ为拉格朗日乘子,m为模型,d为数据,X可推导出该正则化反演方法的模型更新迭代公式为式中J2.1傅氏电位模型视电阻率与电阻率参数对应的雅克比矩阵元素J式中ρ对傅氏电位所满足的线性方程组式(3)两边同时对地下某一单元的电阻率求偏导,可以得到下式为根据互换原理,在模型某一矩形单元节点a、b、c、d处分别供入强度为I2.2极化率反演算法视极化率与极化率参数对应的雅克比矩阵元素表达式为式中η因此,在视电阻率与电阻率参数对应的雅克比矩阵基础上,只需进行一些系数的处理计算,便可较快的获得视极化率与极化率参数对应的雅克比矩阵.32.5维正反演计算速度提高措施3.1mpi标准中多进程并行计算MPI是由全世界工业、科研和政府部门联合建立的一个消息传递编程标准,MPI标准中定义了一组函数接口用于进程间的消息传递,利用这些函数接口从而实现多进程的并行计算(黄易等,2010).3.1.1计算速度的设置在傅氏电位正演求取过程中,由于不同波数的傅氏电位计算相互独立,在每次正演时可以将不同波数的运算工作分配给不同进程,这样可以提升计算速度.根据波数个数和机器配置情况开设合理的进程个数(如N个).首先主进程(0)进行正演前的数据准备工作,进入并行阶段时,主进程(0)将所需数据传递给其他N-1个进程,每个进程计算分配波数的傅氏电位计算工作.当所有波数计算完后各次进程将计算结果传递回给主进程(0),主进程(0)再根据所有波数的傅氏电位计算结果进行反傅氏变换和后续的计算工作.3.1.2并行时的并行处理由于反演中波数域雅克比矩阵计算也存在一个波数循环,且各波数的计算工作相互独立,因此同正演一样,反演中波数域雅克比矩阵的计算也可以进行并行处理,达到加快计算速度的目的,具体处理过程与正演中傅氏电位并行大致相同,此处不多叙述.此外,对于并行计算,不同进程之间的数据传输效率也起到一定的影响,因此参考其他文献(李焱等,2010),将并行时常用的的数据传送和接收功能函数MPI_SEND、MPI_RECV改为聚合通信函数:MPI_GATHER、MPI_SCATTER和MPI_BCAST,因为这些函数具有更加高效的通信效率.3.2拉格朗日乘子的影响在正则化反演过程中,拉格朗日乘子起到了控制模型圆滑程度的作用,拉格朗日乘子越大,模型的圆滑权重越高(刘俊峰等,2013;吴小平和徐果明,1998).常规搜索方案多采用三点确定极小值区间搜索的方法,大致分为以下三个步骤:(1)确定拉格朗日乘子在前几次模型迭代过程中,数据拟合差通常较大,远高于期望拟合差要求,为了使其降低,需要不断地搜索拉格朗日乘子,找到某一拉格朗日乘子对应的模型数据拟合差低于本次迭代的初始数据拟合差,因此可寻找三个点建立一个极小值区间来确定该区域内存在使模型的拟合差降低的拉格朗日乘子点.(2)格朗日乘子拟合当极小值区间确定后,便需要在该区间内寻找到一个拉格朗日乘子,该点对应的模型数据拟合差为该区间内所有点拟合差的极小值,若该极小值仍高于期望拟合差,则保留该点对应的模型作为初始模型,进入下次迭代.(3)拉格朗日乘子搜索方案优化当某次迭代的拟合差极小值小于期望拟合差时,由于过度的数据拟合同样会造成一些不必要的假异常,因此便不再追求更低的拟合差,而是寻找满足期望拟合差且拉格朗日乘子最大的点所对应的模型,此结果即为我们所要寻找的圆滑模型.在搜索尝试不同的拉格朗日乘子时,首先需要进行大型矩阵求逆工作来求解出该点所对应的模型,再对该模型进行正演来计算数据拟合差,由于常规拉格朗日乘子选取方案需要对其进行多次搜索,因此当搜索次数过多时会增大计算量,若能去除一些不必要的搜索步骤,便能够起到加快反演速度的作用.常规搜索方案首先要确定极小值区间、然后确定区间极小值点、再确定极小值是否小于期望拟合差,最后进入搜索交叉点步骤,如图2中b点所示,所以按照规定步骤要求搜索时便会浪费一些不必要的正演时间.在常规步骤(1)确定极小值区间时,即使极小值区间尚未确定,但在搜索过程中只要找到一个点对应的拟合差小于期望拟合差,便无需再继续进行确定极小值区间和确定极小值的步骤,直接跳至常规搜索方案步骤(3)确定交叉点即可,这样相比原来可节省一定的搜索正演次数.如图3a中a点所示,此时三点尚未形成极小区间,但此时便可跳出步骤(1),进入步骤(3).同理,对于常规搜索方案步骤(2)搜索极小值,当某一点的拟合差小于期望拟合差时,无论该点是否为该区间所对应的极小值,便可以跳出该步骤,进入常规搜索方案步骤(3)搜索交叉点的阶段,这样也可以节省一些正演次数,达到节省时间的目的.如图3b中b点所示,此时b点还不是该区间的极小值点,原方案仍会继续搜索,但此时可以跳出步骤(2),直接进入步骤(3)阶段.4计算测试4.1.1极化率反演理论本算例设计的理论地电模型为一个低阻高极化棱柱体.棱柱体上顶面埋深40m、宽60m、高30m,异常体电阻率为10Ω·m,极化率12%,围岩的电阻率100Ω·m,极化率2%,具体位置如图4中黑虚线所示.利用该模型正演获得偶极-偶极装置数据,o-o′选取了20m、30m…690m共68个极距,每个极距观测数据量逐渐递减,采用了类高密度电法的采集方式,总数据个数为2346个.正演数据中加入1.5%的高斯随机误差,形成理论模型合成数据,电阻率和极化率反演的期望拟合差分别设定为1和1.2,然后对该数据进行反演,获得的电阻率反演结果如图4所示.电阻率反演的初始拟合差为15.16,经8次迭代,最终拟合差降为0.99,反演结束,从上图电阻率的反演结果来看,反演结果模型中低阻异常体的分布范围与真实边界吻合较好,异常体中心部位的电阻率也较为接近真实模型的电阻率,合成数据中存在的数据误差对反演结果并没有造成较大的影响,由于正则化约束的影响,异常体与背景围岩之间的电阻率过渡自然.利用该电阻率模型进行视极化率数据反演,获得的极化率反演结果如图5所示.极化率反演的初始拟合差为32.04经8次迭代,最终拟合差降为1.19,反演结束,极化率的反演结果中高极化率异常体的几何形态与位置与真实模型比较接近,反演结果较为圆滑,异常体与周围背景围岩极化率之间同样存在过渡的过程.4.1.2反演算例并行效率分析模型算例的网格大小为132×48,选取了18个波数.测试机配置为,Intel(R)Core(TM)i7-950CPU(3.07GHz),8G内存,在Linux系统下开设了3个进程进行计算.下面对并行计算所花的时间与非并行计算进行对比分析.正演计算对比:由于在一次反演过程中需要进行多次正演计算工作,因此即使一次正演花费的时间可以接受,但多次累积下来也会拖慢反演的速度,所以加快正演的计算速度是很有必要的.对模型算例中未加入并行处理的程序进行多次正演运算计时,平均计算时间为35s,而加入并行处理的程序每个进程计算6个波数的计算工作,平均正演时间为14s,可以看出并行处理能够节省一半多的正演时间.反演时波数域雅克比矩阵计算对比:由于波数域雅克比矩阵的计算工作中也存在3重循环,其计算量比正演的计算量大很多,因此该部分的计算速度对反演整体速度影响比重较大.对未加入并行处理的程序进行波数域雅克比矩阵运算计时,平均计算时间约为921s,近15分钟,而加入并行处理的程序也同样每个进程计算6个波数的计算工作,平均计算的时间为314s,近5分钟,节省了约三分之二左右的计算时间,因此对波数域雅克比矩阵计算加入并行处理的效果是显而易见的.由以上结果可以看出,对波数域雅克比矩阵和傅氏电位的计算进行并行处理可大幅度提升计算效率.4.1.3拉格朗日乘子搜索方案改进测试在对拉格朗日乘子的搜索过程中,每次找到一个新的拉格朗日乘子点,都需要进行求解新模型、模型正演、计算新模型数据拟合差等工作,因此省略一次搜索步骤不仅仅只是省略一次正演的工作,可以节省很多计算量.视极化率数据的反演加入了正则化约束,对该反演过程中拉格朗日乘子搜索次数进行计数,常规流程中拉格朗日乘子总搜索次数为61次,其中开始寻找交叉点前的搜索次数为55次.对改进方案的拉格朗日乘子搜索次数进行计数,其中总搜索次数为56次,开始寻找交叉点前的搜索次数为50次.因此,极化率反演改进方案可以减少5次(约占总搜索次数的8.2%)的搜索计算工作.视极化率数据的反演工作建立在视电阻率反演模型之上,对视电阻率数据的反演同样加入了正则化约束,对该反演中的拉格朗日乘子搜索次数进行计数,常规流程的搜索总搜索次数为121次,其中开始寻找交叉点前的搜索次数为114次.同样,对改进方案进行计数,总搜索次数为104次,而寻找交叉点前的搜索次数为93次.可以看出,此电阻率反演过程中可节省17次(约占总搜索次数的14%)的搜索计算工作.4.2.1起伏地表模型反演此算例的模型设计为地形起伏条件下,存在一个棱柱异常体,该异常体是电阻率为10Ω·m、极化率12%的低阻高极化异常体,上顶面距离地表30m、宽40m、高15m,具体位置如图6中黑虚线所示.地形为高16m宽85m的山峰.围岩的电阻率为100Ω·m,极化率为2%.同上一个算例相似,利用该模型正演获得偶极-偶极装置数据,由于模型中包含了起伏地形,需要网格剖分更加精细,因此采集数据量有所减少,o-o’选取了30m、50m…770m共38个极距,采集的总数据点个数为741.正演数据同样加入1.5%的高斯随机误差,形成理论模型合成数据,电阻率和极化率反演的期望拟合差都设定为1,然后对该数据进行反演,获得的电阻率反演结果如图6所示.电阻率反演的初始拟合差为7.72,经7次迭代,最终拟合差降为1.00,反演结束,从图中可以看出,反演结果相对准确,低阻异常体的形态和位置也能被清晰的刻画出来,电阻率也较为接近真实模型的电阻率,但地形附近的网格单元电阻率出现了一些假异常.利用该电阻率模型进行视极化率数据反演,获得的极化率反演结果如图7所示.极化率反演的初始拟合差为18.8,经7次迭代,最终拟合差降为0.99,反演结束,在此算例反演结果中,可以将高极化率异常体从背景围岩中准确定位,异常体的几何形态也与真实模型较为近似,反演结果较为圆滑,因此可以说明本文的反演程序在应对地形起伏的影响时,可以较好的排除干扰,准确的反演出真实的极化率异常体位置.4.2.2算例并行效率分析由于模型中包含起伏地形,因此地表附近的网格单元需要精细剖分,此算例的网格大小为191×61,共选取了15个波数.在Linux系统下开设了4个进程进行计算,其中3个进程计算4个波数的计算量,而另一个进程计算3个波数的计算量.下面对并行计算所花的时间与非并行计算进行对比分析.正演计算对比:对模型算例中未加入并行处理的程序进行多次正演运算计时,平均计算时间为107s,而加入并行处理的程序,平均正演时间为35s,可以看出并行处理能够节省三分之二的正演时间.反演时波数域雅克比矩阵计算对比:对模型算例中未加入并行处理的程序进行波数域雅克比矩阵运算计时,平均计算时间约为641s,近10min,而加入并行处理的程序,平均计算的时间为173s,约3min,节省了近三分之二多的时间,因此本算例对波数域雅克比矩阵计算加入并行处理的效果同样很突出.4.2.3拉格朗日乘子搜索方案改进测试同上例一样,对极化率反演过程中拉格朗日乘子搜索次数进行计数,常规流程中总搜索次数为55次,其中进行寻找交叉点前的搜索次数为49次.对改进方案的拉格朗日乘子搜索次数进行计数,其中总搜索次数为50次,开始寻找交叉点前的搜索次数为44次.对于此算例的极化