時(shí)間域激電法25維正演算法的研究

時(shí)間域激電法25維正演算法的研究

0時(shí)間域激電法2.5維正則化反演的效率2.5維雍正算法相對(duì)成熟。時(shí)間間隔激發(fā)現(xiàn)象法數(shù)據(jù)的反演通常采用最小衰減二乘法，只選擇最小調(diào)整誤差，反演結(jié)果的分辨率較高。有些采用了正式的反演方案，但當(dāng)實(shí)測(cè)數(shù)據(jù)量大、正反演網(wǎng)格剖面精細(xì)且反演波數(shù)選擇較多時(shí)，偽反演法的計(jì)算效率較低，這將影響算法的實(shí)用性。因此，如何提高2.5維雍正反演的計(jì)算效率是一個(gè)需要進(jìn)一步研究的問(wèn)題。對(duì)于2.5維地電斷面的正反演計(jì)算問(wèn)題,需要對(duì)波數(shù)進(jìn)行循環(huán)計(jì)算,當(dāng)波數(shù)較多時(shí)計(jì)算量很大(羅延鐘和孟永良,1986).對(duì)于這種大計(jì)算量問(wèn)題,可分為兩種方案節(jié)省時(shí)間,其一是在滿足計(jì)算精度的同時(shí)盡可能減少波數(shù)個(gè)數(shù)來(lái)縮小計(jì)算量,從而實(shí)現(xiàn)快速正反演.此外也可通過(guò)提升計(jì)算速度來(lái)達(dá)到減少計(jì)算時(shí)間的目的.本文利用MPI并行方法進(jìn)行處理,在滿足精度的同時(shí),加快計(jì)算速度,達(dá)到了減少計(jì)算時(shí)間的目的.在正則化反演中,由于每次拉格朗日乘子搜索都需要求解新模型、正演、計(jì)算擬合差等工作,因此當(dāng)搜索次數(shù)較多時(shí)會(huì)減慢反演速度,本文對(duì)常規(guī)搜索方案修改優(yōu)化,使得反演速度得以提升.本文對(duì)時(shí)間域激電法2.5維正反演算法的諸多細(xì)節(jié)進(jìn)行了深入研究,實(shí)現(xiàn)了基于MPI并行計(jì)算的時(shí)間域激發(fā)極化法2.5維正則化反演.經(jīng)理論模型算例測(cè)試,反演效果良好,達(dá)到了提高正反演計(jì)算效率的目的.12.5維沖突法1.1域的電位u時(shí)間域激發(fā)極化法2.5維正演是基于“等效電阻率”的概念(傅良魁,1991),在波數(shù)域采用有限單元法計(jì)算傅氏電位V,再通過(guò)反余弦變換得到空間域的電位U(羅延鐘和張桂青,1987;麻昌英等,2014;劉云等,2014).視極化率數(shù)據(jù)根據(jù)兩次電阻率法2.5維正演結(jié)果為式中η電阻率法正演采用二次場(chǎng)算法(魯杏等,2014).將復(fù)雜的電阻率模型分解成背景電阻率模型和剩余電阻率模型.背景電阻率模型通常取均勻半空間或?qū)訝罱橘|(zhì),其電位值可通過(guò)解析式或數(shù)字濾波法求得.復(fù)雜電阻率模型對(duì)應(yīng)總電位V,背景電阻率模型對(duì)應(yīng)一次電位V可以整理出二次電位所滿足的方程(4)為當(dāng)一次電位V1.2電阻率法正演精度.垂直接觸帶模型有解析解,因此可利用垂直接觸帶模型來(lái)檢驗(yàn)有限元正演精度.垂直接觸帶左側(cè)電阻率為10Ω·m,右側(cè)電阻率為100Ω·m.圖1為三極剖面法視電阻率有限元結(jié)果與解析解對(duì)比圖.三極裝置中AM=100m,MN=10m.可以看出有限元正演結(jié)果與解析解吻合較好,表明電阻率法正演精度是可靠的.激發(fā)極化法正演是基于“等效電阻率”概念通過(guò)兩次電阻率法正演得到的,因而激發(fā)極化法正演精度也是可靠的.2極化率模型的反演對(duì)于時(shí)間域激發(fā)極化法數(shù)據(jù)的反演,本文采用分步反演的思路,先反演視電阻率數(shù)據(jù)獲得電阻率模型,固定電阻率模型再反演視極化率數(shù)據(jù),最后獲得極化率模型.正則化反演是對(duì)模型參數(shù)的分布加入某種約束的反演方法.使用較多的方案是對(duì)模型的圓滑程度做約束,獲得模型的圓滑解(Constableetal,1987;deGroot-Hedlin,1990;何梅興等,2011).本文采用的正則化反演目標(biāo)函數(shù)定義為式中U為目標(biāo)函數(shù),μ為拉格朗日乘子,m為模型,d為數(shù)據(jù),X可推導(dǎo)出該正則化反演方法的模型更新迭代公式為式中J2.1傅氏電位模型視電阻率與電阻率參數(shù)對(duì)應(yīng)的雅克比矩陣元素J式中ρ對(duì)傅氏電位所滿足的線性方程組式(3)兩邊同時(shí)對(duì)地下某一單元的電阻率求偏導(dǎo),可以得到下式為根據(jù)互換原理,在模型某一矩形單元節(jié)點(diǎn)a、b、c、d處分別供入強(qiáng)度為I2.2極化率反演算法視極化率與極化率參數(shù)對(duì)應(yīng)的雅克比矩陣元素表達(dá)式為式中η因此,在視電阻率與電阻率參數(shù)對(duì)應(yīng)的雅克比矩陣基礎(chǔ)上,只需進(jìn)行一些系數(shù)的處理計(jì)算,便可較快的獲得視極化率與極化率參數(shù)對(duì)應(yīng)的雅克比矩陣.32.5維正反演計(jì)算速度提高措施3.1mpi標(biāo)準(zhǔn)中多進(jìn)程并行計(jì)算MPI是由全世界工業(yè)、科研和政府部門(mén)聯(lián)合建立的一個(gè)消息傳遞編程標(biāo)準(zhǔn),MPI標(biāo)準(zhǔn)中定義了一組函數(shù)接口用于進(jìn)程間的消息傳遞,利用這些函數(shù)接口從而實(shí)現(xiàn)多進(jìn)程的并行計(jì)算(黃易等,2010).3.1.1計(jì)算速度的設(shè)置在傅氏電位正演求取過(guò)程中,由于不同波數(shù)的傅氏電位計(jì)算相互獨(dú)立,在每次正演時(shí)可以將不同波數(shù)的運(yùn)算工作分配給不同進(jìn)程,這樣可以提升計(jì)算速度.根據(jù)波數(shù)個(gè)數(shù)和機(jī)器配置情況開(kāi)設(shè)合理的進(jìn)程個(gè)數(shù)(如N個(gè)).首先主進(jìn)程(0)進(jìn)行正演前的數(shù)據(jù)準(zhǔn)備工作,進(jìn)入并行階段時(shí),主進(jìn)程(0)將所需數(shù)據(jù)傳遞給其他N-1個(gè)進(jìn)程,每個(gè)進(jìn)程計(jì)算分配波數(shù)的傅氏電位計(jì)算工作.當(dāng)所有波數(shù)計(jì)算完后各次進(jìn)程將計(jì)算結(jié)果傳遞回給主進(jìn)程(0),主進(jìn)程(0)再根據(jù)所有波數(shù)的傅氏電位計(jì)算結(jié)果進(jìn)行反傅氏變換和后續(xù)的計(jì)算工作.3.1.2并行時(shí)的并行處理由于反演中波數(shù)域雅克比矩陣計(jì)算也存在一個(gè)波數(shù)循環(huán),且各波數(shù)的計(jì)算工作相互獨(dú)立,因此同正演一樣,反演中波數(shù)域雅克比矩陣的計(jì)算也可以進(jìn)行并行處理,達(dá)到加快計(jì)算速度的目的,具體處理過(guò)程與正演中傅氏電位并行大致相同,此處不多敘述.此外,對(duì)于并行計(jì)算,不同進(jìn)程之間的數(shù)據(jù)傳輸效率也起到一定的影響,因此參考其他文獻(xiàn)(李焱等,2010),將并行時(shí)常用的的數(shù)據(jù)傳送和接收功能函數(shù)MPI_SEND、MPI_RECV改為聚合通信函數(shù):MPI_GATHER、MPI_SCATTER和MPI_BCAST,因?yàn)檫@些函數(shù)具有更加高效的通信效率.3.2拉格朗日乘子的影響在正則化反演過(guò)程中,拉格朗日乘子起到了控制模型圓滑程度的作用,拉格朗日乘子越大,模型的圓滑權(quán)重越高(劉俊峰等,2013;吳小平和徐果明,1998).常規(guī)搜索方案多采用三點(diǎn)確定極小值區(qū)間搜索的方法,大致分為以下三個(gè)步驟:(1)確定拉格朗日乘子在前幾次模型迭代過(guò)程中,數(shù)據(jù)擬合差通常較大,遠(yuǎn)高于期望擬合差要求,為了使其降低,需要不斷地搜索拉格朗日乘子,找到某一拉格朗日乘子對(duì)應(yīng)的模型數(shù)據(jù)擬合差低于本次迭代的初始數(shù)據(jù)擬合差,因此可尋找三個(gè)點(diǎn)建立一個(gè)極小值區(qū)間來(lái)確定該區(qū)域內(nèi)存在使模型的擬合差降低的拉格朗日乘子點(diǎn).(2)格朗日乘子擬合當(dāng)極小值區(qū)間確定后,便需要在該區(qū)間內(nèi)尋找到一個(gè)拉格朗日乘子,該點(diǎn)對(duì)應(yīng)的模型數(shù)據(jù)擬合差為該區(qū)間內(nèi)所有點(diǎn)擬合差的極小值,若該極小值仍高于期望擬合差,則保留該點(diǎn)對(duì)應(yīng)的模型作為初始模型,進(jìn)入下次迭代.(3)拉格朗日乘子搜索方案優(yōu)化當(dāng)某次迭代的擬合差極小值小于期望擬合差時(shí),由于過(guò)度的數(shù)據(jù)擬合同樣會(huì)造成一些不必要的假異常,因此便不再追求更低的擬合差,而是尋找滿足期望擬合差且拉格朗日乘子最大的點(diǎn)所對(duì)應(yīng)的模型,此結(jié)果即為我們所要尋找的圓滑模型.在搜索嘗試不同的拉格朗日乘子時(shí),首先需要進(jìn)行大型矩陣求逆工作來(lái)求解出該點(diǎn)所對(duì)應(yīng)的模型,再對(duì)該模型進(jìn)行正演來(lái)計(jì)算數(shù)據(jù)擬合差,由于常規(guī)拉格朗日乘子選取方案需要對(duì)其進(jìn)行多次搜索,因此當(dāng)搜索次數(shù)過(guò)多時(shí)會(huì)增大計(jì)算量,若能去除一些不必要的搜索步驟,便能夠起到加快反演速度的作用.常規(guī)搜索方案首先要確定極小值區(qū)間、然后確定區(qū)間極小值點(diǎn)、再確定極小值是否小于期望擬合差,最后進(jìn)入搜索交叉點(diǎn)步驟,如圖2中b點(diǎn)所示,所以按照規(guī)定步驟要求搜索時(shí)便會(huì)浪費(fèi)一些不必要的正演時(shí)間.在常規(guī)步驟(1)確定極小值區(qū)間時(shí),即使極小值區(qū)間尚未確定,但在搜索過(guò)程中只要找到一個(gè)點(diǎn)對(duì)應(yīng)的擬合差小于期望擬合差,便無(wú)需再繼續(xù)進(jìn)行確定極小值區(qū)間和確定極小值的步驟,直接跳至常規(guī)搜索方案步驟(3)確定交叉點(diǎn)即可,這樣相比原來(lái)可節(jié)省一定的搜索正演次數(shù).如圖3a中a點(diǎn)所示,此時(shí)三點(diǎn)尚未形成極小區(qū)間,但此時(shí)便可跳出步驟(1),進(jìn)入步驟(3).同理,對(duì)于常規(guī)搜索方案步驟(2)搜索極小值,當(dāng)某一點(diǎn)的擬合差小于期望擬合差時(shí),無(wú)論該點(diǎn)是否為該區(qū)間所對(duì)應(yīng)的極小值,便可以跳出該步驟,進(jìn)入常規(guī)搜索方案步驟(3)搜索交叉點(diǎn)的階段,這樣也可以節(jié)省一些正演次數(shù),達(dá)到節(jié)省時(shí)間的目的.如圖3b中b點(diǎn)所示,此時(shí)b點(diǎn)還不是該區(qū)間的極小值點(diǎn),原方案仍會(huì)繼續(xù)搜索,但此時(shí)可以跳出步驟(2),直接進(jìn)入步驟(3)階段.4計(jì)算測(cè)試4.1.1極化率反演理論本算例設(shè)計(jì)的理論地電模型為一個(gè)低阻高極化棱柱體.棱柱體上頂面埋深40m、寬60m、高30m,異常體電阻率為10Ω·m,極化率12%,圍巖的電阻率100Ω·m,極化率2%,具體位置如圖4中黑虛線所示.利用該模型正演獲得偶極-偶極裝置數(shù)據(jù),o-o′選取了20m、30m…690m共68個(gè)極距,每個(gè)極距觀測(cè)數(shù)據(jù)量逐漸遞減,采用了類(lèi)高密度電法的采集方式,總數(shù)據(jù)個(gè)數(shù)為2346個(gè).正演數(shù)據(jù)中加入1.5%的高斯隨機(jī)誤差,形成理論模型合成數(shù)據(jù),電阻率和極化率反演的期望擬合差分別設(shè)定為1和1.2,然后對(duì)該數(shù)據(jù)進(jìn)行反演,獲得的電阻率反演結(jié)果如圖4所示.電阻率反演的初始擬合差為15.16,經(jīng)8次迭代,最終擬合差降為0.99,反演結(jié)束,從上圖電阻率的反演結(jié)果來(lái)看,反演結(jié)果模型中低阻異常體的分布范圍與真實(shí)邊界吻合較好,異常體中心部位的電阻率也較為接近真實(shí)模型的電阻率,合成數(shù)據(jù)中存在的數(shù)據(jù)誤差對(duì)反演結(jié)果并沒(méi)有造成較大的影響,由于正則化約束的影響,異常體與背景圍巖之間的電阻率過(guò)渡自然.利用該電阻率模型進(jìn)行視極化率數(shù)據(jù)反演,獲得的極化率反演結(jié)果如圖5所示.極化率反演的初始擬合差為32.04經(jīng)8次迭代,最終擬合差降為1.19,反演結(jié)束,極化率的反演結(jié)果中高極化率異常體的幾何形態(tài)與位置與真實(shí)模型比較接近,反演結(jié)果較為圓滑,異常體與周?chē)尘皣鷰r極化率之間同樣存在過(guò)渡的過(guò)程.4.1.2反演算例并行效率分析模型算例的網(wǎng)格大小為132×48,選取了18個(gè)波數(shù).測(cè)試機(jī)配置為,Intel(R)Core(TM)i7-950CPU(3.07GHz),8G內(nèi)存,在Linux系統(tǒng)下開(kāi)設(shè)了3個(gè)進(jìn)程進(jìn)行計(jì)算.下面對(duì)并行計(jì)算所花的時(shí)間與非并行計(jì)算進(jìn)行對(duì)比分析.正演計(jì)算對(duì)比:由于在一次反演過(guò)程中需要進(jìn)行多次正演計(jì)算工作,因此即使一次正演花費(fèi)的時(shí)間可以接受,但多次累積下來(lái)也會(huì)拖慢反演的速度,所以加快正演的計(jì)算速度是很有必要的.對(duì)模型算例中未加入并行處理的程序進(jìn)行多次正演運(yùn)算計(jì)時(shí),平均計(jì)算時(shí)間為35s,而加入并行處理的程序每個(gè)進(jìn)程計(jì)算6個(gè)波數(shù)的計(jì)算工作,平均正演時(shí)間為14s,可以看出并行處理能夠節(jié)省一半多的正演時(shí)間.反演時(shí)波數(shù)域雅克比矩陣計(jì)算對(duì)比:由于波數(shù)域雅克比矩陣的計(jì)算工作中也存在3重循環(huán),其計(jì)算量比正演的計(jì)算量大很多,因此該部分的計(jì)算速度對(duì)反演整體速度影響比重較大.對(duì)未加入并行處理的程序進(jìn)行波數(shù)域雅克比矩陣運(yùn)算計(jì)時(shí),平均計(jì)算時(shí)間約為921s,近15分鐘,而加入并行處理的程序也同樣每個(gè)進(jìn)程計(jì)算6個(gè)波數(shù)的計(jì)算工作,平均計(jì)算的時(shí)間為314s,近5分鐘,節(jié)省了約三分之二左右的計(jì)算時(shí)間,因此對(duì)波數(shù)域雅克比矩陣計(jì)算加入并行處理的效果是顯而易見(jiàn)的.由以上結(jié)果可以看出,對(duì)波數(shù)域雅克比矩陣和傅氏電位的計(jì)算進(jìn)行并行處理可大幅度提升計(jì)算效率.4.1.3拉格朗日乘子搜索方案改進(jìn)測(cè)試在對(duì)拉格朗日乘子的搜索過(guò)程中,每次找到一個(gè)新的拉格朗日乘子點(diǎn),都需要進(jìn)行求解新模型、模型正演、計(jì)算新模型數(shù)據(jù)擬合差等工作,因此省略一次搜索步驟不僅僅只是省略一次正演的工作,可以節(jié)省很多計(jì)算量.視極化率數(shù)據(jù)的反演加入了正則化約束,對(duì)該反演過(guò)程中拉格朗日乘子搜索次數(shù)進(jìn)行計(jì)數(shù),常規(guī)流程中拉格朗日乘子總搜索次數(shù)為61次,其中開(kāi)始尋找交叉點(diǎn)前的搜索次數(shù)為55次.對(duì)改進(jìn)方案的拉格朗日乘子搜索次數(shù)進(jìn)行計(jì)數(shù),其中總搜索次數(shù)為56次,開(kāi)始尋找交叉點(diǎn)前的搜索次數(shù)為50次.因此,極化率反演改進(jìn)方案可以減少5次(約占總搜索次數(shù)的8.2%)的搜索計(jì)算工作.視極化率數(shù)據(jù)的反演工作建立在視電阻率反演模型之上,對(duì)視電阻率數(shù)據(jù)的反演同樣加入了正則化約束,對(duì)該反演中的拉格朗日乘子搜索次數(shù)進(jìn)行計(jì)數(shù),常規(guī)流程的搜索總搜索次數(shù)為121次,其中開(kāi)始尋找交叉點(diǎn)前的搜索次數(shù)為114次.同樣,對(duì)改進(jìn)方案進(jìn)行計(jì)數(shù),總搜索次數(shù)為104次,而尋找交叉點(diǎn)前的搜索次數(shù)為93次.可以看出,此電阻率反演過(guò)程中可節(jié)省17次(約占總搜索次數(shù)的14%)的搜索計(jì)算工作.4.2.1起伏地表模型反演此算例的模型設(shè)計(jì)為地形起伏條件下,存在一個(gè)棱柱異常體,該異常體是電阻率為10Ω·m、極化率12%的低阻高極化異常體,上頂面距離地表30m、寬40m、高15m,具體位置如圖6中黑虛線所示.地形為高16m寬85m的山峰.圍巖的電阻率為100Ω·m,極化率為2%.同上一個(gè)算例相似,利用該模型正演獲得偶極-偶極裝置數(shù)據(jù),由于模型中包含了起伏地形,需要網(wǎng)格剖分更加精細(xì),因此采集數(shù)據(jù)量有所減少,o-o’選取了30m、50m…770m共38個(gè)極距,采集的總數(shù)據(jù)點(diǎn)個(gè)數(shù)為741.正演數(shù)據(jù)同樣加入1.5%的高斯隨機(jī)誤差,形成理論模型合成數(shù)據(jù),電阻率和極化率反演的期望擬合差都設(shè)定為1,然后對(duì)該數(shù)據(jù)進(jìn)行反演,獲得的電阻率反演結(jié)果如圖6所示.電阻率反演的初始擬合差為7.72,經(jīng)7次迭代,最終擬合差降為1.00,反演結(jié)束,從圖中可以看出,反演結(jié)果相對(duì)準(zhǔn)確,低阻異常體的形態(tài)和位置也能被清晰的刻畫(huà)出來(lái),電阻率也較為接近真實(shí)模型的電阻率,但地形附近的網(wǎng)格單元電阻率出現(xiàn)了一些假異常.利用該電阻率模型進(jìn)行視極化率數(shù)據(jù)反演,獲得的極化率反演結(jié)果如圖7所示.極化率反演的初始擬合差為18.8,經(jīng)7次迭代,最終擬合差降為0.99,反演結(jié)束,在此算例反演結(jié)果中,可以將高極化率異常體從背景圍巖中準(zhǔn)確定位,異常體的幾何形態(tài)也與真實(shí)模型較為近似,反演結(jié)果較為圓滑,因此可以說(shuō)明本文的反演程序在應(yīng)對(duì)地形起伏的影響時(shí),可以較好的排除干擾,準(zhǔn)確的反演出真實(shí)的極化率異常體位置.4.2.2算例并行效率分析由于模型中包含起伏地形,因此地表附近的網(wǎng)格單元需要精細(xì)剖分,此算例的網(wǎng)格大小為191×61,共選取了15個(gè)波數(shù).在Linux系統(tǒng)下開(kāi)設(shè)了4個(gè)進(jìn)程進(jìn)行計(jì)算,其中3個(gè)進(jìn)程計(jì)算4個(gè)波數(shù)的計(jì)算量,而另一個(gè)進(jìn)程計(jì)算3個(gè)波數(shù)的計(jì)算量.下面對(duì)并行計(jì)算所花的時(shí)間與非并行計(jì)算進(jìn)行對(duì)比分析.正演計(jì)算對(duì)比:對(duì)模型算例中未加入并行處理的程序進(jìn)行多次正演運(yùn)算計(jì)時(shí),平均計(jì)算時(shí)間為107s,而加入并行處理的程序,平均正演時(shí)間為35s,可以看出并行處理能夠節(jié)省三分之二的正演時(shí)間.反演時(shí)波數(shù)域雅克比矩陣計(jì)算對(duì)比:對(duì)模型算例中未加入并行處理的程序進(jìn)行波數(shù)域雅克比矩陣運(yùn)算計(jì)時(shí),平均計(jì)算時(shí)間約為641s,近10min,而加入并行處理的程序,平均計(jì)算的時(shí)間為173s,約3min,節(jié)省了近三分之二多的時(shí)間,因此本算例對(duì)波數(shù)域雅克比矩陣計(jì)算加入并行處理的效果同樣很突出.4.2.3拉格朗日乘子搜索方案改進(jìn)測(cè)試同上例一樣,對(duì)極化率反演過(guò)程中拉格朗日乘子搜索次數(shù)進(jìn)行計(jì)數(shù),常規(guī)流程中總搜索次數(shù)為55次,其中進(jìn)行尋找交叉點(diǎn)前的搜索次數(shù)為49次.對(duì)改進(jìn)方案的拉格朗日乘子搜索次數(shù)進(jìn)行計(jì)數(shù),其中總搜索次數(shù)為50次,開(kāi)始尋找交叉點(diǎn)前的搜索次數(shù)為44次.對(duì)于此算例的極化