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

    總強度磁異常各階垂向?qū)?shù)換算新方法

    2011-01-04 07:58:08翟國君卞光浪黃謨濤
    測繪學(xué)報 2011年6期
    關(guān)鍵詞:樣條二階導(dǎo)數(shù)

    翟國君,卞光浪,2,黃謨濤

    1.海軍海洋測繪研究所,天津300061;2.大連艦艇學(xué)院海測工程系,遼寧大連116018

    總強度磁異常各階垂向?qū)?shù)換算新方法

    翟國君1,卞光浪1,2,黃謨濤1

    1.海軍海洋測繪研究所,天津300061;2.大連艦艇學(xué)院海測工程系,遼寧大連116018

    總強度磁異常垂向?qū)?shù)換算在磁性目標(biāo)解釋推斷過程中具有重要意義。分析頻率域總強度磁異常各階垂向?qū)?shù)轉(zhuǎn)換因子的濾波特性,指出在求解總強度磁異常垂向?qū)?shù)時,常規(guī)傅里葉變換法會明顯放大觀測數(shù)據(jù)中高頻噪聲的影響,甚至可能淹沒掉真實信息。給出一個重要定理:總強度磁異常沿垂直方向的積分及各階偏導(dǎo)數(shù)均為準(zhǔn)調(diào)和函數(shù),在此基礎(chǔ)上,提出基于拉普拉斯方程的逐步求解法,聯(lián)合空間域和頻率域換算總強度磁異常沿垂直方向的各階導(dǎo)數(shù)。同時,采用球體磁場模型驗證所提方法的有效性,并推導(dǎo)地磁場方向和磁化強度方向不一致時,球體總強度磁異常沿垂直方向的一階和二階導(dǎo)數(shù)表達式。結(jié)論表明:利用樣條法計算的垂向?qū)?shù)結(jié)果精度明顯優(yōu)于常規(guī)傅里葉變換法計算結(jié)果,且具有較強的抗噪能力。

    總強度磁異常;垂向?qū)?shù);空間域和頻率域;拉普拉斯方程;球體磁場模型

    1 引 言

    磁測方法早期主要應(yīng)用于地球物理和磁力勘探領(lǐng)域。近年來,隨著人類開發(fā)與利用海洋資源進程的加快,該方法以成本低、快速高效等突出優(yōu)點,在水下小尺度磁性目標(biāo)精密探測方面?zhèn)涫荜P(guān)注,特別是沉船、鐵錨、海底電纜、未爆炸物以及水下潛航器的定位與識別方面更具明顯優(yōu)勢[1-2]。

    磁異常導(dǎo)數(shù)被廣泛應(yīng)用于磁性目標(biāo)磁異常解釋,是壓制區(qū)域場、圈定局部場、分離疊加異常以及異常反演中常用的方法[3]。例如,利用解析信號法和歐拉反褶積法確定磁性目標(biāo)平面位置,或者采用泰勒級數(shù)法實現(xiàn)磁場空間穩(wěn)健延拓,都需要以高精度磁異常水平和垂直方向的一階、二階甚至高階導(dǎo)數(shù)作為基礎(chǔ)數(shù)據(jù)。現(xiàn)階段,獲取磁場垂向?qū)?shù)主要有兩種方式,其一是將磁力傳感器設(shè)計為磁陣列模式[4],利用微分原理直接計算磁場的垂向?qū)?shù),該模式對載體運行姿態(tài)要求較高,海洋上受風(fēng)、流、浪等因素影響,實現(xiàn)起來較為困難;其二是將采集的剖面或平面磁場信息進行數(shù)據(jù)轉(zhuǎn)換處理,這種轉(zhuǎn)換既可在空間域內(nèi)進行,也可在頻率域內(nèi)進行。在空間域換算方面,文獻[5]嘗試將樣條函數(shù)插值法引入位場垂向?qū)?shù)解算中,但計算精度有限。目前普遍采用的磁異常垂向?qū)?shù)換算方法是傅里葉變換法[6](簡記為“常規(guī)FFT法”),該方法理論較為完善,但其頻率響應(yīng)函數(shù)相當(dāng)于高頻放大器,對觀測數(shù)據(jù)中的噪聲具有顯著的放大作用,當(dāng)觀測數(shù)據(jù)含有一定誤差時,計算精度較低,且存在較強的邊界效應(yīng)。為避免復(fù)雜的復(fù)數(shù)運算,文獻[7]提出采用離散余弦變換(DCT)計算磁異常導(dǎo)數(shù),該方法去除了原信號的相關(guān)性,在二度磁性目標(biāo)磁異常垂向?qū)?shù)換算中應(yīng)用效果較好,關(guān)于該方法能否適用于三度磁性目標(biāo)磁異常垂向?qū)?shù)換算還有待于進一步研究。

    針對總強度磁異常沿垂直方向?qū)?shù)的精確求解問題,筆者給出了關(guān)于總強度磁異常的一個重要定理,并綜合利用空間域和頻率域運算,提出一種求解總強度磁異常各階垂向?qū)?shù)的新方法。文中詳細研究了該方法的基本原理和技術(shù)措施,通過仿真球體磁場模型檢驗了方法的有效性,并將其計算結(jié)果與常規(guī)FFT法進行了分析比較。

    2 常規(guī)FFT法

    由傅里葉變換,三度磁性目標(biāo)總強度磁異常Bm的一階垂向?qū)?shù)Bmz與其頻譜關(guān)系[8]為

    根據(jù)位場理論,總強度磁異常近似滿足拉普拉斯方程,利用外部Dirichlet問題條件,其空間延拓公式可表示為

    式中,F(xiàn)Bm(u,v,0)為z=0平面上總強度磁異常Bm(x,y,0)的頻譜。

    將式(2)對z方向求偏導(dǎo)數(shù),得

    由于同一函數(shù)有相同的傅里葉表達式,對比式(1)和式(3)兩端,得

    當(dāng)z=0時,式(4)即為平面上總強度磁異常頻譜向其一階垂向?qū)?shù)頻譜轉(zhuǎn)換的關(guān)系式,同理,還可推導(dǎo)出總強度磁異常沿垂向的第k階導(dǎo)數(shù)頻譜為這里稱為總強度磁異常與其第k階垂向?qū)?shù)間的頻率域轉(zhuǎn)換因子(也稱頻率響應(yīng)函數(shù))。計算時,首先將平面上總強度磁異常格網(wǎng)數(shù)據(jù)利用快速傅里葉變換轉(zhuǎn)換為頻譜FBm(u,v,0),再乘以垂向?qū)?shù)頻率域轉(zhuǎn)換因子獲取第k階垂向?qū)?shù)的頻譜,最后通過反傅里葉變換即可得到總強度磁異常各階垂向?qū)?shù)。

    3 基于拉普拉斯方程的逐步求解法

    為減弱高頻噪聲對總強度磁異常垂向?qū)?shù)計算的影響,這里提出基于拉普拉斯方程的逐步求解法。在求解過程中,將運用總強度磁異常的某些特殊性質(zhì),限于篇幅,略去推導(dǎo)過程,給出如下定理。

    定理:在無源空間內(nèi),總強度磁異常沿垂直方向的積分及任意階偏導(dǎo)數(shù)均為準(zhǔn)調(diào)和函數(shù)。

    結(jié)合上述定理,總強度磁異常各階垂向?qū)?shù)的求解步驟如下:

    (1)對z=0平面上總強度磁異常數(shù)據(jù)進行傅里葉變換,設(shè)其在頻率域頻譜為FBm(u,v,0),由式(4),F(xiàn)Bm(u,v,0)乘以頻率域轉(zhuǎn)換因子即為總強度磁異常沿垂直方向積分的頻譜,再利用反傅里葉變換得到積分值J(x,y,0);

    (2)計算J(x,y,0)沿x方向和y方向的二階偏導(dǎo)數(shù)Jxx(x,y,0)與Jyy(x,y,0),由定理知,總強度磁異常沿垂直方向的積分是準(zhǔn)調(diào)和函數(shù),滿足拉普拉斯方程,則總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)為

    (3)計算Bm在水平方向的兩個二階偏導(dǎo)數(shù)與而總強度磁異常自身也為準(zhǔn)調(diào)和函數(shù),同樣利用拉普拉斯方程,得總強度磁異常沿垂直方向的二階偏導(dǎo)數(shù)為

    (4)由于總強度磁異常沿垂直方向的任意階偏導(dǎo)數(shù)都是準(zhǔn)調(diào)和函數(shù),得第k(k=3,4,…,n)階垂直導(dǎo)數(shù)為

    式中,Bmkz為總強度磁異常沿垂直方向的第k階偏導(dǎo)數(shù)和分別是Bm(k-2)z沿x方向和y方向的二階偏導(dǎo)數(shù)。

    4 調(diào)和函數(shù)二階水平方向?qū)?shù)計算方法

    目前,計算調(diào)和函數(shù)沿水平方向的二階偏導(dǎo)數(shù)通常在頻率域內(nèi)進行轉(zhuǎn)換(計算步驟與常規(guī)FFT法計算垂向?qū)?shù)類似),但該方法存在的不足是會顯著放大觀測數(shù)據(jù)中噪聲的影響,且具有較強的邊界效應(yīng)。在數(shù)值計算中,三點二階中心差分法和三次樣條曲線函數(shù)法經(jīng)常被用于計算函數(shù)在一維方向的二階導(dǎo)數(shù),本文根據(jù)平面磁場數(shù)據(jù)的二維空間分布特點,將兩者引入并應(yīng)用到調(diào)和函數(shù)二階水平方向偏導(dǎo)數(shù)計算中。下面以總強度磁異常為例,給出三點二階中心差分法和雙三次樣條曲線函數(shù)法計算其沿水平方向二階偏導(dǎo)數(shù)的具體方法原理。

    4.1 三點二階中心差分法

    假設(shè)x1,x2,…,xM為格網(wǎng)化數(shù)據(jù)沿x方向的等間距節(jié)點,且x1<x2<…<xM,xi-xi-1=Δx。則在(xi,y)節(jié)點處總強度磁異常沿x方向的二階偏導(dǎo)數(shù)可表示為

    類似地,設(shè)y1,y2,…,yN為格網(wǎng)化數(shù)據(jù)沿y方向的等間距節(jié)點,且y1<y2<…<yN,yiyi-1=Δy。則在(x,yi)節(jié)點處總強度磁異常沿y方向的二階偏導(dǎo)數(shù)可表示為

    值得注意的是,利用式(8)和式(9)求解總強度磁異常沿水平方向二階偏導(dǎo)數(shù)時,邊界處的二階導(dǎo)數(shù)需要附加條件,本文仿真計算時,考慮到僅對比梯度變化較大區(qū)域的計算精度,所以未計算邊界部分的二階偏導(dǎo)數(shù)。

    4.2 雙三次樣條曲線函數(shù)法

    雙三次樣條曲線函數(shù)是對x方向和y方向分別采用三次樣條曲線函數(shù)進行擬合[9]。以x方向測線為例,設(shè)有觀測數(shù)據(jù) {x1,Bm(x1)},{x2,Bm(x2)},…,{xM,Bm(xM)},且x1<x2<…<xM,整個區(qū)間[x1,xM]被節(jié)點分成M-1個小區(qū)間[x1,x2]、[x2,x3]、…、[xM-1,xM],現(xiàn)要求構(gòu)造函數(shù)S(x)滿足以下三個條件:

    (1)插值條件

    (2)在每個小區(qū)間[xi,xi+1](i=1,2,…,M-1)上,S(x)是x的三次多項式,即

    (3)在整個區(qū)間[x1,xM]上,S(x)具有連續(xù)的二階導(dǎo)數(shù),在區(qū)間的每個節(jié)點xi(i=1,2,…,M)上,函數(shù)值、一階導(dǎo)數(shù)和二階導(dǎo)數(shù)連續(xù),即

    式中,i=2,3,…,M-1。

    利用雙三次樣條曲線函數(shù)計算水平方向二階偏導(dǎo)數(shù)時,需要給定兩個邊界約束條件,本文選擇自然三次樣條插值函數(shù),即Sxx(x1)=Sxx(xM)=0,結(jié)合方程組(12)得到的Sxx(xi)即為xi(i=2,3,…,M-1)節(jié)點處總強度磁異常沿x方向的二階偏導(dǎo)數(shù)。按照相同的方法,可進一步得到各測點處總強度磁異常沿y方向的二階偏導(dǎo)數(shù)。

    5 仿真試驗分析

    為檢驗提出的總強度磁異常各階垂向?qū)?shù)換算方法的有效性,仿真三個磁化球體磁異常數(shù)據(jù)進行分析。

    仿真測區(qū)面積為295m×295m,測線和測點間距均為5m。地磁場傾角和偏角分別為45°與5°,為檢驗方法的普遍適用性,這里設(shè)球體的磁化強度方向和地磁場方向不一致,球體空間位置及磁性參數(shù)分別如表1所示。

    表1 仿真球體模型相關(guān)參數(shù)Tab.1 Parameters of synthetic sphere model

    為定量分析總強度磁異常的垂向一階導(dǎo)數(shù)、二階導(dǎo)數(shù)計算結(jié)果精度,在文獻[8]給出的球體總強度磁異常表達式基礎(chǔ)上,筆者進一步推導(dǎo)磁化強度方向和地磁場方向不一致時,磁化球體總強度磁異常垂向一階導(dǎo)數(shù)Bmz和二階導(dǎo)數(shù)Bmzz解析表達式,分別如式(13)與式(14)所示

    式中,I0和A0分別為地磁場傾角和偏角;Bxz、Byz、Bzz和Bxzz、Byzz、Bzzz分別是Bv(v=x,y,z)沿垂直方向的一階偏導(dǎo)數(shù)與二階偏導(dǎo)數(shù),其中,

    根據(jù)球體總強度磁異常以及式(13)和式(14),仿真得到z=0m平面上總強度磁異常及其垂向一階和二階偏導(dǎo)數(shù)等值線如圖1所示。

    圖1 總強度磁異常及其垂向一階和二階偏導(dǎo)數(shù)理論等值線圖Fig.1 Contour map of theoretical TMA as well as its 1st and 2ed vertical derivatives

    圖1(a)中總強度磁異常最大值為48.96nT,最小值為-20.69nT,平均值為1.05nT;圖1(b)中總強度磁異常垂向一階導(dǎo)數(shù)最大值為5.30nT/m,最小值為-2.18nT/m,平均值為0.001 8nT/m;圖1(c)中總強度磁異常垂向二階導(dǎo)數(shù)最大值為0.70nT/m2,最小值為-0.30nT/m2,平均值為-0.000 1nT/m2。

    本文在前面給出了兩種垂向?qū)?shù)計算方法,分別是常規(guī)FFT法和基于拉普拉斯方程的逐步求解法。在利用后者求解垂向?qū)?shù)時,提供了兩種求解水平方向二階偏導(dǎo)數(shù)的算法,簡稱為“差分法”和“樣條法”。為充分驗證提出方法的適用性,這里分無噪聲和有噪聲兩種情況分別進行分析。

    5.1 無噪聲情況

    圖2和圖3是根據(jù)圖1(a)中的總強度磁異常數(shù)據(jù),利用常規(guī)FFT法、差分法和樣條法換算得到總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)和二階偏導(dǎo)數(shù)等值線圖。

    圖2 三種方法計算得到的總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)等值線圖(等值線間隔0.2nT/m)Fig.2 Contour map of calculated the 1st vertical derivative of TMA with three methods(contour interval is 0.2nT/m)

    圖3 三種方法計算得到的總強度磁異常沿垂直方向的二階偏導(dǎo)數(shù)等值線圖(等值線間隔0.05nT/m2)Fig.3 Contour map of calculated the 2ed vertical derivative of TMA with three methods(contour interval is 0.05nT/m2)

    從圖2和圖3可以看出,無噪聲情況下,三種換算方法得到的垂向一階偏導(dǎo)數(shù)和二階偏導(dǎo)數(shù)等值線圖與相應(yīng)的理論等值線圖實現(xiàn)了很好的吻合,但常規(guī)FFT法隨著導(dǎo)數(shù)階數(shù)的增加,邊界效應(yīng)越來越明顯。為定量計算總強度磁異常各階垂向?qū)?shù)換算精度,并排除邊界效應(yīng)對統(tǒng)計結(jié)果的影響,選取上述各圖中虛線框D域內(nèi)梯度變化較大的格網(wǎng)點數(shù)據(jù)進行分析,即D={(x,y)|25m≤x≤270m,25m≤y≤270m},采用的精度評價指標(biāo)公式為式(15)。

    式中,^Bmkz(i,j)和Bmkz(i,j)(k為正整數(shù))分別表示格網(wǎng)點上總強度磁異常第k階垂向?qū)?shù)計算值與理論值;Nx和Ny分別為D域內(nèi)x方向和y方向采樣點數(shù);ρk為第k階垂向?qū)?shù)計算結(jié)果中誤差。

    三種方法換算結(jié)果統(tǒng)計情況如表2所示。

    表2 三種方法換算的垂向一階和二階導(dǎo)數(shù)結(jié)果統(tǒng)計Tab.2 Statistics of the 1st and 2ed vertical derivative results with three methods

    從表2可以看到,三種垂向?qū)?shù)換算方法計算結(jié)果均具有較高的精度,其中,樣條法和常規(guī)FFT法計算結(jié)果更加接近于理論值。由于式(15)給出的精度評價指標(biāo)屬于絕對誤差,其大小與各階垂向?qū)?shù)的振幅緊密相關(guān),垂向?qū)?shù)的振幅越大,絕對誤差越大;反之,振幅越小,絕對誤差則越小。因此,該指標(biāo)在一定程度上無法全面反映換算結(jié)果的優(yōu)劣,故進一步采用中誤差與理論值振幅之比作為評價指標(biāo)

    以樣條法為例,總強度磁異常垂向一階導(dǎo)數(shù)和二階導(dǎo)數(shù)理論值振幅分別為7.48nT/m和1nT/m2,計算得到的垂向一階、二階導(dǎo)數(shù)中誤差和振幅之比分別為0.11%與0.17%,驗證了該方法具有很高的垂向?qū)?shù)換算精度,且二階垂向?qū)?shù)換算結(jié)果的相對精度要稍低于一階垂向?qū)?shù)換算結(jié)果。因此,當(dāng)實測資料中無噪聲或噪聲微小時,既可采用常規(guī)FFT法,也可采用樣條法計算總強度磁異常的垂向一階導(dǎo)數(shù)和二階導(dǎo)數(shù)。

    考慮到利用泰勒級數(shù)法進行總強度磁異常空間延拓時,往往還需要獲取總強度磁異常的高階垂向?qū)?shù),本文進一步對比分析了常規(guī)FFT法和樣條法的垂向三階、四階和五階導(dǎo)數(shù)換算結(jié)果。結(jié)果表明:常規(guī)FFT法計算結(jié)果存在較強的邊界效應(yīng),且隨著垂向?qū)?shù)階數(shù)的增加,其邊界效應(yīng)更加明顯,到第五階垂向?qū)?shù)時,誤差已非常大;而樣條法雖然也存在一定的邊界效應(yīng),但其換算效果要明顯優(yōu)于常規(guī)FFT法結(jié)果。因此,實踐中建議采用樣條法計算總強度磁異常高階垂向?qū)?shù)。

    5.2 含有噪聲情況

    為更接近實際情況,進一步檢驗提出方法的穩(wěn)健性,在z=0m平面上總強度磁異常理論值中加入2%高斯白噪聲(振幅為1.39nT),仿真得到含噪聲的總強度磁異常等值線如圖4所示。

    圖4 含有高斯噪聲的總強度磁異常等值線間隔(等值線間隔2nT)Fig.4 Contour map of TMA corrupted with Gauss noise(contour interval is 2nT)

    圖5和圖6是用含有高斯噪聲的總強度磁異常數(shù)據(jù)換算得到其沿垂直方向的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)等值線圖。

    如圖5和圖6所示,當(dāng)總強度磁異常數(shù)據(jù)中加入2%高斯白噪聲后,垂向一階導(dǎo)數(shù)和二階導(dǎo)數(shù)計算結(jié)果精度要明顯低于無噪聲情況下的結(jié)果,表明噪聲對垂向?qū)?shù)換算結(jié)果存在較大影響。在三種方法中,樣條法計算結(jié)果精度最高,根據(jù)D域內(nèi)格網(wǎng)點數(shù)據(jù),利用式(15)和式(16)計算得到的ρ1=0.07nT/m,η1=0.93%,ρ2=0.05nT/m2,η2=5%;差分法計算結(jié)果精度次之,對應(yīng)ρ1=0.13nT/m,η1=1.7%,ρ2=0.11nT/m2,η2=11%;常規(guī)FFT法計算結(jié)果精度最低,其中ρ1=0.21nT/m,η1=2.8%,ρ2=0.13nT/m2,η2=13%。從相對精度指標(biāo)變化情況來看,垂向二階導(dǎo)數(shù)較一階導(dǎo)數(shù)受噪聲的影響更加嚴(yán)重。

    圖5 三種方法計算得到的總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)等值線圖(等值線間隔0.2nT/m)Fig.5 Contour map of calculated the 1st vertical derivative of TMA with three methods(contour interval is 0.2nT/m)

    這里列出零噪聲、1%、2%、3%、5%和10%等不同噪聲水平情況下,三種垂向?qū)?shù)方法計算得到的垂向一階、二階導(dǎo)數(shù)結(jié)果精度。

    圖6 三種方法計算得到的總強度磁異常沿垂直方向的二階偏導(dǎo)數(shù)等值線圖(等值線間隔0.05nT/m2)Fig.6 Contour map of calculated the 2ed vertical derivative of TMA with three methods(contour interval is 0.05nT/m2)

    從表3可以看出,當(dāng)總強度磁異常觀測數(shù)據(jù)中含有噪聲時,同一噪聲水平情況下,樣條法計算結(jié)果精度最高,差分法次之,常規(guī)FFT法精度最低;而對于同一種方法,垂向二階導(dǎo)數(shù)計算結(jié)果相對精度要低于一階導(dǎo)數(shù)計算結(jié)果相對精度。同時,隨著噪聲水平的逐漸增大,各階垂向?qū)?shù)換算結(jié)果精度也不斷降低。

    表3 不同噪聲水平下三種方法換算得到的垂向?qū)?shù)中誤差Tab.3 Mean square error of three methods with different noise level

    6 結(jié) 論

    (1)分析了常規(guī)FFT法在換算總強度磁異常各階垂向?qū)?shù)過程中存在的不足;給出了一個重要定理:在無源空間內(nèi),總強度磁異常沿垂直方向的積分和各階偏導(dǎo)數(shù)同樣為準(zhǔn)調(diào)和函數(shù),在此基礎(chǔ)上,提出了一種新的總強度磁異常各階垂向?qū)?shù)換算方法。

    (2)圍繞調(diào)和函數(shù)的二階水平方向偏導(dǎo)數(shù)求解問題,給出了三點二階中心差分法和雙三次樣條曲線函數(shù)法。通過仿真試驗對三點二階中心差分法和雙三次樣條曲線函數(shù)法應(yīng)用效果進行了比較分析,結(jié)果表明,利用樣條法換算得到的各階垂向?qū)?shù)精度較高,且具有較強的抗噪聲能力。

    (3)當(dāng)總強度磁異常觀測數(shù)據(jù)不含噪聲時,常規(guī)FFT法、差分法和樣條法換算的總強度磁異常垂向一階導(dǎo)數(shù)與二階導(dǎo)數(shù)結(jié)果精度大體相當(dāng),但隨著其垂向?qū)?shù)階數(shù)的增加,常規(guī)FFT法邊界效應(yīng)現(xiàn)象較其他兩種方法嚴(yán)重得多。當(dāng)總強度磁異常觀測數(shù)據(jù)含有噪聲時,隨著噪聲水平的不斷提高,垂向?qū)?shù)換算結(jié)果精度逐漸降低。但總體而言,樣條法垂向?qū)?shù)換算結(jié)果精度要明顯優(yōu)于常規(guī)FFT法換算結(jié)果。

    [1] GAO Ping,LESLIE M C.A Two-dimensional Generalized Likelihood Ratio Test for Land Mine and Small Unexploded Ordnance Detection[J].Signal Processing,2000,80(8):1669-1686.

    [2] YU Bo,ZHAI Guojun,LIU Yanchun,et al.The Downward Continuation Method of Aeromagnetic Data to the Sea Level[J].Acta Geodaetica et Cartographica Sinica,2009,38(3):202-209.(于波,翟國君,劉雁春,等.利用航磁數(shù)據(jù)向下延拓得到海面磁場的方法[J].測繪學(xué)報,2009,38(3):202-209.)

    [3] WEN B D,SHU K H,YI C Y.A Derivative-based Interpretation Approach to Estimating Source Parameters of Simple 2DMagnetic Sources from Euler Deconvolution,the Analytic-signal Method and Analytical Expressions of the Anomalies[J].Geophysical Prospecting,2007,55(2):255-264.

    [4] SALEM A,HAMADA T,ASAHINA J K.Detection of Unexploded Ordnance(UXO)Using Marine Magnetic Gradiometer Data[J].Exploration Geophysics,2005,58(1):97-103.

    [5] WANG Bingzhu.Computing the Vertical Second Derivative and Upward Continuation of Gravity Anomaly by Spline Function Method[J].Oil Geophysical Prospecting,1996,31(3):415-422.(汪炳柱.用樣條函數(shù)法求重力異常垂向二階導(dǎo)數(shù)和向上延拓計算[J].石油地球物理勘探,1996,31(3):415-422.)

    [6] BLAKELY R J.Potential Theory in Gravity and Magnetic Applications[M].Cambridge:Cambridge University Press,1995.

    [7] ZHANG Fengxu,ZHANG Fengqin,MENG Lingshun,et al.Magnetic Potential Spectrum Analysis and Calculating Method of Magnetic Anomaly Derivatives Based on Discrete Cosine Transform[J].Chinese Journal of Geophysics,2007,50(1):297-304.(張鳳旭,張鳳琴,孟令順,等.基于離散余弦變換的磁位譜分析及磁異常導(dǎo)數(shù)計算方法[J].地球物理學(xué)報,2007,50(1):297-304.)

    [8] GUAN Zhining.Magnetic Field and Magnetic Exploration[M].Beijing:Geological Publishing House,2005.(管志寧.地磁場與磁力勘探[M].北京:地質(zhì)出版社,2005.)

    [9] YAO Changli,HUANG Weining,GUAN Zhining.Fast Splines Conversion of Curved-surface Potential Field and Vertical Gradient Data into Horizontal-plane Data[J].Oil Geophysical Prospecting,1997,32(2):229-236.(姚長利,黃衛(wèi)寧,管志寧.綜合利用位場及其垂直梯度的快速樣條曲化平方法[J].石油地球物理勘探,1997,32(2):229-236.)

    A New Method to Calculate the Vertical Derivatives of Total Field Magnetic Anoma-

    lyZHAI Guojun1,BIAN Guanglang1,2,HUANG Motao1
    1.Naval Institute of Hydrographic Surveying and Charting,Tianjin300061,China;2.Department of Hydrography and Cartography,Dalian Naval Academy,Dalian 116018,China

    The vertical derivative of the total field magnetic anomaly(TMA)plays an important role in the interpretation process of magnetic objects.Owing to the filter characteristics of vertical derivative operator in frequency domain,the conversion of vertical derivative is inherently unstable and any high-frequency noise presented in the surveying data gets strongly magnified in the transformed map in such a way to mask any useful signal.Based on the harmonic properties of the vertical integral and derivatives of TMA,an algorithm is presented to perform calculation of vertical derivatives using both frequency and space domain transformations.The effectiveness of the suggested techniques has been illustrated by synthetic sphere magnetic model whose TMA’s first and second vertical derivative expressions are deduced when the geomagnetic direction is different from magnetization direction.The conclusion indicates that the presented vertical derivatives calculation method provides better results and allows a lower degrading of the signal-to-noise ratio than the standard Fourier method.

    total field magnetic anomaly;vertical derivative;space and frequency domain;Laplace equation;sphere magnetic model

    ZHAI Guojun(1961—),male,PhD,senior engineer,PhD supervisor,majors in technology and application of bathymetry and marine geodesy.

    1001-1595(2011)06-0671-08

    P229

    A

    國家自然科學(xué)基金(40671161)

    雷秀麗)

    2010-12-28

    2011-03-01

    翟國君(1961—),男,博士,高級工程師,博士生導(dǎo)師,主要從事海底地形測量和海洋大地測量技術(shù)與應(yīng)用研究。

    E-mail:zhaigj@163.com

    猜你喜歡
    樣條二階導(dǎo)數(shù)
    一元五次B樣條擬插值研究
    解導(dǎo)數(shù)題的幾種構(gòu)造妙招
    一類二階迭代泛函微分方程的周期解
    一類二階中立隨機偏微分方程的吸引集和擬不變集
    二階線性微分方程的解法
    一類二階中立隨機偏微分方程的吸引集和擬不變集
    三次參數(shù)樣條在機床高速高精加工中的應(yīng)用
    三次樣條和二次刪除相輔助的WASD神經(jīng)網(wǎng)絡(luò)與日本人口預(yù)測
    軟件(2017年6期)2017-09-23 20:56:27
    基于樣條函數(shù)的高精度電子秤設(shè)計
    關(guān)于導(dǎo)數(shù)解法
    男男h啪啪无遮挡| 看片在线看免费视频| 一本综合久久免费| 老熟妇仑乱视频hdxx| 欧美在线黄色| 丝袜美腿诱惑在线| 免费不卡黄色视频| av国产精品久久久久影院| av免费在线观看网站| av天堂在线播放| 99精品欧美一区二区三区四区| 久久人妻熟女aⅴ| 亚洲人成电影免费在线| 91老司机精品| 日韩中文字幕欧美一区二区| 久久天堂一区二区三区四区| 亚洲av熟女| 三级毛片av免费| 成人影院久久| 一级a爱视频在线免费观看| 成人国产一区最新在线观看| 亚洲熟妇熟女久久| 久久精品成人免费网站| 午夜成年电影在线免费观看| 国产精品乱码一区二三区的特点 | 国产91精品成人一区二区三区| 日本免费a在线| 国产欧美日韩综合在线一区二区| 99精品在免费线老司机午夜| 亚洲成人久久性| 国产av一区二区精品久久| 久久人妻福利社区极品人妻图片| 亚洲欧美激情综合另类| 精品国产乱子伦一区二区三区| 黄色a级毛片大全视频| 丝袜美足系列| 99久久久亚洲精品蜜臀av| 午夜福利一区二区在线看| 日韩视频一区二区在线观看| 日韩有码中文字幕| 十分钟在线观看高清视频www| 黑人猛操日本美女一级片| 9191精品国产免费久久| 最好的美女福利视频网| 久久中文字幕一级| 亚洲成国产人片在线观看| 99精品在免费线老司机午夜| 伊人久久大香线蕉亚洲五| 亚洲久久久国产精品| 99riav亚洲国产免费| 久久久国产一区二区| 美国免费a级毛片| 嫩草影视91久久| 日本黄色日本黄色录像| 精品熟女少妇八av免费久了| 精品高清国产在线一区| 午夜免费成人在线视频| videosex国产| 国产熟女午夜一区二区三区| 男男h啪啪无遮挡| 久久久久久人人人人人| 国产免费现黄频在线看| 99香蕉大伊视频| 黑人欧美特级aaaaaa片| 国产高清视频在线播放一区| 五月开心婷婷网| 国产av又大| 国产精品一区二区免费欧美| 黄网站色视频无遮挡免费观看| 在线观看免费视频网站a站| 午夜亚洲福利在线播放| 成熟少妇高潮喷水视频| 曰老女人黄片| 99re在线观看精品视频| 亚洲五月天丁香| 成熟少妇高潮喷水视频| 精品少妇一区二区三区视频日本电影| 18禁国产床啪视频网站| 香蕉国产在线看| 最新美女视频免费是黄的| 亚洲精品一卡2卡三卡4卡5卡| 久久国产精品男人的天堂亚洲| 国产精品av久久久久免费| 人人妻人人澡人人看| 亚洲色图综合在线观看| 午夜免费成人在线视频| 午夜福利,免费看| 少妇裸体淫交视频免费看高清 | 欧美性长视频在线观看| 精品一区二区三卡| 一级片免费观看大全| 精品日产1卡2卡| 高清欧美精品videossex| av中文乱码字幕在线| 国产精品自产拍在线观看55亚洲| 欧美丝袜亚洲另类 | 亚洲精品成人av观看孕妇| 精品人妻在线不人妻| 电影成人av| 欧美一级毛片孕妇| 制服人妻中文乱码| www日本在线高清视频| 在线视频色国产色| 久久99一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 欧美乱码精品一区二区三区| 精品一区二区三卡| 亚洲精华国产精华精| 一进一出抽搐gif免费好疼 | 91九色精品人成在线观看| 国产深夜福利视频在线观看| 免费观看人在逋| www国产在线视频色| 国产免费男女视频| 国产高清激情床上av| 夜夜爽天天搞| 亚洲情色 制服丝袜| 丰满饥渴人妻一区二区三| 男女床上黄色一级片免费看| 国产成人一区二区三区免费视频网站| 久久久久久久久中文| 99国产极品粉嫩在线观看| 90打野战视频偷拍视频| 桃色一区二区三区在线观看| 国产激情欧美一区二区| 久久久久久久久中文| 少妇裸体淫交视频免费看高清 | 日韩免费av在线播放| 国产精品国产av在线观看| 51午夜福利影视在线观看| 亚洲av第一区精品v没综合| 视频在线观看一区二区三区| 真人一进一出gif抽搐免费| www国产在线视频色| 我的亚洲天堂| 免费日韩欧美在线观看| 亚洲精品国产区一区二| 国产一区在线观看成人免费| 黄色视频不卡| 久久久久久久久久久久大奶| 精品久久久久久电影网| 搡老乐熟女国产| 高清欧美精品videossex| 欧美日韩亚洲高清精品| 亚洲精品久久成人aⅴ小说| 露出奶头的视频| 中国美女看黄片| 高清av免费在线| 精品久久久精品久久久| 国产不卡一卡二| 国产精品久久视频播放| 看片在线看免费视频| 老熟妇仑乱视频hdxx| 嫩草影院精品99| 久久精品国产99精品国产亚洲性色 | 男女之事视频高清在线观看| 久久人妻av系列| 免费在线观看影片大全网站| 亚洲欧美一区二区三区黑人| 久久国产精品人妻蜜桃| 国产高清视频在线播放一区| 男女下面插进去视频免费观看| 男女高潮啪啪啪动态图| 波多野结衣高清无吗| 国内毛片毛片毛片毛片毛片| 丰满迷人的少妇在线观看| 国产高清视频在线播放一区| 国产成人精品无人区| 伦理电影免费视频| 精品人妻1区二区| 欧美成人性av电影在线观看| 脱女人内裤的视频| 色在线成人网| 国产伦人伦偷精品视频| 精品久久久久久久久久免费视频 | 色精品久久人妻99蜜桃| 99国产精品99久久久久| 在线永久观看黄色视频| 日本一区二区免费在线视频| 国产精品1区2区在线观看.| 亚洲久久久国产精品| 精品一区二区三区四区五区乱码| 侵犯人妻中文字幕一二三四区| 欧美黑人欧美精品刺激| 久久精品国产99精品国产亚洲性色 | 黄色怎么调成土黄色| 人人妻,人人澡人人爽秒播| 国产有黄有色有爽视频| 最新美女视频免费是黄的| 久久人人97超碰香蕉20202| 国产精品av久久久久免费| 老司机福利观看| 天堂俺去俺来也www色官网| 视频区欧美日本亚洲| 国产单亲对白刺激| 黄色 视频免费看| av在线天堂中文字幕 | 久久欧美精品欧美久久欧美| av超薄肉色丝袜交足视频| 精品一区二区三区av网在线观看| 亚洲av成人不卡在线观看播放网| 亚洲一区二区三区欧美精品| 91在线观看av| 天堂动漫精品| 大码成人一级视频| ponron亚洲| 狠狠狠狠99中文字幕| 久久人妻av系列| 亚洲精品在线观看二区| 亚洲五月婷婷丁香| 又黄又爽又免费观看的视频| 视频区欧美日本亚洲| 神马国产精品三级电影在线观看 | 免费在线观看视频国产中文字幕亚洲| 国产精品二区激情视频| 国产亚洲欧美精品永久| 99久久人妻综合| 母亲3免费完整高清在线观看| 50天的宝宝边吃奶边哭怎么回事| 美女午夜性视频免费| 成年人黄色毛片网站| 精品高清国产在线一区| 亚洲精品一二三| 国产精品秋霞免费鲁丝片| 在线观看舔阴道视频| 国产主播在线观看一区二区| 桃色一区二区三区在线观看| 亚洲五月色婷婷综合| 国产一区二区激情短视频| 老司机福利观看| 麻豆久久精品国产亚洲av | 不卡av一区二区三区| 午夜免费观看网址| 亚洲欧美精品综合久久99| 亚洲精品成人av观看孕妇| 亚洲情色 制服丝袜| 日本五十路高清| 一级,二级,三级黄色视频| 91精品国产国语对白视频| 日韩 欧美 亚洲 中文字幕| 波多野结衣一区麻豆| 欧美国产精品va在线观看不卡| 精品乱码久久久久久99久播| 性少妇av在线| 精品少妇一区二区三区视频日本电影| 欧美日韩亚洲综合一区二区三区_| 国产高清videossex| 韩国精品一区二区三区| 1024视频免费在线观看| 欧美国产精品va在线观看不卡| 亚洲av成人av| 国产人伦9x9x在线观看| 久久国产精品男人的天堂亚洲| 波多野结衣一区麻豆| 老司机亚洲免费影院| 国产一区二区三区综合在线观看| 成人影院久久| 国产不卡一卡二| 免费在线观看黄色视频的| 97人妻天天添夜夜摸| 国产一区二区激情短视频| 日韩欧美在线二视频| 韩国精品一区二区三区| 国产欧美日韩一区二区精品| 在线观看66精品国产| 精品久久久久久久毛片微露脸| 久久久久精品国产欧美久久久| 国产av在哪里看| 日韩高清综合在线| 欧美日韩av久久| 久久精品人人爽人人爽视色| 久久午夜综合久久蜜桃| 在线观看一区二区三区| 成人国语在线视频| 999精品在线视频| 精品福利永久在线观看| 亚洲国产欧美网| 亚洲第一青青草原| 免费久久久久久久精品成人欧美视频| 无遮挡黄片免费观看| 69精品国产乱码久久久| 激情视频va一区二区三区| 久久精品国产亚洲av高清一级| 啦啦啦在线免费观看视频4| 国产精品久久久人人做人人爽| 久久这里只有精品19| 亚洲成人免费电影在线观看| 国产精品久久久久久人妻精品电影| 免费一级毛片在线播放高清视频 | av网站免费在线观看视频| 美女福利国产在线| 国产一区二区三区综合在线观看| 两个人看的免费小视频| 国产精品久久久人人做人人爽| 女人精品久久久久毛片| www.自偷自拍.com| 中文字幕人妻熟女乱码| 亚洲精品在线美女| 中国美女看黄片| 两性午夜刺激爽爽歪歪视频在线观看 | 久久香蕉激情| 在线观看66精品国产| 99国产综合亚洲精品| 精品国产乱子伦一区二区三区| 久久亚洲精品不卡| 欧美中文日本在线观看视频| 日本免费一区二区三区高清不卡 | 法律面前人人平等表现在哪些方面| 欧美+亚洲+日韩+国产| 亚洲五月婷婷丁香| 国产精品亚洲一级av第二区| 天堂中文最新版在线下载| 欧美中文日本在线观看视频| 亚洲成人久久性| 久久亚洲真实| 欧美性长视频在线观看| 中文字幕另类日韩欧美亚洲嫩草| 国内久久婷婷六月综合欲色啪| 一级黄色大片毛片| 亚洲aⅴ乱码一区二区在线播放 | 免费少妇av软件| 在线观看免费午夜福利视频| 丝袜在线中文字幕| 日韩人妻精品一区2区三区| 精品久久久精品久久久| 欧美日韩中文字幕国产精品一区二区三区 | 久久香蕉激情| 亚洲精品在线美女| 国产精品成人在线| 久热爱精品视频在线9| 一夜夜www| 91精品国产国语对白视频| 91麻豆精品激情在线观看国产 | 欧美黑人欧美精品刺激| 亚洲欧洲精品一区二区精品久久久| 亚洲av成人不卡在线观看播放网| 亚洲国产精品一区二区三区在线| 亚洲欧美精品综合久久99| 国产乱人伦免费视频| 99国产精品免费福利视频| 免费久久久久久久精品成人欧美视频| 久久香蕉精品热| 久久久国产精品麻豆| 高清在线国产一区| 欧美黑人精品巨大| 亚洲成国产人片在线观看| 亚洲av成人不卡在线观看播放网| 精品熟女少妇八av免费久了| 麻豆久久精品国产亚洲av | 三上悠亚av全集在线观看| 黄片大片在线免费观看| 国产精品亚洲一级av第二区| 麻豆av在线久日| 99国产精品一区二区三区| 露出奶头的视频| 脱女人内裤的视频| 国产精品一区二区三区四区久久 | 久久亚洲真实| 欧美 亚洲 国产 日韩一| 女生性感内裤真人,穿戴方法视频| 两个人看的免费小视频| 午夜免费观看网址| 男女高潮啪啪啪动态图| 丰满的人妻完整版| 日本三级黄在线观看| 男女午夜视频在线观看| 成人黄色视频免费在线看| 久久国产亚洲av麻豆专区| 国产精品久久久av美女十八| 国产麻豆69| 麻豆久久精品国产亚洲av | 韩国精品一区二区三区| 国产极品粉嫩免费观看在线| 男人舔女人的私密视频| 99精品久久久久人妻精品| 欧美日韩精品网址| 久久久久久久久久久久大奶| 手机成人av网站| 国产一区二区激情短视频| 国产单亲对白刺激| 午夜两性在线视频| av在线播放免费不卡| av电影中文网址| 悠悠久久av| 一二三四社区在线视频社区8| 亚洲成a人片在线一区二区| 黑人巨大精品欧美一区二区mp4| 亚洲精品国产精品久久久不卡| 久久久久久亚洲精品国产蜜桃av| 老司机在亚洲福利影院| 视频区图区小说| 国产伦一二天堂av在线观看| 精品人妻在线不人妻| 午夜精品在线福利| 欧美激情久久久久久爽电影 | 久久狼人影院| 亚洲人成电影免费在线| 欧美日韩精品网址| 丁香欧美五月| 久久久国产成人免费| 亚洲va日本ⅴa欧美va伊人久久| 国产精品久久久人人做人人爽| 在线观看日韩欧美| 黄色毛片三级朝国网站| 村上凉子中文字幕在线| 亚洲五月天丁香| 咕卡用的链子| 99re在线观看精品视频| 交换朋友夫妻互换小说| 丁香六月欧美| 亚洲av成人不卡在线观看播放网| 99国产精品一区二区蜜桃av| 亚洲在线自拍视频| 伊人久久大香线蕉亚洲五| 欧美日韩乱码在线| 国产日韩一区二区三区精品不卡| 新久久久久国产一级毛片| 18禁美女被吸乳视频| 久久久久久久久免费视频了| 久久人人精品亚洲av| 一级,二级,三级黄色视频| 精品高清国产在线一区| 免费av毛片视频| 国产麻豆69| 女同久久另类99精品国产91| 精品国产乱码久久久久久男人| 最新美女视频免费是黄的| 女人高潮潮喷娇喘18禁视频| 如日韩欧美国产精品一区二区三区| 中文亚洲av片在线观看爽| 中亚洲国语对白在线视频| 精品少妇一区二区三区视频日本电影| a在线观看视频网站| 亚洲五月色婷婷综合| 黄片大片在线免费观看| 妹子高潮喷水视频| 亚洲五月婷婷丁香| 三上悠亚av全集在线观看| 欧美日韩一级在线毛片| 大型黄色视频在线免费观看| 麻豆一二三区av精品| 国产主播在线观看一区二区| av天堂久久9| 宅男免费午夜| 我的亚洲天堂| 国产成人精品在线电影| 亚洲专区字幕在线| 日韩av在线大香蕉| 麻豆成人av在线观看| 国产精品久久视频播放| 久久婷婷成人综合色麻豆| 久9热在线精品视频| 女人被躁到高潮嗷嗷叫费观| 夜夜夜夜夜久久久久| 亚洲国产精品sss在线观看 | 悠悠久久av| 免费在线观看黄色视频的| ponron亚洲| 黄色片一级片一级黄色片| 神马国产精品三级电影在线观看 | 午夜福利欧美成人| 国产高清videossex| 国产欧美日韩一区二区三| 久久国产精品人妻蜜桃| 99在线人妻在线中文字幕| 丁香六月欧美| 免费女性裸体啪啪无遮挡网站| 免费日韩欧美在线观看| 久久久久久久久久久久大奶| 不卡一级毛片| 三级毛片av免费| 午夜影院日韩av| 高清av免费在线| 亚洲片人在线观看| 搡老乐熟女国产| 啦啦啦 在线观看视频| 啦啦啦在线免费观看视频4| 欧美日韩中文字幕国产精品一区二区三区 | 国产精品亚洲av一区麻豆| 久久国产乱子伦精品免费另类| 国产成人欧美在线观看| 久久精品aⅴ一区二区三区四区| 亚洲精品国产色婷婷电影| 99热只有精品国产| 免费一级毛片在线播放高清视频 | 日本wwww免费看| 欧美日本亚洲视频在线播放| 色综合欧美亚洲国产小说| 99国产极品粉嫩在线观看| 母亲3免费完整高清在线观看| 久久精品亚洲熟妇少妇任你| 亚洲性夜色夜夜综合| 精品久久久久久久久久免费视频 | 午夜亚洲福利在线播放| 午夜免费鲁丝| 日韩欧美一区二区三区在线观看| 婷婷丁香在线五月| 久久久久久人人人人人| 一级片免费观看大全| 午夜福利影视在线免费观看| 日本黄色日本黄色录像| 国产精品永久免费网站| 老司机深夜福利视频在线观看| 欧美日本中文国产一区发布| 久久精品国产综合久久久| 日本a在线网址| 久热爱精品视频在线9| 欧美乱妇无乱码| 久久人人97超碰香蕉20202| 国产成人欧美在线观看| 1024香蕉在线观看| 黄色怎么调成土黄色| 久久中文字幕人妻熟女| 欧美日韩精品网址| 国产色视频综合| 女同久久另类99精品国产91| 最近最新免费中文字幕在线| 精品卡一卡二卡四卡免费| 欧美激情极品国产一区二区三区| 成人影院久久| 亚洲第一av免费看| 日本一区二区免费在线视频| 精品国产超薄肉色丝袜足j| 国产精品 国内视频| 一本综合久久免费| 久久久国产精品麻豆| 黄色怎么调成土黄色| av在线播放免费不卡| 中文字幕精品免费在线观看视频| 啦啦啦 在线观看视频| 夜夜躁狠狠躁天天躁| 午夜福利一区二区在线看| 国产单亲对白刺激| 亚洲av成人不卡在线观看播放网| 成人黄色视频免费在线看| 黄网站色视频无遮挡免费观看| 美女高潮喷水抽搐中文字幕| 日本三级黄在线观看| 80岁老熟妇乱子伦牲交| 日日夜夜操网爽| 久久人妻熟女aⅴ| 亚洲五月色婷婷综合| 亚洲色图 男人天堂 中文字幕| 亚洲伊人色综图| 亚洲成国产人片在线观看| 在线永久观看黄色视频| 午夜日韩欧美国产| 欧美激情高清一区二区三区| 国产精品98久久久久久宅男小说| 一进一出抽搐动态| av在线播放免费不卡| 国产成人av激情在线播放| 香蕉丝袜av| 一级a爱视频在线免费观看| 欧美色视频一区免费| 亚洲欧美日韩无卡精品| 国产人伦9x9x在线观看| 露出奶头的视频| 亚洲欧美精品综合一区二区三区| 国产熟女午夜一区二区三区| 亚洲第一av免费看| 亚洲 欧美一区二区三区| 99精国产麻豆久久婷婷| 精品少妇一区二区三区视频日本电影| 久久久精品国产亚洲av高清涩受| 国产高清激情床上av| 国产成人免费无遮挡视频| 久久人妻熟女aⅴ| 动漫黄色视频在线观看| 成年女人毛片免费观看观看9| a级片在线免费高清观看视频| 女人爽到高潮嗷嗷叫在线视频| 老鸭窝网址在线观看| 欧美亚洲日本最大视频资源| 99久久综合精品五月天人人| 亚洲精品美女久久久久99蜜臀| 欧美人与性动交α欧美软件| 真人一进一出gif抽搐免费| 色在线成人网| 精品国产乱子伦一区二区三区| 精品国产乱码久久久久久男人| 亚洲三区欧美一区| 少妇的丰满在线观看| 欧美黑人欧美精品刺激| 中文字幕精品免费在线观看视频| tocl精华| 国产又爽黄色视频| 久久 成人 亚洲| av网站免费在线观看视频| 国产深夜福利视频在线观看| 黄色丝袜av网址大全| 国产一区二区三区在线臀色熟女 | 后天国语完整版免费观看| 丁香欧美五月| 中文字幕av电影在线播放| 久久久久久久精品吃奶| www.自偷自拍.com| 国产片内射在线| 五月开心婷婷网| 国产精品1区2区在线观看.| 国产成人av教育| 免费av毛片视频| 亚洲专区字幕在线| 中文欧美无线码| 巨乳人妻的诱惑在线观看| 亚洲 欧美一区二区三区| 欧美在线一区亚洲| 波多野结衣高清无吗| 伦理电影免费视频| 欧美日韩乱码在线| 国产无遮挡羞羞视频在线观看| 亚洲第一欧美日韩一区二区三区| 国产一区二区激情短视频| 97人妻天天添夜夜摸| 女人精品久久久久毛片|