時間域激電法25維正演算法的研究

時間域激電法25維正演算法的研究

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