• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    考慮電流型畸變的大地電磁三維反演

    2016-07-28 07:04:10李鑫白登海閆永利馬曉冰孔祥儒
    地球物理學(xué)報(bào) 2016年6期

    李鑫, 白登海, 閆永利, 馬曉冰, 孔祥儒

    1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所, 北京 100029 2 中國(guó)科學(xué)院大學(xué), 北京 100039

    ?

    考慮電流型畸變的大地電磁三維反演

    李鑫1,2, 白登海1, 閆永利1, 馬曉冰1, 孔祥儒1

    1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所, 北京1000292 中國(guó)科學(xué)院大學(xué), 北京100039

    摘要大地電磁測(cè)深(MT)的觀測(cè)數(shù)據(jù)易受到由近地表小尺度非均勻體或地形起伏引起的電流型畸變干擾,消除或壓制這種干擾對(duì)獲取可靠的深部電性結(jié)構(gòu)至關(guān)重要.當(dāng)區(qū)域結(jié)構(gòu)為二維時(shí),電流型畸變可采用張量分解等方法予以消除或壓制.當(dāng)區(qū)域結(jié)構(gòu)為三維時(shí),畸變問(wèn)題更加復(fù)雜和嚴(yán)重,傳統(tǒng)張量分解方法往往效果不佳或無(wú)效,嚴(yán)重地制約了MT三維反演技術(shù)的實(shí)用性.對(duì)此,本文提出一種考慮電流型畸變的MT三維反演算法,將完整的電流型畸變參數(shù)引入到目標(biāo)函數(shù),并采用非線性共軛梯度法與電阻率參數(shù)同時(shí)反演,從而達(dá)到壓制畸變的目的.該算法有兩個(gè)關(guān)鍵點(diǎn):一是通過(guò)分析實(shí)測(cè)數(shù)據(jù)所遭受畸變的分布特征,在目標(biāo)函數(shù)中對(duì)其進(jìn)行有效約束;二是在迭代過(guò)程中,通過(guò)自適應(yīng)地調(diào)整雙正則化因子保障算法的穩(wěn)定和效率.理論模型測(cè)試結(jié)果顯示,常規(guī)三維反演算法不能合理解釋數(shù)據(jù)中的畸變成分,而只能通過(guò)引入虛假異常體強(qiáng)制地?cái)M合受畸變數(shù)據(jù),從而造成電阻率模型嚴(yán)重失真.與之相比,本文算法能夠在反演中自動(dòng)求解各測(cè)點(diǎn)所受到的畸變,獲得更接近真實(shí)的電阻率模型.

    關(guān)鍵詞電流型畸變; 靜態(tài)效應(yīng); 三維反演; 大地電磁測(cè)深

    1引言

    大地電磁測(cè)深(MT)是一種探測(cè)地球深部電性結(jié)構(gòu)的有效方法,它通過(guò)在地球表面同步觀測(cè)天然電、磁場(chǎng)分量,并在頻率域根據(jù)各分量之間的線性關(guān)系獲得多種響應(yīng)函數(shù),如阻抗張量(視電阻率/相位)、傾子等,進(jìn)而通過(guò)反演手段尋找能夠滿足觀測(cè)數(shù)據(jù)并符合物理規(guī)律的電阻率模型(Simpson and Bahr, 2005).然而,實(shí)際觀測(cè)的MT數(shù)據(jù)極易受到電流型畸變(Galvanic distortion)的干擾.若不對(duì)該畸變進(jìn)行妥善處理而貿(mào)然開(kāi)展反演,將難以獲得真實(shí)的地電模型,甚至因此導(dǎo)致錯(cuò)誤的解釋(Jiracek, 1990; 晉光文和孔祥儒, 2006).

    過(guò)去數(shù)十年間,國(guó)內(nèi)外同行針對(duì)如何消除或壓制電流型畸變已開(kāi)展了大量的理論分析和技術(shù)實(shí)踐,提出了諸多解決方案,其中被廣泛使用的是阻抗張量分解方法(Swift, 1967; Bahr, 1988; Groom and Bailey, 1989,1991; Groom and Bahr, 1992; McNeice and Jones,2001).這類方法通常假設(shè)區(qū)域結(jié)構(gòu)可簡(jiǎn)化為一、二維結(jié)構(gòu),進(jìn)而根據(jù)阻抗張量/畸變張量的數(shù)學(xué)/物理特征將畸變從阻抗張量中分解出來(lái),獲得基本不受畸變干擾的正常阻抗張量.然而,在實(shí)際運(yùn)用中,這類方法受到兩方面的嚴(yán)重制約:首先是其僅適用于較簡(jiǎn)單區(qū)域結(jié)構(gòu)情況(一、二維),其次是其僅能求解部分畸變參數(shù),而對(duì)于另外一部分關(guān)鍵的畸變參數(shù)(如靜位移因子)則無(wú)法求取.靜態(tài)位移效應(yīng)會(huì)造成觀測(cè)視電阻率曲線發(fā)生整體性偏移,一般可采用曲線平移校正、標(biāo)志層校正、空間濾波、輔助瞬變電磁或大極距直流電測(cè)深等方法予以應(yīng)對(duì)(Jones, 1988; Pellerin and Hohmann, 1990; Sternberg et al., 1988).這些方法在特定環(huán)境下具有一定效果,但是均存在不足.比如,很難找到一個(gè)普遍適用準(zhǔn)則來(lái)進(jìn)行曲線平移,而標(biāo)志層校正只適用于盆地等特定的地質(zhì)構(gòu)造環(huán)境,其他空間濾波、瞬變電磁或大極距直流電測(cè)深等均需要增加額外的觀測(cè).為此,DeGroot-Hedlin (1995)創(chuàng)造性地將靜位移和構(gòu)造走向參數(shù)同時(shí)引入到二維OCCAM反演中,與電阻率模型一并求解以消除其干擾,取得了較好的效果.

    近年來(lái)MT三維正反演技術(shù)發(fā)展迅速,逐漸開(kāi)始被用于一些實(shí)際資料解釋中(Zhang et al., 2015;Xiao et al., 2015).在三維區(qū)域結(jié)構(gòu)下,電流型畸變對(duì)MT觀測(cè)數(shù)據(jù)所造成的影響更加復(fù)雜和嚴(yán)重,上述現(xiàn)有的張量分解或靜校正方法都不再適用,因此發(fā)展一套能夠配合MT三維反演的畸變消除技術(shù)很有必要.已有一些學(xué)者對(duì)此開(kāi)展了研究,如Utada和Munekane(2000)利用水平電場(chǎng)與垂直磁場(chǎng)的空間導(dǎo)數(shù)關(guān)系壓制三維區(qū)域結(jié)構(gòu)下的電流型畸變干擾.Garcia和Jones (2002)在GB分解的基礎(chǔ)上(Groom and Bailey, 1989),假設(shè)相鄰測(cè)點(diǎn)遭受不同的畸變影響但反映相同的區(qū)域電性結(jié)構(gòu),利用牛頓法求解三維區(qū)域結(jié)構(gòu)中的畸變參數(shù).Bibby等 (2005)利用不受電流型畸變影響的相位張量獲得地下維度信息,進(jìn)而采用低維度(一,二維)數(shù)據(jù)求解畸變參數(shù).Patro等(2013)直接三維反演了相位張量數(shù)據(jù),有效壓制了畸變干擾,然而由于相位張量數(shù)據(jù)中不包含幅值信息導(dǎo)致初始電阻率模型的選取對(duì)反演結(jié)果影響很大.Sasaki和Meju (2006)參考DeGroot-Hedlin (1995)的思路,試圖將靜位移參數(shù)與三維電阻率結(jié)果同時(shí)反演.然而,在三維區(qū)域結(jié)構(gòu)中,由于阻抗張量各分量的模值和相位均受到畸變干擾,僅僅采用靜位移參數(shù)來(lái)描述電流畸變是不夠的 (Jones,2011).針對(duì)這一缺陷,Avdeeva等(2015)將完整的電流型畸變參數(shù)引入到三維反演中,采用擬牛頓法將其與電阻率參數(shù)同步反演,取得了較好的效果.

    本文首先給出了三維區(qū)域結(jié)構(gòu)下電流型畸變作用于區(qū)域阻抗張量的數(shù)學(xué)表達(dá)形式,并依據(jù)大量實(shí)測(cè)數(shù)據(jù)分析了電流畸變各參數(shù)的大致分布特征.根據(jù)上述信息,建立了新的目標(biāo)函數(shù),進(jìn)一步采用非線性共軛梯度算法同時(shí)反演了三維電阻率結(jié)構(gòu)和電流型畸變參數(shù).三維理論模型測(cè)試表明,該算法可以有效壓制電流型畸變的影響,獲得更接近真實(shí)的電阻率模型.

    2電流型畸變

    2.1電流型畸變的物理涵義

    電流型畸變是由近地表小尺度非均勻體電性界面處的電荷積累所產(chǎn)生的散射場(chǎng)與區(qū)域結(jié)構(gòu)所產(chǎn)生的正常場(chǎng)矢量疊加所造成的(Berdichevsky and Dmitriev,1976).假設(shè)三維區(qū)域背景的電導(dǎo)率分布為σ0(r),而近地表分布著多個(gè)局部三維散射體,其中第k個(gè)散射體的電導(dǎo)率為σk(r),它與區(qū)域電導(dǎo)率的差異為δσk(r)=σk(r)-σ0(r).Chave和Smith(1994)推導(dǎo)了第k個(gè)散射體內(nèi)外任意位置r處的電場(chǎng)E(r)的積分方程解為

    =E0(r)+eI(r)+eG(r).

    (1)

    式中r′為該散射體內(nèi)任意點(diǎn)的空間位置,E(r′)為r′處的電場(chǎng)矢量,Vk是第k個(gè)散射體的體積,g(r,r′)為標(biāo)量全空間格林函數(shù),積分對(duì)近地表所有局部散射體所占據(jù)的區(qū)域進(jìn)行.右側(cè)第一項(xiàng)E0(r)通常被稱為正常場(chǎng)或區(qū)域場(chǎng),是由我們感興趣的區(qū)域結(jié)構(gòu)所產(chǎn)生的電場(chǎng).右側(cè)剩余的兩項(xiàng)分別為近地表局部非均勻體產(chǎn)生的感應(yīng)型畸變電場(chǎng)eI(r)和電流型畸變電場(chǎng)eG(r).這兩個(gè)畸變場(chǎng)的強(qiáng)度取決于近地表局部異常體的體積、形狀及其與區(qū)域電導(dǎo)率的差異情況.采用波恩近似,可將(1)式積分內(nèi)的近地表散射體內(nèi)部的電場(chǎng)E(r′)替換為區(qū)域電場(chǎng)E0(r),得到觀測(cè)電場(chǎng)E(r)與區(qū)域電場(chǎng)E0(r)之間的關(guān)系:

    =c(r)E0(r).

    (2)

    上式中G(r)為格林并矢張量函數(shù).對(duì)于常規(guī)大地電磁的觀測(cè)頻段(尤其在低頻段),感應(yīng)型畸變迅速減小,通常忽略不計(jì).因此,畸變張量c(r)一般可以簡(jiǎn)化為一個(gè)與頻率無(wú)關(guān)的二階實(shí)數(shù)張量.對(duì)于磁場(chǎng),一般認(rèn)為在低頻段局部電荷積累對(duì)觀測(cè)磁場(chǎng)B(r)的影響一般較小(Chave and Smith, 1994),不予考慮.將(2)式代入阻抗張量的一般表達(dá)式,可以得到觀測(cè)阻抗張量D(r)與區(qū)域阻抗張量Z(r)二者的關(guān)系式如下:

    D(r)=E(r)B(r)-1=c(r)E0(r)B0(r)-1

    =c(r)Z(r),

    (3)

    上式中D(r)和Z(r)均為二階復(fù)數(shù)張量,隨頻率而變化.當(dāng)觀測(cè)數(shù)據(jù)不受電流型畸變影響時(shí),c(r)退化為單位矩陣I(晉光文和孔祥儒, 2006).

    2.2三維區(qū)域結(jié)構(gòu)中的電流型畸變

    當(dāng)區(qū)域電性結(jié)構(gòu)是三維時(shí),區(qū)域阻抗張量Z的對(duì)角、非對(duì)角分量均不為零,將(3)式展開(kāi)有

    (4)

    很明顯,觀測(cè)阻抗張量D各分量的幅值和相位均受到了畸變張量c的影響.上式中包含畸變張量c中4個(gè)和區(qū)域阻抗張量Z中8個(gè)(實(shí)、虛部各4個(gè))共12個(gè)待求解參數(shù),在不作假設(shè)的情況下,無(wú)法根據(jù)觀測(cè)阻抗張量D中8個(gè)(實(shí)、虛部各4個(gè))已知量建立的方程組直接求解這個(gè)非線性的欠定問(wèn)題.

    2.3電流型畸變參數(shù)的分布特征

    2.1節(jié)中已證實(shí)MT數(shù)據(jù)受到的電流型畸變通常只與觀測(cè)位置附近的淺層環(huán)境相關(guān),因此具有很強(qiáng)的隨機(jī)性.Degroot-Hedlin (1995)在將靜位移因子引入到二維反演中時(shí),采用了對(duì)數(shù)域內(nèi)所有測(cè)點(diǎn)所受靜態(tài)位移之和為零作為約束,測(cè)試結(jié)果表明當(dāng)測(cè)點(diǎn)間距較小且數(shù)量較多時(shí),這種約束是大致合理的.Sasaki和Meju (2006)依據(jù)不同地質(zhì)環(huán)境中采集的MT和瞬變電磁(TEM)數(shù)據(jù)分析了靜態(tài)位移參數(shù)的分布特點(diǎn),認(rèn)為其大體呈高斯分布,在高度風(fēng)化或沉積低阻蓋層地區(qū)表現(xiàn)為負(fù)靜態(tài)位移(視電阻率曲線下移),而在火山巖分布地區(qū)表現(xiàn)為正靜態(tài)位移(視電阻率曲線上移).Avdeev等(2015)根據(jù)經(jīng)驗(yàn)假設(shè)電流型畸變張量的對(duì)角分量大體以1為中心呈對(duì)稱分布,非對(duì)角分量大體以0為中心呈對(duì)稱分布.需要指出的是,上述假設(shè)均缺乏充分的實(shí)測(cè)數(shù)據(jù)的證實(shí).

    基于美國(guó)臺(tái)陣項(xiàng)目USARRAY(www.usarray.org)所采集的587個(gè)測(cè)點(diǎn)的高質(zhì)量MT數(shù)據(jù),我們采用相位張量分解方法(Caldwell et al., 2004;Bibby et al., 2005)估算了所有測(cè)點(diǎn)遭受的畸變情況,統(tǒng)計(jì)結(jié)果如圖1所示.這些觀測(cè)數(shù)據(jù)覆蓋范圍廣,地質(zhì)構(gòu)造環(huán)境豐富,具有較好的代表性.分析結(jié)果表明,畸變張量的非對(duì)角分量(c12,c21)取值以0為中心呈大致對(duì)稱;對(duì)角分量(c11,c22)取值主要集中在1的附近,但在1的右側(cè)具有較長(zhǎng)的尾部,在1的左側(cè)出現(xiàn)的次數(shù)迅速下降.

    圖1 USARRAY實(shí)測(cè)數(shù)據(jù)所受的電流型畸變參數(shù)值分布非對(duì)角分量(c12,c21)取值以0為中心呈近對(duì)稱分布;對(duì)角分量(c11,c22)在1附近集中出現(xiàn), 但在1的右側(cè)具有更長(zhǎng)的尾部.Fig.1 Histograms of the galvanic distortion parameters of the USARRAY MT dataThe value of off-diagonal elements symmetrically distributed around zero, while the diagonal elements are concentrated around one but have a longer tail on the right side.

    3考慮電流型畸變的三維反演算法

    3.1目標(biāo)函數(shù)

    與其他地球物理反演問(wèn)題一樣,能夠擬合一組MT觀測(cè)數(shù)據(jù)的電阻率模型和電流型畸變參數(shù)都不是唯一的,因此需要根據(jù)這兩類參數(shù)各自的特征在反演過(guò)程中進(jìn)行有效的約束.基于常規(guī)MT反演的目標(biāo)函數(shù),我們重新定義了一個(gè)考慮電流型畸變的目標(biāo)函數(shù):

    (5)

    右側(cè)第一項(xiàng)Φd為觀測(cè)數(shù)據(jù)擬合項(xiàng),其具體形式為

    (6)

    右側(cè)第二項(xiàng)Φm為電阻率模型約束項(xiàng),定義為

    (7)

    右側(cè)第三項(xiàng)Φc為電流型畸變參數(shù)約束項(xiàng),c與cprior分別為預(yù)測(cè)和先驗(yàn)畸變參數(shù),在反演開(kāi)始前同樣需要為c設(shè)置一個(gè)初始狀態(tài),即初始畸變參數(shù)c0;λ2為該項(xiàng)的正則化因子.從統(tǒng)計(jì)學(xué)的觀點(diǎn)來(lái)看,對(duì)于大致滿足高斯分布的數(shù)據(jù)采用L2范數(shù)解比較合理(王家映, 1998).上文2.3節(jié)分析已表明畸變張量的對(duì)角、非對(duì)角分量大致分別以1和0呈近似對(duì)稱高斯分布,因此我們采用了預(yù)測(cè)和先驗(yàn)畸變參數(shù)二者之差的L2范數(shù)這種約束形式,即

    (8)

    在沒(méi)有精確的先驗(yàn)信息情況下,同樣可根據(jù)畸變張量的分布特征將cprior設(shè)為單位矩陣I.

    3.2非線性共軛梯度反演算法

    共軛梯度法(CG)最早由Stiefel(1951)和Hestenes(1951)分別獨(dú)立提出,Hestenes和Stiefel(1952)給出了用CG方法求解正定系數(shù)矩陣的線性方程組的詳細(xì)過(guò)程.經(jīng)過(guò)Fletcher和Reeves(1964)發(fā)展后,CG被用于解決非線性最優(yōu)化問(wèn)題.Mackie和Madden(1993)首先在三維MT反演中采用了CG方法.Rodi和Mackie (2001)發(fā)展了非線性共軛梯度法(NLCG),并將其成功用于MT二維反演中.由于無(wú)需直接計(jì)算和存儲(chǔ)雅克比矩陣且收斂速度較快,NLCG很快得到了廣泛應(yīng)用.Newman和Alumbaugh(2000,注:該工作實(shí)際在Rodi和Mackie (2001)之后)首先把NLCG方法用于三維MT反演中,有效降低了三維反演的計(jì)算成本.

    本文算法以目標(biāo)函數(shù)(5)式關(guān)于電阻率模型m和電流型畸變參數(shù)c的負(fù)梯度方向分別作為二組參數(shù)的初始搜索方向,之后在每一次迭代中以當(dāng)前的兩組Polak-Ribiere共軛梯度方向作為搜索前進(jìn)方向(Polak and Ribière, 1969).在每一次迭代中,通過(guò)三次多項(xiàng)式插值確定兩組最優(yōu)搜索步長(zhǎng),按照兩組搜索方向和步長(zhǎng)更新兩組參數(shù),并通過(guò)Armijo條件判斷當(dāng)前迭代中目標(biāo)函數(shù)的下降是否充分,若不充分則需要調(diào)整正則化因子,重新計(jì)算目標(biāo)函數(shù)值并將搜索方向置為當(dāng)前的負(fù)梯度方向,重新開(kāi)始迭代直至收斂.算法的詳細(xì)步驟見(jiàn)圖2.

    兩組共軛梯度方向的計(jì)算是上述算法的一個(gè)關(guān)鍵點(diǎn).將(6)—(8)式代入(5)式可得到目標(biāo)函數(shù)的完整表達(dá)形式,將其對(duì)電阻率模型m及畸變參數(shù)c分別求偏導(dǎo),可以得到當(dāng)前的兩組梯度方向:

    (9)

    式中Jm為受到電流型畸變c影響的預(yù)測(cè)阻抗張量dpred關(guān)于電阻率模型m的敏感度矩陣(雅克比矩陣),表達(dá)式如下

    (10)

    式中N=Nsite×Nfreq為數(shù)據(jù)總量(Nsite為測(cè)點(diǎn)數(shù),Nfreq為頻點(diǎn)數(shù)),M為電阻率模型m的電阻率參數(shù)總量.對(duì)于第i個(gè)測(cè)點(diǎn)的第j個(gè)頻點(diǎn)而言,

    (11)

    注意第i個(gè)測(cè)點(diǎn)所受電流型畸變ci與頻率無(wú)關(guān),同等作用于該測(cè)點(diǎn)的所有頻點(diǎn).由于引起電流型畸變的淺層異常體的規(guī)模遠(yuǎn)小于MT信號(hào)波長(zhǎng),MT數(shù)據(jù)的橫、縱向分辨率不足以分辨這些異常體,因此一般認(rèn)為電阻率模型m與電流型畸變c無(wú)關(guān),可得

    圖2 自適應(yīng)壓制電流型畸變的三維MT反演算法流程圖Fig.2 Diagram of the 3D MT inversion algorithm with adaptive galvanic distortion suppression

    (12)

    與常規(guī)的NLCG方法相比,本算法計(jì)算Jm僅需要將各測(cè)點(diǎn)的預(yù)測(cè)畸變參數(shù)與常規(guī)的敏感度矩陣相乘即可.對(duì)于如何高效計(jì)算和存儲(chǔ)常規(guī)的敏感度矩陣,前人已開(kāi)展了大量有意義的工作(Mackie and Madden, 1993; Farquharson and Oldenburg, 1996; 胡祖志等, 2006; Egbert and Kelbert, 2012).本文采用了Egbert和Kelbert (2012)給出的電磁反演問(wèn)題中敏感度矩陣的一般表達(dá)形式,為了縮小計(jì)算量和存儲(chǔ)的規(guī)模,將敏感度矩陣分解為多個(gè)單元分別予以求解.

    (9)式中的Jc是電流型畸變的敏感度矩陣,同樣對(duì)于第i個(gè)測(cè)點(diǎn)的第j個(gè)頻點(diǎn)而言,有

    (13)

    (14)

    3.3自適應(yīng)雙正則化因子策略

    正則化反演的一個(gè)難點(diǎn)是如何確定合理的正則化因子.(5)式中的λ1若取值過(guò)大將會(huì)對(duì)模型約束過(guò)于嚴(yán)厲,難以擬合觀測(cè)數(shù)據(jù),而取值過(guò)小又容易過(guò)度擬合,導(dǎo)致模型的穩(wěn)定性欠佳.本文在引入電流型畸變參數(shù)約束項(xiàng)的同時(shí)引入了一個(gè)新的正則化因子λ2,在反演進(jìn)程中,需要不斷調(diào)整這兩個(gè)正則化因子,保持?jǐn)?shù)據(jù)擬合項(xiàng)(Φd)、電阻率模型約束項(xiàng)(Φm)及畸變參數(shù)約束項(xiàng)(Φc)三者的平衡.

    正則化因子通??刹捎脽o(wú)偏風(fēng)險(xiǎn)預(yù)測(cè)法、廣義交叉檢驗(yàn)法、L-曲線法等方法確定,但這些方法或需要已知數(shù)據(jù)的誤差水平,或需要進(jìn)行大量的實(shí)驗(yàn)性計(jì)算,在三維反演中難以實(shí)現(xiàn).近年來(lái),一種自適應(yīng)正則化因子的方法在MT反演中得到了廣泛應(yīng)用(Zhdanov, 2002; 陳小斌等, 2005; 吳小平和徐果明, 2000).我們將該方法改進(jìn)后拓展到雙正則化因子情況.正則化因子λ1、λ2的初始值仍需在反演前根據(jù)具體情況預(yù)設(shè),通??稍O(shè)置為較大的值以保障在反演前期模型的穩(wěn)定.在反演迭代過(guò)程中,當(dāng)目標(biāo)函數(shù)下降到小于某一閾值時(shí)(通常設(shè)置為一極小值),將正則化因子乘以某一衰減因子獲得到新的正則化因子,重新計(jì)算目標(biāo)函數(shù)值和梯度方向,并將搜索方向置為當(dāng)前的負(fù)梯度方向,以此放寬對(duì)模型的約束,進(jìn)一步擬合觀測(cè)數(shù)據(jù).衰減因子的取值并沒(méi)有明確的標(biāo)準(zhǔn),Zhdanov(2002)認(rèn)為衰減因子的經(jīng)驗(yàn)取值區(qū)間為[0.5,0.9].此外,我們?cè)谶@個(gè)過(guò)程中引入了一組更新效率因子δm和δc:

    (15)

    式中δm和δc分別為電阻率模型和畸變參數(shù)的更新效率因子.當(dāng)算法判斷需要衰減正則化因子進(jìn)一步擬合觀測(cè)數(shù)據(jù)時(shí),通過(guò)比較當(dāng)前的兩組更新效率因子自適應(yīng)地判斷是衰減λ1還是衰減λ2.如果電阻率模型更新效率優(yōu)于畸變參數(shù)(即δm>δc),衰減λ2以放寬對(duì)畸變參數(shù)的約束,反之衰減λ1以放寬對(duì)模型參數(shù)的約束.通過(guò)采用上述自適應(yīng)調(diào)整正則化因子策略,算法既能滿足對(duì)觀測(cè)數(shù)據(jù)的高效擬合,又能保障電阻率模型和畸變參數(shù)的穩(wěn)定.

    4三維理論模型測(cè)試

    4.1模型設(shè)置與響應(yīng)數(shù)據(jù)

    圖3所示為一個(gè)理論三維電阻率模型(mtrue).在100 Ωm的背景半空間中分布著三個(gè)電阻率分別為10 Ωm、500 Ωm和5000 Ωm的異常塊體.采用交錯(cuò)網(wǎng)格有限差分算法(Staggered-Grid Finite Difference method)作正演計(jì)算(Mackie and Madden, 1993),共得到225個(gè)測(cè)點(diǎn)的全阻抗張量數(shù)據(jù)(周期范圍為1~1000 s,周期數(shù)為8).該數(shù)據(jù)不含電流畸變成分,下文中稱為正常數(shù)據(jù).由于理論模型測(cè)試的主要目的是探討新方法在壓制電流型畸變上的有效性,所以該數(shù)據(jù)中未加入任何人為噪聲.

    4.2反演結(jié)果對(duì)比和分析

    4.2.1電阻率模型比較分析

    為了便于比較分析,所有反演方案均采用如下設(shè)置:

    圖3 三維理論模型示意圖(a) 平面視圖; (b) 沿測(cè)線P1的剖面視圖; 白色小圓點(diǎn)為225個(gè)MT測(cè)點(diǎn); 黑色虛線為選取的穿過(guò)模型的測(cè)線P1,白色大圓點(diǎn)S1為選取的一個(gè)代表性測(cè)點(diǎn).Fig.3 Sketch map of the 3-D synthetic model(a) Plan view; (b) Cross-section along the profile P1. Small circles are the location of 225 MT sites; the black dashed line represents a selected profile P1, the bigger circle S1 marks a selected site.

    (1) 正常數(shù)據(jù):由4.1節(jié)所示模型的正演計(jì)算所得,不含電流畸變成分和任何噪音;

    (2) 畸變數(shù)據(jù):將隨機(jī)生成的理論畸變參數(shù)(ctrue)作用于正常數(shù)據(jù)產(chǎn)生的受畸變的數(shù)據(jù)(稱為畸變數(shù)據(jù)).

    兩類反演方法:

    (1) 不考慮電流型畸變因素的常規(guī)三維反演方法(Kelbert et al., 2014)(稱為常規(guī)方法);

    (2) 考慮電流型畸變因素的三維反演算法(本文方法).

    測(cè)試中共設(shè)計(jì)了4種測(cè)試方案,方案1和方案2分別為常規(guī)方法和本文方法對(duì)正常數(shù)據(jù)的反演,目的是檢測(cè)在無(wú)畸變條件下(即均采用正常數(shù)據(jù))本文方法與常規(guī)方法所得結(jié)果的一致性.常規(guī)方法的可靠性已經(jīng)過(guò)了大量的測(cè)試,若方案1和方案2的結(jié)果相同或接近,則表明本文方法在面對(duì)正常數(shù)據(jù)時(shí)是可靠的.方案3和方案4分別為常規(guī)方法和本文方法對(duì)畸變數(shù)據(jù)的反演,目的是表明在面對(duì)畸變數(shù)據(jù)時(shí)的常規(guī)方法的局限性以及本文方法的有效性.四種測(cè)試方案的具體情況見(jiàn)表1.

    方案1:常規(guī)方法反演正常數(shù)據(jù)

    方案1的目的是檢測(cè)常規(guī)方法在數(shù)據(jù)未受到電流型畸變干擾情況下的反演效果.將正常數(shù)據(jù)和先驗(yàn)電阻率模型(mprior)的響應(yīng)數(shù)據(jù)代入(6)式,可以得到數(shù)據(jù)擬合項(xiàng)(Φd)的初值為6.589.經(jīng)過(guò)17次迭代,隨著預(yù)測(cè)電阻率模型逐漸接近理論模型,Φd的值由6.589減小為1.003(圖4a中的黑色實(shí)線),達(dá)到理想擬合值1(圖4a中的黑色虛線),滿足擬合條件,迭代正常終止.電阻率模型約束項(xiàng)(Φm)的值反映了電阻率模型的粗糙度,它的初值為0(初始和先驗(yàn)電阻率模型相同,均為均勻半空間),在迭代過(guò)程中隨著電阻率模型逐漸偏離先驗(yàn)電阻率模型,Φm增大至0.455(圖4b中黑色實(shí)線),趨于其理論值2.761(圖4b中黑色虛線).圖5b、6b分別展示了方案1反演所得電阻率模型沿P1測(cè)線的剖面圖及不同深度的切片圖,與理論模型(圖5a、6a)對(duì)比,除了三個(gè)異常體底部由于MT方法的特性無(wú)法得到精確約束外,電阻率大小和尺度都得到了較好的還原,尤其是各塊體之間的橫向電性分界面清晰可辨,表明常規(guī)方法對(duì)正常數(shù)據(jù)的反演是正確有效的.

    表1 四種反演方案下目標(biāo)函數(shù)中各項(xiàng)的變化情況

    方案2:本文方法反演正常數(shù)據(jù)

    采用本文方法同時(shí)反演電阻率模型參數(shù)與電流型畸變參數(shù)時(shí),所面臨的一種風(fēng)險(xiǎn)是可能會(huì)將一些觀測(cè)數(shù)據(jù)中包含地下電性結(jié)構(gòu)的信息誤判為畸變作用,從而導(dǎo)致真實(shí)的電性結(jié)構(gòu)被掩蓋或者扭曲,這相當(dāng)于方法本身帶入了某種不確定的噪聲.因此,有必要驗(yàn)證當(dāng)數(shù)據(jù)未受畸變干擾時(shí)本文方法是否能夠同樣穩(wěn)定地還原真實(shí)電性結(jié)構(gòu).將λ2的初值設(shè)為600,初始畸變張量(c0)和先驗(yàn)畸變張量(cprior)均設(shè)為單位矩陣I,其余參數(shù)與方案1保持一致.經(jīng)過(guò)30次迭代后,該項(xiàng)由6.589減小為1.036(圖4a中藍(lán)色實(shí)線),同樣達(dá)到預(yù)期目標(biāo),反演終止.電阻率模型約束項(xiàng)(Φm)的終值為0.508(圖4b中藍(lán)色實(shí)線),接近方案1結(jié)果.由于數(shù)據(jù)未受到任何畸變影響,整個(gè)反演過(guò)程中畸變參數(shù)約束項(xiàng)(Φc)保持為一接近于0的極小值(圖4c中的藍(lán)色實(shí)線),符合實(shí)際.電阻率模型反演結(jié)果如圖5c、6c所示,與方案1反演結(jié)果(圖5b、6b)基本一致,三個(gè)異常體均得以成功還原,說(shuō)明在沒(méi)有畸變時(shí)本文方法與常規(guī)方法是一致的,沒(méi)有帶入附加噪聲,也沒(méi)有丟失有效信息.

    圖4 四種反演方案下Φd,Φm和Φc三項(xiàng)隨迭代次數(shù)的變化情況(a) 數(shù)據(jù)擬合項(xiàng)(Φd); (b) 電阻率模型約束項(xiàng)(Φm);(c) 電流型畸變參數(shù)約束項(xiàng)(Φc).黑色虛線為各項(xiàng)的理想值.Fig.4 Variation of Φd, Φm and Φc versus iteration number for the four different inversion strategies(a) Data misfit term (Φd); (b) Resistivity model constraint term (Φm); (c) Galvanic distortion constraint term (Φc). The black dashed lines represent the theoretical values of these terms.

    圖5 比較四種反演方案所得模型沿測(cè)線P1的電阻率剖面圖(a) 理論模型; (b) 方案1所得反演結(jié)果,較好還原了三個(gè)異常體; (c) 方案2所得反演結(jié)果,與方案1基本一致; (d) 方案3所得反演結(jié)果,嚴(yán)重偏離理論模型(如圖中的虛假低阻異常A); (e) 方案4所得反演結(jié)果, 與方案3相比更接近理論模型, 畸變得到有效壓制.Fig.5 Comparison of resistivity cross-sections along the profile P1 derived from the four inversion strategies(a) Synthetic model; (b) Inversion result derived from strategy 1 (conventional inversion of the undistorted data); (c) Inversion result derived from strategy 2 (new algorithm inversion of the undistorted data); (d) Inversion result derived from strategy 3 (conventional inversion of the distorted data); (e) Inversion result derived from strategy 4 (new algorithm inversion of the distorted data).

    方案3:常規(guī)方法反演畸變數(shù)據(jù)

    將畸變數(shù)據(jù)和先驗(yàn)電阻率模型(100 Ωm的均勻半空間)的響應(yīng)數(shù)據(jù)代人(6)式,得到數(shù)據(jù)擬合項(xiàng)(Φd)的初值為10.061.其他參數(shù)與方案1相同.經(jīng)過(guò)84次迭代后陷入局部最優(yōu),反演被強(qiáng)行終止.數(shù)據(jù)擬合項(xiàng)(Φd)最終減小至6.475(圖4a中紅色實(shí)線),遠(yuǎn)大于理想值1(圖4a中黑色虛線),數(shù)據(jù)擬合質(zhì)量差;電阻率模型約束項(xiàng)(Φm)的終值為4.651(圖4b中紅色實(shí)線),遠(yuǎn)大于理論值2.761(圖4b中黑色虛線),表明反演過(guò)程中常規(guī)算法對(duì)電阻率模型的約束基本失效.同時(shí),三個(gè)異常塊體的形態(tài)和電阻率幅值均遭受到嚴(yán)重的畸變,深部背景空間中出現(xiàn)了大范圍的虛假低阻區(qū)域(如圖5d中的A),并且模型淺層出現(xiàn)了大量虛假的小異常體(圖6d).這說(shuō)明,當(dāng)觀測(cè)數(shù)據(jù)含有電流型畸變時(shí),常規(guī)方法的反演結(jié)果偏離真實(shí)結(jié)構(gòu).

    圖6 比較四種反演方案所得電阻率模型在不同深度的水平切片圖(a)—(e)與圖5相同.Fig.6 Comparison of horizontal resistivity slices at different depths derived from the four inversion strategies(a)—(e) same as in Fig.5.

    圖7 比較所有測(cè)點(diǎn)受到的電流畸變理論值與方案4所得預(yù)測(cè)值(a) 電流型畸變參數(shù)理論值(ctrue); (b) 本文方法預(yù)測(cè)的電流型畸變參數(shù)(cpred).Fig.7 Comparison of theoretical and predicted galvanic distortion values(a) Theoretical galvanic distortion values (ctrue); (b) Galvanic distortion values predicted by our new inversion algorithm (cpred).

    方案4:本文方法反演畸變數(shù)據(jù)

    反演參數(shù)與方案3完全相同,采用本文方法經(jīng)過(guò)77次迭代后,數(shù)據(jù)擬合項(xiàng)(Φd)由10.061降低為2.210(圖4a中綠色實(shí)線),雖然未達(dá)到理論目標(biāo)值1,但相較于方案3所得最終值(6.475),數(shù)據(jù)的擬合程度得到了顯著改善;電阻率模型約束項(xiàng)(Φm)(圖4b中綠色實(shí)線)最終值為0.855,畸變參數(shù)模型約束項(xiàng)(Φc)(圖4c中綠色實(shí)線)的最終值為8.635×10-2,均趨近于理論值,表明電阻率模型和電流型畸變參數(shù)都得到了有效地約束.從反演的電阻率模型來(lái)看,盡管數(shù)據(jù)受到了畸變的干擾,本文方法仍較好地還原了理論模型中的三個(gè)異常體.與常規(guī)方法(方案3)的結(jié)果(圖5d, 圖6d)比較,本文方法的結(jié)果(圖5e,圖6e)更接近理論模型,表明反演質(zhì)量得到了明顯改善.

    4.2.2電流型畸變參數(shù)比較分析

    圖7比較了所有測(cè)點(diǎn)所遭受的理論電流型畸變參數(shù)(ctrue)與本文方法預(yù)測(cè)的電流型畸變參數(shù)(cpred),二者整體上具有較好的一致性,說(shuō)明本文方法有效地還原了電流型畸變參數(shù).這一點(diǎn)對(duì)于評(píng)價(jià)本文方法的有效性相當(dāng)重要,在反演過(guò)程中,電流型畸變成分如果未能得到有效的識(shí)別和還原,勢(shì)必將參與到電阻率模型的還原中,從而導(dǎo)致電阻率模型的畸變.

    4.3數(shù)據(jù)擬合分析

    僅根據(jù)數(shù)據(jù)擬合項(xiàng)(Φd)的取值情況并不能全面、真實(shí)地反映數(shù)據(jù)的擬合程度.圖8、圖9分別給出了所有測(cè)點(diǎn)阻抗張量的非對(duì)角和對(duì)角分量在周期1.389 s處的視電阻率和相位平面圖.對(duì)比正常數(shù)據(jù)(圖8a,圖9a)和畸變數(shù)據(jù)(圖8b,圖9b),可以看出當(dāng)三維區(qū)域阻抗張量受到畸變時(shí),對(duì)角和非對(duì)角分量的視電阻率和相位均受到不同程度的影響,特別是對(duì)角分量受到的影響更嚴(yán)重.比較圖8b和圖8a, 非對(duì)角分量的視電阻率受到畸變的干擾比較明顯,而相位受到的影響較小,這說(shuō)明即使在三維結(jié)構(gòu)中,非對(duì)角分量仍主要表現(xiàn)為靜態(tài)效應(yīng)特征(即幅值變化大,而相位不變或變化很小).比較圖9b和圖9a, 對(duì)角分量的模值和相位同時(shí)受到嚴(yán)重的畸變干擾,這是由于對(duì)角分量的模值較小,更容易受到電流型畸變的干擾.采用常規(guī)算法反演畸變數(shù)據(jù)時(shí)(方案3),數(shù)據(jù)的擬合程度很差(圖8c對(duì)比圖8a;圖9c對(duì)比圖9a).面對(duì)同樣的畸變數(shù)據(jù),本文方法顯著改善了數(shù)據(jù)的擬合質(zhì)量(對(duì)比圖8d和圖8a).即使對(duì)于受畸變影響嚴(yán)重的對(duì)角分量,本文方法的擬合質(zhì)量也比常規(guī)方法好很多(對(duì)比圖9d和圖9a).

    為了分析數(shù)據(jù)各頻點(diǎn)的擬合情況,圖10給出了S1點(diǎn)的全阻抗張量數(shù)據(jù)隨周期的變化曲線.從圖中可以看出,在計(jì)算的頻段范圍內(nèi),與常規(guī)方法(方案3)(圖10中的紅色圓點(diǎn))比較,本文方法(方案4)計(jì)算的阻抗張量(圖10中的綠色三角)更接近理論值(圖10中的黑色實(shí)線).進(jìn)一步表明在全頻段范圍內(nèi),本文方法的擬合質(zhì)量明顯優(yōu)于常規(guī)方法.

    圖8 比較周期1.389 s處方案3和方案4預(yù)測(cè)模型響應(yīng)的非對(duì)角分量擬合情況(a) 正常數(shù)據(jù); (b) 畸變數(shù)據(jù); (c) 方案3預(yù)測(cè)的模型響應(yīng); (d) 方案4預(yù)測(cè)的模型響應(yīng).Fig.8 Comparison of the off-diagonal data fits for strategy 3 and 4 at period 1.389 s(a) Synthetic data; (b) Distorted data; (c) Predicted responses by strategy 3; (d) Predicted responses by strategy 4.

    5結(jié)論與討論

    在不考慮外界人為噪聲的情況下,對(duì)大地電磁觀測(cè)信號(hào)影響最嚴(yán)重的地質(zhì)噪音當(dāng)屬電流型畸變(包含靜位移效應(yīng)).數(shù)學(xué)上來(lái)說(shuō),解決該問(wèn)題的關(guān)鍵在于求解一個(gè)非線性的欠定問(wèn)題.當(dāng)區(qū)域結(jié)構(gòu)可簡(jiǎn)化為一、二維時(shí),通過(guò)阻抗張量分解等方法可以在反演前將部分畸變成分分解出來(lái),在一定程度上緩解畸變帶來(lái)的危害.然而,當(dāng)區(qū)域結(jié)構(gòu)具有較強(qiáng)的三維性時(shí),目前尚無(wú)有效措施應(yīng)對(duì)該畸變.本文試圖發(fā)展一種考慮電流型畸變的三維反演算法,在反演過(guò)程中自適應(yīng)地同步擬合電阻率模型和電流型畸變參數(shù),從而達(dá)到壓制畸變的影響、提升電阻率模型的可信度的目的.

    與常規(guī)反演電阻率模型一樣,在反演電流型畸變參數(shù)時(shí)同樣面臨非唯一性問(wèn)題.如何在反演過(guò)程中對(duì)其進(jìn)行有效約束是解決該問(wèn)題的關(guān)鍵.基于統(tǒng)計(jì)分析大量實(shí)測(cè)MT數(shù)據(jù)所遭受的畸變分布特征,我們認(rèn)為,在無(wú)法通過(guò)其他手段獲得更精確的先驗(yàn)畸變信息的情況下,采用單位矩陣作為先驗(yàn)畸變參數(shù)并且根據(jù)預(yù)測(cè)畸變參數(shù)和先驗(yàn)畸變參數(shù)之差的L2范數(shù)來(lái)衡量畸變模型的復(fù)雜程度(粗糙度)是合理、有效的.

    為了保障反演算法的穩(wěn)定和效率,我們將二維反演中常用的自適應(yīng)正則化因子策略改進(jìn)后推廣到雙正則化因子情況中,并且通過(guò)引入新的更新效率因子,依據(jù)電阻率模型和畸變兩組參數(shù)的收斂情況,在兩個(gè)正則化因子之間自主選擇需要衰減的因子.

    圖9 比較周期1.389 s處方案3和方案4預(yù)測(cè)模型響應(yīng)的對(duì)角分量擬合情況(a)—(d)與圖8相同.Fig.9 Comparison of the diagonal data fits for strategy 3 and 4 at period 1.389 s.(a)—(d) same as in Fig.8.

    圖10 比較測(cè)點(diǎn)S1處模型理論響應(yīng)和方案3、4所得預(yù)測(cè)響應(yīng)曲線黑色實(shí)線為理論模型響應(yīng),紅色圓點(diǎn)為方案3預(yù)測(cè)模型響應(yīng),綠色正方形為方案4預(yù)測(cè)模型響應(yīng).Fig.10 Comparison of the theoretical and predicted response curves for strategy 3 and 4 at site S1.Solid black line: synthetic model responses; red circle: predicted responses by strategy 3;green square: predicted responses by strategy 4.

    在反演初期,偏重于電阻率模型和畸變兩組參數(shù)的穩(wěn)定性;在反演后期,當(dāng)兩組參數(shù)趨于穩(wěn)健后,適度放寬約束,進(jìn)一步擬合觀測(cè)數(shù)據(jù).需要注意的是兩個(gè)正則化因子的初值仍需要在反演前根據(jù)具體情況設(shè)定并會(huì)對(duì)反演結(jié)果產(chǎn)生一定影響.

    理論模型測(cè)試表明,對(duì)于受到電流型畸變的觀測(cè)數(shù)據(jù),常規(guī)的三維反演方法難以還原真實(shí)的電性結(jié)構(gòu).面對(duì)同等強(qiáng)度的畸變干擾,采用本文方法能夠同時(shí)還原真實(shí)的電性結(jié)構(gòu)和電流型畸變參數(shù).最后需要指出的是,雖然理論模型測(cè)試驗(yàn)證了本文方法在壓制電流型畸變干擾上的有效性,但地球深部的實(shí)際情況要遠(yuǎn)比理論模型復(fù)雜,面對(duì)受到各類噪聲干擾的實(shí)測(cè)數(shù)據(jù),該算法要取得理想的效果仍需要進(jìn)一步的研究和改進(jìn).

    致謝文中所用MT實(shí)測(cè)數(shù)據(jù)來(lái)自USARRAY項(xiàng)目.感謝俄勒岡大學(xué)Gary Egbert教授和Anna Kelbert博士提供MT三維正反演代碼和熱心幫助.匿名審稿專家提出的寶貴意見(jiàn)和建議對(duì)于本文的改進(jìn)和提高非常重要.

    ReferencesAvdeeva A, Moorkamp M, Avdeev D, et al. 2015. Three-dimensional inversion of magnetotelluric impedance tensor data and full distortion matrix.GeophysicalJournalInternational, 202(1): 464-481. Bahr K. 1988. Interpretation of the magnetotelluric impedance tensor: regional induction and local telluric distortion.JournalofGeophysics, 62(2): 119-127. Berdichevsky M N, Dmitriev V I. 1976. Distortion of magnetic and electrical fields by near-surface lateral inhomogeneities.ActaGeodaetica,GeophysicaetMontanistica, 11(3-4): 447-483.

    Bibby H M, Caldwell T G, Brown C. 2005. Determinable and non-determinable parameters of galvanic distortion in magnetotellurics.GeophysicalJournalInternational, 163(3): 915-930. Caldwell T G, Bibby H M, Brown C. 2004. The magnetotelluric phase tensor.GeophysicalJournalInternational, 158(2): 457-469.

    Chave A D, Smith J T. 1994. On electric and magnetic galvanic distortion tensor decompositions.JournalofGeophysicalResearch:SolidEarth, 99(B3): 4669-4682.Chen X B, Zhao G Z, Tang J, et al. 2005. An adaptive regularized inversion algorithm for magnetotelluric data.ChineseJournalofGeophysics(in Chinese), 48(4): 937-946, doi: 10.3321/j.issn:0001-5733.2005.04.029.

    DeGroot-Hedlin C. 1995. Inversion for regional 2-D resistivity structure in the presence of galvanic scatterers.GeophysicalJournalInternational, 122(3): 877-888. Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems.GeophysicalJournalInternational, 189(1): 251-267. Farquharson C G, Oldenburg D W. 1996. Approximate sensitivities for the electromagnetic inverse problem.GeophysicalJournalInternational, 126(1): 235-252.

    Fletcher R, Reeves C M. 1964. Function minimization by conjugate gradients.TheComputerJournal, 7(2): 149-154.

    Garcia X, Jones A G. 2002. Decomposition of three-dimensional magnetotelluric data. ∥Proceedings of the Second International Symposium.Amsterdam:Elsevier, 35: 235-250.

    Groom R W, Bailey R C. 1989. Decomposition of magnetotelluric impedance tensors in the presence of local three-dimensional galvanic distortion.JournalofGeophysicalResearch:SolidEarthandPlanets, 94(B2): 1913-1925.

    Groom R W, Bailey R C. 1991. Analytic investigations of the effects of near-surface three-dimensional galvanic scatterers on MT tensor decompositions.Geophysics, 56(4): 496-518.

    Groom R W, Bahr K. 1992. Corrections for nearsurface effects: Decomposition of the magnetotelluric impedance tensor and scaling corrections for regional resistivities: Atutorial.SurveysinGeophysics, 13(4-5): 341-379.

    Hestenes M R.1951. Iterative methods for solving linear equations. NAML Report 52-9.Los Angeles, CA: National Bureau of Standards.

    Hestenes M R, Stiefel E. 1952. Methods of conjugate gradients for solving linear systems.JournalofResearchoftheNationalBureauofStandards, 49(6): 409-436.

    Hu Z Z, Hu X Y, He Z X. 2006. Pseudo-Three-Dimensional magnetotelluric inversion using nonlinear conjugate gradients.ChineseJournalofGeophysics(in Chinese), 49(4): 1226-1234, doi: 10.3321/j.issn:0001-5733.2006.04.039.

    Jin G W, Kong X R.2006. The Distortion and Decomposition of Magnetotelluric Impedance Tensor (in Chinese). Beijing: Seismological Press.

    Jiracek G R. 1990. Near-surface and topographic distortions in electromagnetic induction.SurveysinGeophysics, 11(2-3): 163-203.

    Jones A G. 1988. Static shift of magnetotelluric data and its removal in a sedimentary basin environment.Geophysics, 53(7): 967-978.

    Jones A G. 2011. Three-dimensional galvanic distortion of three-dimensional regional conductivity structures: Comment on “Three-dimensional joint inversion for magnetotelluric resistivity and static shift distributions in complex media” by Yutaka Sasaki and Max A. Meju.JournalofGeophysicalResearch:SolidEarth, 116(B12), doi: 10.1029/2011JB008665.

    Kelbert A, Meqbel N, Egbert G D, et al. 2014. ModEM: A modular system for inversion of electromagnetic geophysical data.Computers&Geosciences, 66: 40-53.

    Mackie R L, Madden T R. 1993. Three-dimensional magnetotelluric inversion using conjugate gradients.GeophysicalJournalInternational, 115(1): 215-229. McNeice G W, Jones A G. 2001. Multisite, multifrequency tensor decomposition of magnetotelluric data.Geophysics, 66(1): 158-173. Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients.GeophysicalJournalInternational, 140(2): 410-424. Patro P K, Uyeshima M, Siripunvaraporn W. 2013. Three-dimensional inversion of magnetotelluric phase tensor data.GeophysicalJournalInternational, 192(1): 58-66. Pellerin L, Hohmann G W. 1990. Transient electromagnetic inversion: Aremedy for magnetotelluric static shifts.Geophysics, 55(9): 1242-1250. Polak E, Ribière G. 1969. Note sur la convergence de méthodes de directions conjuguées. Revue Fran?aise d′Informatique et de Recherche Opérationnelle, 16: 35-43.

    Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion.Geophysics, 66(1): 174-187.

    Sasaki Y, Meju M A. 2006. Three-dimensional joint inversion for magnetotelluric resistivity and static shift distributions in complex media.JournalofGeophysicalResearch:SolidEarth, 111(B5): B05101.Simpson F, Bahr K. 2005. Practical Magnetotellurics. Cambridge:Cambridge University Press.

    Sternberg B K, Washburne J C, Pellerin L. 1988. Correction for the static shift in magnetotellurics using transient electromagnetic soundings.Geophysics, 53(11): 1459-1468.

    Stiefel E.1951. Some special methods of relaxation technique.∥Peoceedings of the symposium on simultaneous linear equations and the determination of eigenvalues.Washington, D.C.:U.S. Government Printing Office.

    Swift C M Jr.1967. A magnetotelluric investigation of an electrical conductivity anomaly in the Southwestern United States[Ph. D. thesis].Massachusetts:Massachusetts Institute of Technology.

    Utada H, Munekane H. 2000. On galvanic distortion of regional

    three-dimensional magnetotelluric impedances.GeophysicalJournalInternational, 140(2): 385-398. Wang J Y. 1998. Inverse Theory in Geophysics. Wuhan: China University of Geosciences press.

    Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method.ChineseJournalofGeophysics(in Chinese), 43(3): 420-427, doi: 10.3321/j.issn:0001-5733.2000.03.016.

    Xiao Q B, Shao G H, Jing L Z, et al. 2015. Eastern termination of the altyn tagh fault, western China: Constraints from a magnetotelluric survey.JournalofGeophysicalResearch:SolidEarth, 120(5): 2838-2858.

    Zhang L T, Unsworth M, Jin S, et al. 2015. Structure of the Central Altyn Tagh Fault revealed by magnetotelluric data: New insights into the structure of the northern margin of the India-Asia collision.EarthandPlanetaryScienceLetters, 415: 67-79.

    Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. Amsterdam: Elsevier Science.

    附中文參考文獻(xiàn)

    陳小斌, 趙國(guó)澤, 湯吉等. 2005. 大地電磁自適應(yīng)正則化反演算法. 地球物理學(xué)報(bào),48(4): 937-946, doi: 10.3321/j.issn:0001-5733.2005.04.029.

    胡祖志, 胡祥云, 何展翔. 2006. 大地電磁非線性共軛梯度擬三維反演. 地球物理學(xué)報(bào),49(4): 1226-1234, doi: 10.3321/j.issn:0001-5733.2006.04.039.

    晉光文, 孔祥儒. 2006. 大地電磁阻抗張量的畸變與分解. 北京: 地震出版社.

    王家映. 1998. 地球物理反演理論. 北京: 中國(guó)地質(zhì)大學(xué)出版社.

    吳小平, 徐果明. 2000. 利用共軛梯度法的電阻率三維反演研究. 地球物理學(xué)報(bào), 43(3): 420-427, doi: 10.3321/j.issn:0001-5733.2000.03.016.

    (本文編輯胡素芳)

    基金項(xiàng)目國(guó)家自然科學(xué)基金(41330212,41274102,41474056)資助.

    作者簡(jiǎn)介李鑫,男,1986年生,在讀博士研究生,主要研究方向?yàn)榇蟮仉姶?MT)三維反演技術(shù)和青藏高原深部電性結(jié)構(gòu)和動(dòng)力學(xué). E-mail:lixin342@mail.iggcas.ac.cn

    doi:10.6038/cjg20160632 中圖分類號(hào)P319

    收稿日期2015-12-17,2016-03-14收修定稿

    Three-dimensional inversion of magnetotelluric resistivity model with galvanic distortion

    LI Xin1,2, BAI Deng-Hai1, YAN Yong-Li1, MA Xiao-Bing1, KONG Xiang-Ru1

    1InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China

    AbstractGalvanic distortion, which is caused by shallow small-scale inhomogeneities or topography, has long been recognized as an obstacle in magnetotelluric (MT) sounding studies. Removal of these distortions from the field-observed MT data is essential to obtain a reliable model of the subsurface electrical structure. In two-dimensional (2-D) situations, a number of impedance decomposition methods have been widely used for distortion removal. For data from three-dimensional (3-D) structures, however, the effects of galvanic distortion become far more complex and cannot be solved by these traditional methods. This paper proposes a method that employs the Nonlinear Conjugate Gradient algorithm (NLCG) to solve both the regional resistivity structures and the parameters of full galvanic distortion during iterations. The statistical distribution of galvanic distortion is firstly estimated by using the phase tensor method and suggests that the off-diagonal elements of the distortion tensor are quasi-symmetrically distributed around zero, and the diagonal elements are centralized around one. Using the prior information contained in the field data, we build a new objective function that takes the galvanic effects into consideration and adds an extra constraint on the distortion parameters. Furthermore, to ensure stability and efficiency of the algorithm, an adaptive regularized strategy is used to determine the trade-off factors between the data misfit term and the roughness terms of the resistivity model and galvanic distortion respectively. The synthetic model study suggests that conventional inversion approach cannot reasonably explain the distorted part of the synthetic data and will inevitably bring erroneous structures to the resistivity model. In contrast, our new inversion algorithm can reliably retrieve both the regional resistivity structures and galvanic distortion, irrespective of whether galvanic distortion is present or not.

    KeywordsGalvanic distortion; Static shift; 3-D inversion; Magnetotelluric

    李鑫, 白登海, 閆永利. 2016. 考慮電流型畸變的大地電磁三維反演. 地球物理學(xué)報(bào),59(6):2302-2315,doi:10.6038/cjg20160632.

    Li X, Bai D H, Yan Y L. 2016. Three-dimensional inversion of magnetotelluric resistivity model with galvanic distortion.ChineseJ.Geophys. (in Chinese),59(6):2302-2315,doi:10.6038/cjg20160632.

    国产免费av片在线观看野外av| 精品国产国语对白av| 精品久久久精品久久久| 亚洲精品中文字幕在线视频| 在线永久观看黄色视频| 丝袜美腿诱惑在线| 欧美丝袜亚洲另类 | 亚洲欧美精品综合一区二区三区| 999久久久国产精品视频| 精品人妻在线不人妻| 啦啦啦 在线观看视频| av超薄肉色丝袜交足视频| 欧美色欧美亚洲另类二区 | 长腿黑丝高跟| 日韩高清综合在线| 成人18禁高潮啪啪吃奶动态图| 国产精品乱码一区二三区的特点 | 成人国语在线视频| 99国产精品99久久久久| 日本免费一区二区三区高清不卡 | 国产1区2区3区精品| 欧美日韩瑟瑟在线播放| 亚洲国产精品sss在线观看| 成人亚洲精品一区在线观看| 国产成人av教育| 亚洲精品国产色婷婷电影| 两性午夜刺激爽爽歪歪视频在线观看 | 少妇 在线观看| 久久久久久免费高清国产稀缺| 精品一品国产午夜福利视频| 久久亚洲真实| 国产精品永久免费网站| www.自偷自拍.com| 精品乱码久久久久久99久播| 午夜精品在线福利| 亚洲欧美日韩高清在线视频| 少妇粗大呻吟视频| 亚洲av日韩精品久久久久久密| 别揉我奶头~嗯~啊~动态视频| 欧美一级a爱片免费观看看 | 亚洲一区中文字幕在线| 天天躁狠狠躁夜夜躁狠狠躁| 色播亚洲综合网| 非洲黑人性xxxx精品又粗又长| 亚洲 欧美一区二区三区| 999久久久精品免费观看国产| 亚洲成人国产一区在线观看| 国产精华一区二区三区| 国产激情欧美一区二区| 久久久久亚洲av毛片大全| 女生性感内裤真人,穿戴方法视频| 午夜福利免费观看在线| 国产欧美日韩综合在线一区二区| 国产精品一区二区免费欧美| 国产精品亚洲一级av第二区| 亚洲av成人av| 亚洲人成伊人成综合网2020| 国产成人精品久久二区二区免费| 国产亚洲精品第一综合不卡| 亚洲欧美一区二区三区黑人| 日本在线视频免费播放| 欧美日韩黄片免| 国产精品二区激情视频| 亚洲色图 男人天堂 中文字幕| 黄片播放在线免费| 日本黄色视频三级网站网址| 久久精品影院6| 嫁个100分男人电影在线观看| 天堂影院成人在线观看| 亚洲国产欧美日韩在线播放| 露出奶头的视频| 久久这里只有精品19| 波多野结衣巨乳人妻| 久久人人97超碰香蕉20202| 久久久久久久午夜电影| 老司机午夜福利在线观看视频| 老司机深夜福利视频在线观看| 国产三级在线视频| av在线天堂中文字幕| 亚洲成a人片在线一区二区| 国产熟女xx| 亚洲国产精品999在线| 欧美成人性av电影在线观看| 色尼玛亚洲综合影院| 亚洲欧美日韩另类电影网站| 亚洲第一青青草原| 国产国语露脸激情在线看| 9191精品国产免费久久| 一区福利在线观看| 亚洲色图 男人天堂 中文字幕| 日日夜夜操网爽| 亚洲专区中文字幕在线| 极品教师在线免费播放| 国产单亲对白刺激| 欧美日韩瑟瑟在线播放| 国产免费av片在线观看野外av| 日韩高清综合在线| 精品卡一卡二卡四卡免费| 涩涩av久久男人的天堂| 中亚洲国语对白在线视频| 伦理电影免费视频| 久久精品91蜜桃| 老汉色av国产亚洲站长工具| 桃色一区二区三区在线观看| 欧美 亚洲 国产 日韩一| 国产一区在线观看成人免费| 欧美乱妇无乱码| 亚洲 欧美一区二区三区| 神马国产精品三级电影在线观看 | 十八禁网站免费在线| 久久久久九九精品影院| 黄片播放在线免费| 在线视频色国产色| 青草久久国产| 亚洲 欧美一区二区三区| 婷婷精品国产亚洲av在线| 制服丝袜大香蕉在线| 熟女少妇亚洲综合色aaa.| 麻豆成人av在线观看| 亚洲,欧美精品.| 99香蕉大伊视频| 午夜日韩欧美国产| 日韩欧美一区视频在线观看| 男人的好看免费观看在线视频 | 欧洲精品卡2卡3卡4卡5卡区| 黄色 视频免费看| 琪琪午夜伦伦电影理论片6080| 两人在一起打扑克的视频| 色婷婷久久久亚洲欧美| 在线免费观看的www视频| 在线观看日韩欧美| 好看av亚洲va欧美ⅴa在| 国产又色又爽无遮挡免费看| 一级作爱视频免费观看| 一本大道久久a久久精品| 一区福利在线观看| 欧美xxxx性猛交bbbb| 国产伦在线观看视频一区| 亚洲欧美日韩无卡精品| 亚洲精品久久国产高清桃花| 日韩在线高清观看一区二区三区 | 一级黄片播放器| 成人永久免费在线观看视频| 人妻夜夜爽99麻豆av| 人妻夜夜爽99麻豆av| 日韩av在线大香蕉| 国产精品嫩草影院av在线观看 | 亚洲成a人片在线一区二区| 国产在线精品亚洲第一网站| 亚洲精品一卡2卡三卡4卡5卡| 亚洲精品在线观看二区| 3wmmmm亚洲av在线观看| 丰满的人妻完整版| 两个人视频免费观看高清| 国产淫片久久久久久久久| 欧美日韩亚洲国产一区二区在线观看| 久久久久久久亚洲中文字幕| netflix在线观看网站| 搡老熟女国产l中国老女人| 国产男人的电影天堂91| 琪琪午夜伦伦电影理论片6080| 1024手机看黄色片| 亚洲av不卡在线观看| 十八禁国产超污无遮挡网站| 在线国产一区二区在线| 成人午夜高清在线视频| 中文字幕精品亚洲无线码一区| 亚洲久久久久久中文字幕| 国产精品亚洲美女久久久| 日本五十路高清| 噜噜噜噜噜久久久久久91| 免费在线观看成人毛片| 国产精品野战在线观看| 国内精品一区二区在线观看| 亚洲精品在线观看二区| 九九热线精品视视频播放| 日本一二三区视频观看| 成人亚洲精品av一区二区| 国产精品乱码一区二三区的特点| 97人妻精品一区二区三区麻豆| 偷拍熟女少妇极品色| 午夜爱爱视频在线播放| 亚洲中文日韩欧美视频| 村上凉子中文字幕在线| av在线天堂中文字幕| 真人一进一出gif抽搐免费| 亚洲成人久久爱视频| 久久久精品大字幕| 国产老妇女一区| 91av网一区二区| 校园人妻丝袜中文字幕| 国内久久婷婷六月综合欲色啪| 午夜福利高清视频| 国产成人影院久久av| 成人精品一区二区免费| avwww免费| 国产欧美日韩一区二区精品| 两性午夜刺激爽爽歪歪视频在线观看| 国产黄a三级三级三级人| 久久精品国产99精品国产亚洲性色| 欧美区成人在线视频| 亚洲av日韩精品久久久久久密| 内地一区二区视频在线| 国产男靠女视频免费网站| 中亚洲国语对白在线视频| 真人做人爱边吃奶动态| 男女做爰动态图高潮gif福利片| 人妻少妇偷人精品九色| 欧美区成人在线视频| 听说在线观看完整版免费高清| 天堂影院成人在线观看| eeuss影院久久| 别揉我奶头~嗯~啊~动态视频| 欧美日韩瑟瑟在线播放| 黄色视频,在线免费观看| 高清在线国产一区| 亚洲国产精品sss在线观看| 毛片女人毛片| 亚洲专区国产一区二区| 亚洲七黄色美女视频| 大又大粗又爽又黄少妇毛片口| 一级a爱片免费观看的视频| 久久午夜亚洲精品久久| 小说图片视频综合网站| 女人十人毛片免费观看3o分钟| 亚洲三级黄色毛片| 久久精品国产清高在天天线| 人妻久久中文字幕网| 人妻夜夜爽99麻豆av| 99国产精品一区二区蜜桃av| 18+在线观看网站| 我的女老师完整版在线观看| 啦啦啦啦在线视频资源| 亚洲 国产 在线| 一边摸一边抽搐一进一小说| 国产精品亚洲美女久久久| 国产高清视频在线观看网站| 国内久久婷婷六月综合欲色啪| 精品人妻1区二区| 性色avwww在线观看| 欧美+日韩+精品| 99在线人妻在线中文字幕| 中文字幕免费在线视频6| 91午夜精品亚洲一区二区三区 | 亚洲va日本ⅴa欧美va伊人久久| 亚洲av美国av| 精品日产1卡2卡| 99久久久亚洲精品蜜臀av| 亚洲av熟女| 中文字幕av在线有码专区| 久久久久久伊人网av| 99热这里只有是精品在线观看| 亚洲中文字幕一区二区三区有码在线看| 91在线精品国自产拍蜜月| 日韩一区二区视频免费看| 亚洲美女视频黄频| 国产v大片淫在线免费观看| 久久热精品热| 99在线视频只有这里精品首页| 国产黄片美女视频| 国产高清有码在线观看视频| 亚洲无线观看免费| 亚洲av日韩精品久久久久久密| 老师上课跳d突然被开到最大视频| 亚洲精品日韩av片在线观看| x7x7x7水蜜桃| 成人毛片a级毛片在线播放| 99国产精品一区二区蜜桃av| 老司机深夜福利视频在线观看| 亚洲av成人av| 高清在线国产一区| 国产高清激情床上av| 在线看三级毛片| 精品一区二区三区人妻视频| 免费观看的影片在线观看| 天天一区二区日本电影三级| 男女边吃奶边做爰视频| 欧美3d第一页| 观看美女的网站| 国产人妻一区二区三区在| 看黄色毛片网站| eeuss影院久久| 男插女下体视频免费在线播放| 成人二区视频| 国产白丝娇喘喷水9色精品| 久久精品国产清高在天天线| 久久久久国内视频| 欧美成人免费av一区二区三区| 一个人看的www免费观看视频| 国产精品精品国产色婷婷| 久久久精品大字幕| 国产精品人妻久久久影院| 亚洲国产精品久久男人天堂| 欧美bdsm另类| 91麻豆av在线| 日韩欧美精品v在线| 国产精品一区www在线观看 | 久久久久九九精品影院| 99视频精品全部免费 在线| 女人十人毛片免费观看3o分钟| 久久久久免费精品人妻一区二区| 噜噜噜噜噜久久久久久91| 国产精品电影一区二区三区| 变态另类丝袜制服| 淫妇啪啪啪对白视频| 国产免费av片在线观看野外av| 成人二区视频| 婷婷精品国产亚洲av在线| 午夜爱爱视频在线播放| 天堂√8在线中文| 蜜桃亚洲精品一区二区三区| 国产精品福利在线免费观看| 美女高潮的动态| 欧美+亚洲+日韩+国产| 国内精品宾馆在线| 日本黄大片高清| 嫁个100分男人电影在线观看| 国产精品98久久久久久宅男小说| 亚洲黑人精品在线| 日韩精品中文字幕看吧| av视频在线观看入口| 亚洲,欧美,日韩| 丰满乱子伦码专区| 色综合婷婷激情| АⅤ资源中文在线天堂| 人妻丰满熟妇av一区二区三区| 高清日韩中文字幕在线| 真实男女啪啪啪动态图| 在线免费观看的www视频| 在线播放国产精品三级| 精品一区二区免费观看| 级片在线观看| 国产91精品成人一区二区三区| 欧美高清性xxxxhd video| 一级a爱片免费观看的视频| 美女高潮的动态| 最好的美女福利视频网| 99热这里只有是精品在线观看| 国产三级中文精品| 免费在线观看成人毛片| 午夜影院日韩av| 久久99热6这里只有精品| 此物有八面人人有两片| 久久久久性生活片| 亚洲成人精品中文字幕电影| 91麻豆av在线| 搡女人真爽免费视频火全软件 | 精品久久久久久久久久免费视频| 日韩欧美一区二区三区在线观看| 久久精品国产亚洲av天美| 欧美bdsm另类| 久久久久久久久久久丰满 | 中出人妻视频一区二区| 五月伊人婷婷丁香| 精品99又大又爽又粗少妇毛片 | 精品久久久久久久久久免费视频| 高清日韩中文字幕在线| 淫妇啪啪啪对白视频| 88av欧美| 亚洲三级黄色毛片| 亚洲av日韩精品久久久久久密| 淫妇啪啪啪对白视频| 欧美日韩精品成人综合77777| 亚洲黑人精品在线| 综合色av麻豆| 中文在线观看免费www的网站| 午夜福利18| 免费观看的影片在线观看| 亚洲经典国产精华液单| 99九九线精品视频在线观看视频| 女的被弄到高潮叫床怎么办 | 亚州av有码| 麻豆av噜噜一区二区三区| av福利片在线观看| 波多野结衣高清作品| 麻豆av噜噜一区二区三区| 亚洲精品一卡2卡三卡4卡5卡| 精品久久久久久久久av| xxxwww97欧美| 国产视频一区二区在线看| av在线天堂中文字幕| 精品久久久久久久久久免费视频| 久久国产乱子免费精品| 麻豆国产av国片精品| 婷婷亚洲欧美| 国产一区二区亚洲精品在线观看| 中文字幕免费在线视频6| 日韩一本色道免费dvd| 欧美xxxx黑人xx丫x性爽| 91狼人影院| 九色成人免费人妻av| 九九爱精品视频在线观看| 男人舔奶头视频| 99久久精品国产国产毛片| 日韩精品中文字幕看吧| 嫩草影院精品99| 欧美一级a爱片免费观看看| 婷婷色综合大香蕉| 色av中文字幕| 国产蜜桃级精品一区二区三区| 精品人妻一区二区三区麻豆 | 欧美一区二区国产精品久久精品| 久久久久免费精品人妻一区二区| 女人十人毛片免费观看3o分钟| 国产白丝娇喘喷水9色精品| av视频在线观看入口| 久久精品国产亚洲av涩爱 | 人人妻人人澡欧美一区二区| 国产精品久久久久久av不卡| 国产欧美日韩精品一区二区| 欧美bdsm另类| 女人被狂操c到高潮| 国产乱人视频| 99久久精品国产国产毛片| 一个人看的www免费观看视频| 国产伦人伦偷精品视频| 国产午夜精品论理片| 综合色av麻豆| 亚洲av中文字字幕乱码综合| 天堂√8在线中文| 久久精品国产亚洲av天美| 国产精品福利在线免费观看| 欧美xxxx性猛交bbbb| ponron亚洲| 男人和女人高潮做爰伦理| 亚洲国产精品合色在线| 在线免费观看的www视频| 中文在线观看免费www的网站| 国内揄拍国产精品人妻在线| 精品午夜福利视频在线观看一区| 最好的美女福利视频网| 成年免费大片在线观看| 午夜福利成人在线免费观看| 亚洲成人久久爱视频| 婷婷精品国产亚洲av| 免费看美女性在线毛片视频| 大又大粗又爽又黄少妇毛片口| 国产一区二区三区av在线 | 麻豆久久精品国产亚洲av| 嫩草影院入口| 伦精品一区二区三区| 丰满乱子伦码专区| 亚洲欧美日韩高清在线视频| 久久久久久久久大av| 欧美另类亚洲清纯唯美| av中文乱码字幕在线| 成人av在线播放网站| 全区人妻精品视频| 丝袜美腿在线中文| 精品久久久久久久久久免费视频| 小说图片视频综合网站| 亚洲av熟女| 国产高清不卡午夜福利| 乱系列少妇在线播放| 免费看a级黄色片| 国产精品久久电影中文字幕| 五月玫瑰六月丁香| 亚洲午夜理论影院| 一个人观看的视频www高清免费观看| 亚洲中文字幕日韩| 精品无人区乱码1区二区| 精品久久久久久久久亚洲 | 成人特级av手机在线观看| 欧美+日韩+精品| 亚洲精华国产精华精| 我的老师免费观看完整版| 99热6这里只有精品| 精品人妻视频免费看| 久久久久精品国产欧美久久久| 亚洲18禁久久av| 免费一级毛片在线播放高清视频| 国内毛片毛片毛片毛片毛片| 日本成人三级电影网站| 99精品久久久久人妻精品| 无遮挡黄片免费观看| 尤物成人国产欧美一区二区三区| 久久久久久久午夜电影| 亚洲欧美日韩东京热| 嫁个100分男人电影在线观看| 日韩一本色道免费dvd| 欧美成人免费av一区二区三区| 国产精品免费一区二区三区在线| 国产精品人妻久久久久久| 丝袜美腿在线中文| 我的老师免费观看完整版| 啦啦啦啦在线视频资源| 久久中文看片网| 18+在线观看网站| 99久久无色码亚洲精品果冻| av在线蜜桃| 久久99热6这里只有精品| 国产免费av片在线观看野外av| 欧美潮喷喷水| 久久久久久久精品吃奶| 动漫黄色视频在线观看| 国产日本99.免费观看| 热99在线观看视频| 日日夜夜操网爽| 制服丝袜大香蕉在线| 三级男女做爰猛烈吃奶摸视频| 乱系列少妇在线播放| 日韩中字成人| 免费人成在线观看视频色| 嫩草影院新地址| 在线观看66精品国产| 日韩 亚洲 欧美在线| 在线观看一区二区三区| 大又大粗又爽又黄少妇毛片口| 制服丝袜大香蕉在线| 一区二区三区激情视频| 亚洲成a人片在线一区二区| 亚洲综合色惰| 看黄色毛片网站| 欧美成人a在线观看| 日本免费一区二区三区高清不卡| 中文字幕免费在线视频6| 天堂√8在线中文| 亚洲欧美日韩东京热| 久久久色成人| 日韩中字成人| 国产又黄又爽又无遮挡在线| 五月玫瑰六月丁香| 精品99又大又爽又粗少妇毛片 | 看十八女毛片水多多多| 天美传媒精品一区二区| 伦精品一区二区三区| 一夜夜www| 黄色配什么色好看| 校园人妻丝袜中文字幕| 99久久九九国产精品国产免费| 一个人观看的视频www高清免费观看| 国产精品亚洲美女久久久| 精品国产三级普通话版| 亚洲avbb在线观看| 99久久精品热视频| 99热6这里只有精品| 美女免费视频网站| 中文在线观看免费www的网站| 久久精品国产清高在天天线| 久9热在线精品视频| 国内少妇人妻偷人精品xxx网站| x7x7x7水蜜桃| 国产精品综合久久久久久久免费| 日韩欧美精品免费久久| 日韩中文字幕欧美一区二区| 禁无遮挡网站| 日韩欧美 国产精品| 97超级碰碰碰精品色视频在线观看| 日韩欧美精品v在线| 亚洲avbb在线观看| 国产真实伦视频高清在线观看 | 国产一区二区三区av在线 | 久久国产乱子免费精品| 亚洲人成网站在线播| 国产主播在线观看一区二区| 亚洲精品乱码久久久v下载方式| 国产精品久久久久久精品电影| 观看免费一级毛片| 变态另类丝袜制服| 欧美日韩国产亚洲二区| 色播亚洲综合网| 亚洲av五月六月丁香网| a在线观看视频网站| 12—13女人毛片做爰片一| 在线看三级毛片| 精品久久久久久久久亚洲 | 久久精品国产清高在天天线| 淫妇啪啪啪对白视频| 成年女人永久免费观看视频| 日韩在线高清观看一区二区三区 | 成熟少妇高潮喷水视频| www.色视频.com| 精品久久久噜噜| 淫秽高清视频在线观看| 欧美绝顶高潮抽搐喷水| 99国产精品一区二区蜜桃av| 在线a可以看的网站| 男女视频在线观看网站免费| 日韩欧美国产一区二区入口| 国产色爽女视频免费观看| 狠狠狠狠99中文字幕| 国产精品综合久久久久久久免费| 久久久久久久久久久丰满 | a级一级毛片免费在线观看| 午夜福利在线观看吧| 亚洲美女黄片视频| 99精品在免费线老司机午夜| 男女视频在线观看网站免费| 最好的美女福利视频网| 免费观看精品视频网站| 性插视频无遮挡在线免费观看| 中文字幕av成人在线电影| 一本精品99久久精品77| 一个人免费在线观看电影| ponron亚洲| 国产aⅴ精品一区二区三区波| 日本爱情动作片www.在线观看 | 听说在线观看完整版免费高清| 熟女电影av网| 成熟少妇高潮喷水视频| 国产一区二区三区在线臀色熟女| 日本-黄色视频高清免费观看| 久久精品国产99精品国产亚洲性色| 长腿黑丝高跟| 一级黄色大片毛片| 欧美丝袜亚洲另类 | 尾随美女入室| 国产淫片久久久久久久久| 黄色女人牲交| 两性午夜刺激爽爽歪歪视频在线观看| 一进一出抽搐动态| 我的女老师完整版在线观看| 国产毛片a区久久久久| 国内精品久久久久久久电影| 国产亚洲91精品色在线|