一維電測深正反演程序

1、一、 電阻率測深法原理電阻率法是通過觀測地表的電場來了解地下介質電性分布的，要探測一定深度的地層的存在，必須使其明顯的影響到地表的電場分布，也就是要求其對觀測點處的電場有明顯的擾動，而要做到這一點，就要求流入相應深度的電流份額足夠多。因為在電阻率法中都是使用點電流源，因此需要考察距離點電流源不同距離的時透入某一給定深度以下供電電流所占比例的變化規(guī)律。在相距2L的兩個異性點電流源AB之間的中垂面上任意一點上的電流密度為：J=IL1(L2+y2+z2)32式中，y為觀測點距AB連線的水平距離；z為深度；I為供電電流強度。透入給定深度z以下的相對電流強度為：IzI=L-zdydz(L2+y2+z2)

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

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

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

5、，假定地面水平,在地下有n層水平層狀地層，各層的電阻率分別為1，2，n；厚度分別為h1,h2,hn，其中hn；每層地面到地面的距離為H1，H2，Hn,其中Hn。在A點有一電流源供電，電流強度為I。 用柱坐標系，將原點設在A點，z軸垂直向下。由于問題的解具有軸對稱性，與 無關,因此電位分布滿足以下形式的拉普拉斯方程2Ur2+1rUr+2Uz2=0 （3-1） 以及如下邊界條件（1） 電源點附近，趨于地面點電流源的正常電位，即U1R=r2+z20I12R （3-2）（2） 在地面處電流密度法向分量為零，即 11U1zz=0=0 （3-3） （3）除場源點外，電位處處有限，且無窮遠處電位為零，即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）可用分離變量法求解，設U=RrZz （3-6）經(jīng)分離變量后得到兩個二階常微分方程2Rr2+1rRr+m2R=0 （3-7）2Zz2-m2Z=0 （3-8） 其中，式（3-7）的解為第一類和第二類零階貝塞爾函數(shù)J0(mr)和Y0(mr)。第二類零階貝塞爾函數(shù)Y0(mr)在r=0的Z軸上趨于無窮，這種特征不符合點場源特征，因此應舍去，這樣式（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。 在第一層介質中的電位，還應附加電源點電位，即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層介質中么，當z時，電位有限，因此Bnm=0， Unr,z=I20Anme-mzJ0mrdm （3-14） 代入邊界條件后可得到地面電位的表達式U1r=I20T1mJ0mrdm （3-15） 式中，T1為電阻率轉換函數(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）水平地層上視電阻率表達式及濾波算法式（3-15）對微分，得電場強度 （3-17）由此得到時的對稱四極裝置視電阻率表達式 （3-18）其中是供電極距（AB/2

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

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

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

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

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

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

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

16、正演 已知模型參數(shù)：三層水平層狀模型，第一層電阻率為50，地層厚度為5；第二層電阻率為10，地層厚度為10；第三層電阻率為100。在程序中用dat.txt文件給出，文件內容（3 50 10 100 5 10）第一項表示3層數(shù)，后面3項為各層的電阻率、最后2項為各層厚度。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%水平層狀大地模型電阻率測深正演%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;運行程序所得電測深曲線如下：如圖所示，為取了18個不同的供電極距AB/2所得的離散化的視電阻率圖形。2） 反演用上面已知模型modelfile.txt所包含的模型參數(shù)3 50 10 100 5 10經(jīng)過正演所得的數(shù)據(jù)（保存在Ps.mat中）作為反演過程中已知 的觀測視電阻率數(shù)據(jù)，進行自動反演運算，反演

19、過程中設定的初始 模型保存在modelfile1.txt文件中，用行向量表示為3 3 3 3 3 3第1項為地層層數(shù)，第2、3、4項分別為各層電阻率，第5、6 項分別為各層厚度。MATLAB編程如下：%三層大地模型電測深自動反演% 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 err>0.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(i)/Ps1(i)-1);%求取系數(shù)矩陣 end end b=(