一維電測(cè)深正反演程序

1、一、 電阻率測(cè)深法原理電阻率法是通過(guò)觀測(cè)地表的電場(chǎng)來(lái)了解地下介質(zhì)電性分布的，要探測(cè)一定深度的地層的存在，必須使其明顯的影響到地表的電場(chǎng)分布，也就是要求其對(duì)觀測(cè)點(diǎn)處的電場(chǎng)有明顯的擾動(dòng)，而要做到這一點(diǎn)，就要求流入相應(yīng)深度的電流份額足夠多。因?yàn)樵陔娮杪史ㄖ卸际鞘褂命c(diǎn)電流源，因此需要考察距離點(diǎn)電流源不同距離的時(shí)透入某一給定深度以下供電電流所占比例的變化規(guī)律。在相距2L的兩個(gè)異性點(diǎn)電流源AB之間的中垂面上任意一點(diǎn)上的電流密度為：J=IL1(L2+y2+z2)32式中，y為觀測(cè)點(diǎn)距AB連線的水平距離；z為深度；I為供電電流強(qiáng)度。透入給定深度z以下的相對(duì)電流強(qiáng)度為：IzI=L-zdydz(L2+y2+z2)

2、32=1-2arctanzL （1-1）下圖所示為透入深度z以下空間的電流Iz/I隨L/z變化的情況。從圖中可以看出，當(dāng)L/z值比較小時(shí)，透入深度z以下空間的電流Iz/I比例也小，只能探測(cè)到近地表的情況；增大電測(cè)深的供電電極距L時(shí)，透入某一給定深度z以下的供電電流比例將隨之增大。當(dāng)L/z較大時(shí)，就可以探測(cè)到較深的部位。 在研究地下介質(zhì)電阻率的垂向變化時(shí)，希望盡量減小橫向電阻率變化的影響。如前所述，移動(dòng)測(cè)量電極MN對(duì)地下介質(zhì)電阻率的橫向變化反映非常明顯，而移動(dòng)供電電極AB對(duì)地下介質(zhì)電阻率的橫向變化反映則遠(yuǎn)沒(méi)有那么明顯。為了減少橫向電阻率變化的影響，應(yīng)該采用一種測(cè)量電極MN基本保持不動(dòng)，主要移動(dòng)供

3、電電極AB的裝置。在實(shí)際工作中，一般采用對(duì)稱四極測(cè)深裝置，在施工條件限制時(shí)，也可采用三極測(cè)深裝置，其他裝置則很少使用。二、 對(duì)稱四極測(cè)深裝置簡(jiǎn)介對(duì)稱四極測(cè)深裝置野外工作布置如下圖所示，供電電極AB和測(cè)量電極MN都以測(cè)點(diǎn)O為中心對(duì)稱布置在一條直線上。最初的供電電極距僅數(shù)米，逐步取一系列的遞增值，每個(gè)數(shù)量級(jí)距離供電極距改變約56次，各供電極距AB/2在對(duì)數(shù)軸上應(yīng)均勻分布（大致按照相同的倍數(shù)增大）。每一個(gè)供電極距與前一個(gè)供電極距的比值大約為1.21.5左右。選擇供電極距時(shí)，要求最小的極距應(yīng)能反映地表淺層電阻率，最大的極距則能滿足勘探深度要求，并保證測(cè)深曲線尾支的完整，不妨礙解釋最后一個(gè)電性層。從勘探

4、深度方面考慮，供電電極距AB/2應(yīng)從最小勘探深度的一半到最大勘探深度的5倍左右。測(cè)量電極MN開(kāi)始是固定的，例如取0.5m；直到（隨著供電電極距的加大）電壓過(guò)小時(shí)，才去另一增大值，例如3m，以此類推，一般MN的大小大約為AB的1/31/30 。在改變MN時(shí)一般要求有2個(gè)供電極距以2組MN極距觀測(cè)。因?yàn)樵龃鬁y(cè)量電極距MN會(huì)降低勘探深度，因此增大測(cè)量電極距時(shí)，s曲線通常會(huì)出現(xiàn)脫節(jié)現(xiàn)象。另外還有一種特殊的對(duì)稱四極裝置，它是始終保持MN=AB/3，稱為Wenner裝置，西方國(guó)家用的較多，這種裝置的s曲線是光滑的，沒(méi)有脫節(jié)問(wèn)題。三、 多層水平地層上視電阻率表達(dá)式1） 水平地層上地面點(diǎn)電流源的電場(chǎng)如下圖所示

5、，假定地面水平,在地下有n層水平層狀地層，各層的電阻率分別為1，2，n；厚度分別為h1,h2,hn，其中hn；每層地面到地面的距離為H1，H2，Hn,其中Hn。在A點(diǎn)有一電流源供電，電流強(qiáng)度為I。 用柱坐標(biāo)系，將原點(diǎn)設(shè)在A點(diǎn)，z軸垂直向下。由于問(wèn)題的解具有軸對(duì)稱性，與 無(wú)關(guān),因此電位分布滿足以下形式的拉普拉斯方程2Ur2+1rUr+2Uz2=0 （3-1） 以及如下邊界條件（1） 電源點(diǎn)附近，趨于地面點(diǎn)電流源的正常電位，即U1R=r2+z20I12R （3-2）（2） 在地面處電流密度法向分量為零，即 11U1zz=0=0 （3-3） （3）除場(chǎng)源點(diǎn)外，電位處處有限，且無(wú)窮遠(yuǎn)處電位為零，即Un

6、z0 （4） 在巖層分界面處電位連續(xù)，即Uiz=Hi=Ui+1z=Hi i=1，2，n-1 （3-4）（5） 在巖層分界面處電流密度法向分量連續(xù)，即1iUizz=Hi=1i+1Ui+1zz=Hi i=1，2，n-1 （3-5）式（3-1）可用分離變量法求解，設(shè)U=RrZz （3-6）經(jīng)分離變量后得到兩個(gè)二階常微分方程2Rr2+1rRr+m2R=0 （3-7）2Zz2-m2Z=0 （3-8） 其中，式（3-7）的解為第一類和第二類零階貝塞爾函數(shù)J0(mr)和Y0(mr)。第二類零階貝塞爾函數(shù)Y0(mr)在r=0的Z軸上趨于無(wú)窮，這種特征不符合點(diǎn)場(chǎng)源特征，因此應(yīng)舍去，這樣式（3-7）的解為第一類零

7、階貝塞爾函數(shù)J0mr。式（3-8）的解為 Zz=Ae-mz+Bemz （3-9） 由此可寫出各層電位積分形式的通解Uir,z=I20Aime-mz+BimemzJ0mrdm （3-10） 式中，Ai(m)和Bi(m)為待定函數(shù)；i=1,2，n。 在第一層介質(zhì)中的電位，還應(yīng)附加電源點(diǎn)電位，即U1r,z=I12R+I20A1me-mz+B1memzJ0mrdm （3-11）由U1(r,z)z=I1z2(r2+z2)32+I20-A1me-mz+B1memzmJ0mrdm 據(jù)邊界條件式（3-3），可得A1m=B1m 又根據(jù)李普希積分：1R=0e-mzJ0mrdm （3-12）得U1r,z=I201+

8、B1me-mz+B1memzJ0mrdm （3-13） 在第n層介質(zhì)中么，當(dāng)z時(shí)，電位有限，因此Bnm=0， Unr,z=I20Anme-mzJ0mrdm （3-14） 代入邊界條件后可得到地面電位的表達(dá)式U1r=I20T1mJ0mrdm （3-15） 式中，T1為電阻率轉(zhuǎn)換函數(shù)。注意到Tn=n，可得到T1的遞推公式為 Tn=n Ti=iTi+1+i+(Ti+1-i)e-2mhiTi+1+i-(Ti+1-i)e-2mhi （3-16） 2）水平地層上視電阻率表達(dá)式及濾波算法式（3-15）對(duì)微分，得電場(chǎng)強(qiáng)度 （3-17）由此得到時(shí)的對(duì)稱四極裝置視電阻率表達(dá)式 （3-18）其中是供電極距（AB/2

9、）在電測(cè)深視電阻率的表達(dá)式中的被積函數(shù)可以分為兩部分，一是包含地下各層所有信息（厚度及電阻率）的電阻率轉(zhuǎn)換函數(shù)，二是與地層參數(shù)無(wú)關(guān)的貝塞爾函數(shù)，雖然沒(méi)有解析計(jì)算結(jié)果，但可用線性濾波方法進(jìn)行計(jì)算。根據(jù)電阻率轉(zhuǎn)換函數(shù)的變化規(guī)律，對(duì)的抽樣取對(duì)數(shù)等間隔比較合適，因此，首先令，則電阻率的表達(dá)式變?yōu)?（3-19）根據(jù)采樣定理，一個(gè)函數(shù)可以用它在等間距離散抽樣點(diǎn)上的抽樣值表達(dá) （3-20）將電阻率裝換函數(shù)用它在數(shù)軸上的離散抽樣值表達(dá)為 （3-21）記 （3-22）則 （3-23）將預(yù)先計(jì)算出來(lái)，實(shí)際上取有限項(xiàng)求和就可以達(dá)到足夠的精度了。從頻譜分析的觀點(diǎn)看，當(dāng)電阻率裝換函數(shù)用它在數(shù)軸上的離散抽樣值表達(dá)式時(shí)，相

10、當(dāng)于濾取了頻率高于的諧波成分，因此這種計(jì)算視電阻率的方法稱為濾波計(jì)算方法，稱為濾波系數(shù)。為了提高線性濾波計(jì)算的精度，減少濾波系數(shù)的數(shù)目，需要對(duì)的抽樣點(diǎn)進(jìn)行位移，實(shí)際使用的線性濾波計(jì)算視電阻率的公式為 （3-24）式中，為供電極距為時(shí)的視電阻率；為時(shí)的電阻率轉(zhuǎn)換函數(shù)；為第個(gè)濾波系數(shù)；為抽樣間隔；為位移系數(shù)。 實(shí)際線性濾波計(jì)算常用的抽樣間隔有兩種。一種為六點(diǎn)是抽樣間隔，即，直流電阻率測(cè)深曲線一般都采用這種抽樣間隔進(jìn)行計(jì)算。另一種為十點(diǎn)式抽樣間隔，即，電磁測(cè)深曲線一般都采用這種抽樣間隔進(jìn)行計(jì)算。下表為常用的一套六點(diǎn)式抽樣間隔的濾波系數(shù)，共有20個(gè)系數(shù)：，其位移系數(shù)，。用上述公式就可計(jì)算某個(gè)供電極距的

11、視電阻率，只要計(jì)算20個(gè)對(duì)應(yīng)于這個(gè)供電極距的不同值的電阻率轉(zhuǎn)換函數(shù)，將其與下表中對(duì)應(yīng)的20個(gè)濾波系數(shù)相乘再求和就可以了。因?yàn)橐?jì)算一條視電阻率測(cè)深曲線就需要計(jì)算多個(gè)供電極距的視電阻率，為減少計(jì)算量，取，則，這樣，計(jì)算不同供電極距的視電阻率所需要的電阻率轉(zhuǎn)換函數(shù)大多數(shù)可以共用。這樣取值時(shí)，線性濾波計(jì)算視電阻率的公式可以寫為 （3-25）下面是用數(shù)值濾波法計(jì)算視電阻率測(cè)深曲線的計(jì)算機(jī)流程：（1） 輸入層參數(shù)，包括層數(shù)、各層的層厚度和電阻率；（2） 輸入要計(jì)算的供電極距范圍，并由此得到中的變化范圍：；（3） 計(jì)算的變化范圍：；（4） 用電阻率轉(zhuǎn)換函數(shù)遞推公式，循環(huán)計(jì)算，時(shí)的；（5）用濾波公式2-19

12、，循環(huán)計(jì)算，時(shí)的。四、自動(dòng)反演法理論反演是由響應(yīng)到模型的過(guò)程，是正演的逆過(guò)程，對(duì)水平層狀大地對(duì)稱四極電阻率測(cè)深曲線的解釋就是反演問(wèn)題，主要解釋方法有四種：特征點(diǎn)法、直接反演法、正演擬合法，自動(dòng)反演法。由于計(jì)算機(jī)的便利，我們通常采用自動(dòng)反演法。電阻率測(cè)深數(shù)據(jù)的自動(dòng)反演方法可歸結(jié)為尋找模型使下面的目標(biāo)函數(shù)趨于極?。?（4-1）式中， 是第個(gè)極距的實(shí)測(cè)視電阻率，是由模型正演計(jì)算所得的第個(gè)極距的理論視電阻率，是模型參數(shù)（地層的電阻率或厚度），是模型參數(shù)個(gè)數(shù)，是供電極距數(shù)，是范數(shù)，當(dāng)時(shí)即為最小二乘法。將模型在預(yù)測(cè)模型處展開(kāi)，并忽略二次項(xiàng)以上的項(xiàng)，式4-1表達(dá)可改為求預(yù)測(cè)模型修改量使目標(biāo)函數(shù)趨于極?。?（

13、4-2） 式4-2要趨于極小，則對(duì)于各供電極距要滿足下面的線性方程： （4-3） 解上述線性方程組，就可得到預(yù)測(cè)模型的修改量，于是可得到新的預(yù)測(cè)模型，計(jì)算新模型的理論視電阻率，與實(shí)測(cè)視電阻率進(jìn)行對(duì)比，如果精度滿足要求，則新的預(yù)測(cè)模型即為反演結(jié)果；否則重新展開(kāi)，計(jì)算模型修改量，直到滿足要求。下面是對(duì)反演過(guò)程的三個(gè)具體問(wèn)題的具體處理方法。（1）偏導(dǎo)數(shù)的計(jì)算?？捎貌罘址椒▉?lái)計(jì)算，取，則 （4-4） （2）模型參數(shù)的處理。模型參數(shù)中有不同量綱的電阻率和厚度，尤其對(duì)電阻率參數(shù)，變化范圍很大，如果直接由上述方法求解，不但會(huì)導(dǎo)致方程4-3嚴(yán)重病態(tài)，而且電阻率和厚度參數(shù)的修改量也不會(huì)正確，從而導(dǎo)致反演方法不收

14、斂。為了解決這個(gè)問(wèn)題，可對(duì)方程4-3無(wú)量綱化。將式4-4的偏導(dǎo)數(shù)代入方程式4-3中，有（4-5）上式兩端同時(shí)除以得 （4-6）令 這樣方程組4-6可寫為矩陣形式 （4-7）解線性方程組4-7可得新模型參數(shù)為 （4-8）另外為了防止模型參數(shù)修改過(guò)量，實(shí)際過(guò)程中又可作如下規(guī)定：時(shí)?。粫r(shí)取，即每次修改量不超過(guò)原有模型參數(shù)值的一半，保證收斂穩(wěn)定。（3）線性方程組的求解。用奇異值分解法求解，其基本思想是建立在奇異值分解定理上即任意階矩陣均可分解為，這里為階正交陣和為階正交陣，為的奇異值組成的階對(duì)角陣。當(dāng)非奇異時(shí)，奇異值較大，方程組3-7有廣義逆解，當(dāng)接近奇異時(shí)，有的奇異值就較小，此時(shí)由于中系數(shù)過(guò)大，上述

15、解的誤差就較大。為了解決這個(gè)問(wèn)題，維根斯（Wiggins）提出用最接近的矩陣來(lái)代替，其中是將中小的奇異值用零代替，因此有較精確的廣義逆解。 下面是計(jì)算機(jī)反演的具體流程：（1） 輸入實(shí)測(cè)視電阻率曲線；（2） 根據(jù)實(shí)測(cè)視電阻率曲線確定預(yù)測(cè)模型的層數(shù)；（3） 根據(jù)曲線特征，初步確定預(yù)測(cè)模型的層參數(shù)；（4） 調(diào)用正演程序計(jì)算預(yù)測(cè)模型的視電阻率曲線；（5） 對(duì)比和曲線，計(jì)算擬合誤差； （4-9）若擬合誤差滿足精度要求時(shí)，輸出預(yù)測(cè)模型參數(shù)，反演結(jié)束，若不滿足精度要求時(shí)，由方程3-6計(jì)算出模型修改量，得到新的預(yù)測(cè)模型，返回第四步，直到滿足精度要求。五、 MATLAB在三層大地模型電阻率測(cè)深正反演中的實(shí)現(xiàn)1）

16、正演 已知模型參數(shù)：三層水平層狀模型，第一層電阻率為50，地層厚度為5；第二層電阻率為10，地層厚度為10；第三層電阻率為100。在程序中用dat.txt文件給出，文件內(nèi)容（3 50 10 100 5 10）第一項(xiàng)表示3層數(shù)，后面3項(xiàng)為各層的電阻率、最后2項(xiàng)為各層厚度。MATLAB正演編程如下：%定義function子程序%function Ps=myfunction(n,P,h)T=P(n);C=0.003042 -0.001198 0.01284 0.0235 0.08688 0.2374 0.6194 1.1817 0.4248 -3.4507 2.7044 -1.1324 0.393

17、-0.1436 0.05812 -0.02521 0.01125 -0.004978 0.002072 -0.000318;ii=1:18;Ps=0;for kk=1:20 jj=(kk-ii(1):-1:(kk-ii(18); m=0.11396*10.(jj/6); for k=(n-1):-1:1 T=P(k)*(T+P(k)+(T-P(k).*exp(-2*m*h(k)./(T+P(k)-(T-P(k).*exp(-2*m*h(k); end Ps=Ps+T*C(kk);end%水平層狀大地模型電阻率測(cè)深正演%load modelfile.txt;M=modelfile;n=M(1);

18、P=M(2:(n+1);h=M(n+2):(2*n);ii=1:18;Ps=myfunction(n,P,h);r=10.(ii/6);loglog(r,Ps,.);xlabel(AB/2);ylabel(Ps);save Ps.mat;運(yùn)行程序所得電測(cè)深曲線如下：如圖所示，為取了18個(gè)不同的供電極距AB/2所得的離散化的視電阻率圖形。2） 反演用上面已知模型modelfile.txt所包含的模型參數(shù)3 50 10 100 5 10經(jīng)過(guò)正演所得的數(shù)據(jù)（保存在Ps.mat中）作為反演過(guò)程中已知 的觀測(cè)視電阻率數(shù)據(jù)，進(jìn)行自動(dòng)反演運(yùn)算，反演過(guò)程中設(shè)定的初始 模型保存在modelfile1.txt文件

19、中，用行向量表示為3 3 3 3 3 3第1項(xiàng)為地層層數(shù)，第2、3、4項(xiàng)分別為各層電阻率，第5、6 項(xiàng)分別為各層厚度。MATLAB編程如下：%三層大地模型電測(cè)深自動(dòng)反演% load modelfile1.txt load Ps.mat r=10.(ii/6); loglog(r,Ps,.);xlabel(AB/2);ylabel(Ps);hold on; M1=modelfile1; n1=M1(1); P1=M1(2:(n1+1); h1=M1(n1+2):(2*n1); Ps1=myfunction(n1,P1,h1); err=sqrt(sum(Ps-Ps1)./Ps).2)/18); while err0.05 for j=1:2*n1-1 mm=M1(2:2*n1); mm(j)=1.1*mm(j); P1=mm(1:n1); h1=mm(n1+1):(2*n1-1); Ps2=myfunction(n1,P1,h1); for i=1:18 a(i,j)=10*(Ps2(