本發(fā)明涉及數(shù)據(jù)軟處理的計(jì)算機(jī)技術(shù)，特別是最小二乘法電阻率三維近似反演技術(shù)。

背景技術(shù)：

近年來電阻率法發(fā)展迅速,在儀器研制及實(shí)際應(yīng)用方面均取得不少成果。然而實(shí)際探測目標(biāo)多表現(xiàn)為三維電性結(jié)構(gòu)，spitzer通過電阻率法靈敏度矩陣揭示出：對三度體用電剖面數(shù)據(jù)做二維解釋不可避免會有偏差，所以，電阻率三維測量及相應(yīng)反演方法的研究成為該領(lǐng)域的前沿課題。

地球物理反問題是不適定的，二、三維地電場反演更是如此。因?yàn)槠湔菔腔谟邢薏罘只蛴邢拊活愑?jì)算方法的數(shù)值解。反演是求數(shù)量巨大的各網(wǎng)格剖分單元電阻率值。因此，反演過程中不可避免地要遇到如何求偏導(dǎo)數(shù)矩陣、如何解決計(jì)算速度極慢以及龐大機(jī)器內(nèi)存需求的問題。另外，由于反演參數(shù)太多導(dǎo)致反演不穩(wěn)定性和高度非唯一性等諸多問題都給反演帶來巨大困難。盡管如此，隨著計(jì)算機(jī)及數(shù)值計(jì)算技術(shù)的發(fā)展，近年來國內(nèi)外已有不少二維地電場反演結(jié)果發(fā)表，至于電阻率三維反演，國外亦起步不久。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于提供一種最小二乘法電阻率三維近似反演技術(shù)，計(jì)算量較低，加快反演計(jì)算速度，節(jié)省機(jī)器內(nèi)存，改善反演的穩(wěn)定性與可靠性，有利于較細(xì)網(wǎng)格參數(shù)單元的復(fù)雜模型反演。

本發(fā)明的目的是這樣實(shí)現(xiàn)的：一種最小二乘法電阻率三維近似反演技術(shù)：

步驟⑴-建立立體單元網(wǎng)格：設(shè)被探測區(qū)域的地下空間為在局部均勻分布的層狀地電斷面，物探測線垂直于構(gòu)造走向，以瞬變電磁測點(diǎn)為中心將整個(gè)被探測區(qū)域表面劃分成n個(gè)nx×ny的二維網(wǎng)格，同時(shí)在垂向上按地下空間最大被探測深度將整個(gè)地下空間劃分成若干段，使整個(gè)地下空間劃分成n個(gè)nx×ny×nz的三維網(wǎng)格，使之形成有限單元法計(jì)算網(wǎng)格；采用高密度電阻率法從深度小于xm的淺部火燒區(qū)及采空區(qū)采集淺部電阻率數(shù)值；nx，ny，nz依次代表在三維正交坐標(biāo)x方向、y方向與z方向上的等分?jǐn)?shù)量值；

步驟⑵-建立線性方程組：根據(jù)從被探測區(qū)域的鉆孔被測得的電阻率數(shù)值，取得被探測區(qū)域內(nèi)某一點(diǎn)s的實(shí)測電阻率值，然后對電阻率數(shù)值進(jìn)行劃分以獲得m個(gè)實(shí)測電阻率數(shù)值，繼而獲取電阻率序列{psi}＝[ps1,ps2……psm]，采用瞬變電磁法測得每個(gè)測點(diǎn)視電導(dǎo)率序列[δ]＝[δ1,δ2……δm]，根據(jù)已知電阻率數(shù)列及視電導(dǎo)率數(shù)列形成對應(yīng)于被探測區(qū)域的m×n階jacobi矩陣；

步驟⑶-⑴采用阻尼最小二乘法與共軛梯度法求解jacobi矩陣，借助公式①及公式②獲得反演模型的光滑解：

公式①：(g^tg+λc^tc)δm＝g^tδd，

公式②：

在公式①或公式②中，g是jacobi矩陣，g^t是對應(yīng)于g的轉(zhuǎn)置矩陣，λ是拉格朗日因子，c是模型光滑矩陣，c^t是對應(yīng)于c的轉(zhuǎn)置矩陣，δm為修改量，δd是殘差向量；在公式②中，ψ是目標(biāo) 函數(shù)，第一項(xiàng)是以差分代替一階微分后模型粗糙度的離散表達(dá)式，粗糙度定義為δm在x、y、z方向一階偏微分的平方，rx、ry、rz分別是反演模型在x、y、z方向的粗糙度矩陣；

步驟⑷-將從jacobi矩陣求解而獲得的三維空間電阻率數(shù)值轉(zhuǎn)換成以nx×ny×nz為單元的電阻率三維數(shù)據(jù)體。

在地球物理反演中，傳統(tǒng)最小二乘迭代方法得到非常廣泛、有效應(yīng)用。電阻率三維反演的一般形式可表示為：δd＝gδm，g為jacobi矩陣；δd為觀測數(shù)據(jù)d^obs與正演理論值d0的殘差向量；δm為初始模型m0的修改向量。對于三維反演問題，可將模型剖分成nx×ny×nz的三維網(wǎng)格，反演要求的參數(shù)就是各網(wǎng)格單元內(nèi)的電導(dǎo)率σ值，三維反演的觀測數(shù)據(jù)則是e-scan測量的單極—單極電位值由于它們變化范圍大，一般用對數(shù)來標(biāo)定反演數(shù)據(jù)及模型參數(shù)，即及m＝lnσ，有利于改善反演的穩(wěn)定性。由于反演參數(shù)太多，傳統(tǒng)的阻尼最小二乘反演往往導(dǎo)致過于復(fù)雜的模型，即產(chǎn)生所謂多余構(gòu)造，它是數(shù)據(jù)本身所不要求的或是不可分辨的構(gòu)造信息，給解釋帶來困難。本發(fā)明主要被用于探測和更準(zhǔn)確地定位煤礦采空區(qū)及火燒區(qū)含水范圍。由于電阻率勘探廣泛應(yīng)用于石油、基礎(chǔ)建設(shè)、金屬礦產(chǎn)等，本發(fā)明還可被應(yīng)用到非煤礦山領(lǐng)域。本發(fā)明能有效達(dá)到近似反演的技術(shù)目標(biāo)，計(jì)算量較低，加快反演計(jì)算速度，節(jié)省機(jī)器內(nèi)存，改善反演的穩(wěn)定性與可靠性，有利于較細(xì)網(wǎng)格參數(shù)單元的復(fù)雜模型反演。

說明書附圖

圖1為對應(yīng)于未經(jīng)本發(fā)明處理過的電阻率數(shù)值的原始可視數(shù)據(jù)例圖；

圖2為對應(yīng)于已經(jīng)本發(fā)明處理過的電阻率數(shù)值的特種可視反演例圖；

圖3為對應(yīng)于已經(jīng)本發(fā)明處理過的電阻率三維數(shù)據(jù)體的可視三維資料解釋例圖。

具體實(shí)施方式

一種最小二乘法電阻率三維近似反演技術(shù)：

步驟⑴-建立立體單元網(wǎng)格：設(shè)被探測區(qū)域的地下空間為在局部均勻分布的層狀地電斷面，物探測線垂直于構(gòu)造走向，以瞬變電磁測點(diǎn)為中心將整個(gè)被探測區(qū)域表面劃分成n個(gè)nx×ny的二維網(wǎng)格，同時(shí)在垂向上按地下空間最大被探測深度將整個(gè)地下空間劃分成若干段，使整個(gè)地下空間劃分成n個(gè)nx×ny×nz的三維網(wǎng)格，使之形成有限單元法計(jì)算網(wǎng)格；采用高密度電阻率法從深度小于xm的淺部火燒區(qū)及采空區(qū)采集淺部電阻率數(shù)值；nx，ny，nz依次代表在三維正交坐標(biāo)x方向、y方向與z方向上的等分?jǐn)?shù)量值；

步驟⑵-建立線性方程組：根據(jù)從被探測區(qū)域的鉆孔被測得的電阻率數(shù)值，取得被探測區(qū)域內(nèi)某一點(diǎn)s的實(shí)測電阻率值，然后對電阻率數(shù)值進(jìn)行劃分以獲得m個(gè)實(shí)測電阻率數(shù)值，繼而獲取電阻率序列{psi}＝[ps1,ps2……psm]，采用瞬變電磁法測得每個(gè)測點(diǎn)視電導(dǎo)率序列[δ]＝[δ1,δ2……δm]，根據(jù)已知電阻率數(shù)列及視電導(dǎo)率數(shù)列形成對應(yīng)于被探測區(qū)域的m×n階jacobi矩陣；

步驟⑶-⑴采用阻尼最小二乘法與共軛梯度方法求解jacobi矩陣，借助公式①及公式②獲得反演模型的光滑解：

公式①：(g^tg+λc^tc)δm＝g^tδd，

公式②：

在公式①或公式②中，g是jacobi矩陣，g^t是對應(yīng)于g的轉(zhuǎn)置矩陣，λ是拉格朗日因子，c是模型光滑矩陣，c^t是對應(yīng)于c的轉(zhuǎn)置矩陣，δm為修改量，δd是殘差向量；在公式②中，ψ是目標(biāo)函數(shù)，第一項(xiàng)是以差分代替一階微分后模型粗糙度(roughness)的離散表達(dá)式，粗糙度定義為δm在x、y、z方向一階偏微分的平方，rx、ry、rz分別是反演模型在x、y、z方向的粗糙度矩陣；

步驟⑷-將從jacobi矩陣求解而獲得的三維空間電阻率數(shù)值轉(zhuǎn)換成以nx×ny×nz為單元的電阻率三維數(shù)據(jù)體。對應(yīng)于已經(jīng)本發(fā)明處理過的電阻率三維數(shù)據(jù)體的可視三維資料解釋例圖如圖3所示。

如圖1所示，圖1的瞬變電磁數(shù)據(jù)在橫向分辨率高，但在縱向上與地質(zhì)構(gòu)造形態(tài)一致性差且深度誤差大。如圖2所示，圖2顯示的深度與已知鉆探深度吻合，地質(zhì)構(gòu)造形態(tài)清晰，水體的位置、形態(tài)均得到良好呈現(xiàn)，水位線深度明顯，但不足之處在于平面分辨率降低。因此，在實(shí)際應(yīng)用中，與原始可視數(shù)據(jù)例圖顯示特性類似的未被反演圖和與特種可視反演例圖顯示特性類似的被反演圖需以相結(jié)合的方式被應(yīng)用。

采用阻尼最小二乘法求解jacobi矩陣，在最小二乘準(zhǔn)則中加入光滑約束，以求得光滑模型，提高解的穩(wěn)定性。模型修改量δm的算法為：(g^tg+λc^tc)δm＝g^tδd，其中c是模型光滑矩陣。然而，傳統(tǒng)三維反演技術(shù)不能避免jacobi矩陣g及大型矩陣逆的巨量計(jì)算，因此，無論從計(jì)算速度還是機(jī)器內(nèi)存需求而言，均只能適宜較少單元的簡單模型反演。共軛梯度(cg)方法是解大型最優(yōu)化問題最有效方法之一。它直接從反演目標(biāo)函數(shù)出發(fā)，無須直接解上式的正則化方程，可避免偏導(dǎo)數(shù)矩陣g的計(jì)算和存儲，提高計(jì)算速度、節(jié)省機(jī)器內(nèi)存。本發(fā)明引入共軛梯度(cg)方法求電阻率三維最小構(gòu)造反演，目標(biāo)函數(shù)ψ為：

目標(biāo)函數(shù)ψ中的第一項(xiàng)是以差分代替一階微分后模型粗糙度(roughness)的離散形式，粗糙度定義為δm在x、y、z方向一階偏微分的平方和rx、ry、rz分別是模型在x、y、z方向的粗糙度矩陣，λ是拉格朗日因子。求目標(biāo)函數(shù)ψ極小實(shí)際是在觀測數(shù)據(jù)的約束下，旨在使模型的粗糙度極小，亦即求模型的光滑解。引入共軛梯度(cg)方法求目標(biāo)函數(shù)ψ極小,可簡單導(dǎo)出求最小構(gòu)造反演流程；令：k＝1，2，……nmax^(inv)(反演迭代)；(1)給定初始模型mk；(2)計(jì)算模型mk正演理論值d0(k)；(3)計(jì)算殘差向量：δdk＝d^(obs)-d0(k)；(4)計(jì)算右端項(xiàng)：bk＝g^tδdk，(只需計(jì)算，而y＝δdk)；(5)cg迭代初始化：δmk(0)＝0；r(0)＝bk；p(0)＝r(0)，再令：ii＝1，2，……nmax^(cg)(cg迭代)，可得：

其中，y＝gx，x＝p⁽ⁱⁱ⁾，

β⁽ⁱⁱ⁾＝(r⁽ⁱⁱ⁺¹⁾，r⁽ⁱⁱ⁺¹⁾)/(r⁽ⁱⁱ⁾，r⁽ⁱⁱ⁾)，

p⁽ⁱⁱ⁺¹⁾＝r⁽ⁱⁱ⁺¹⁾+β⁽ⁱⁱ⁾p⁽ⁱⁱ⁾。

cg迭代結(jié)束后，δmk即可被求得。修改初始模型重新迭代至反演收斂：δmk+1＝mk+δmk，此即為利用共軛梯度(cg)方法的最小構(gòu)造反演算法,其中的nmax^(lnv)、nmax^(gg)分別是給定的最大反演迭代次數(shù)和cg迭代次數(shù)。本發(fā)明無需直接求g及g^tg的逆矩陣，只需g與任一向量x的乘積gx及g^t與任一向量y的乘積針對e-scan數(shù)據(jù)的電阻率三維反演,由cg迭代可求得δmk，然后修改初始模型mk，重新迭代直到滿足收斂標(biāo)準(zhǔn)，反演迭代自此結(jié)束。這樣可實(shí)現(xiàn)每次非線性反演迭代只需一次正演計(jì)算，大大小于傳統(tǒng)反演技術(shù)的計(jì)算量，加快反演計(jì)算速度，避免g及g^tg數(shù)據(jù)存儲量巨大的技術(shù)問題，有利于較細(xì)網(wǎng)格參數(shù)單元的復(fù)雜模型反演，而且在理論上如此計(jì)算出的向量gx及的精度與正演精度一樣，反演的穩(wěn)定性和可靠性也得到保證。