• <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ù)解法
    国产亚洲精品久久久com| 日韩欧美 国产精品| 男女下面进入的视频免费午夜| 亚洲精品美女久久av网站| 国产成人精品久久二区二区91| 午夜激情福利司机影院| 色综合婷婷激情| 欧美一区二区精品小视频在线| 99热只有精品国产| 久久天躁狠狠躁夜夜2o2o| 国产97色在线日韩免费| 午夜日韩欧美国产| 国产精品影院久久| 欧美日韩乱码在线| 亚洲一区高清亚洲精品| 特级一级黄色大片| 亚洲成人精品中文字幕电影| 90打野战视频偷拍视频| 日本黄色片子视频| 国产毛片a区久久久久| 亚洲乱码一区二区免费版| 国产一区二区激情短视频| 男人的好看免费观看在线视频| 亚洲男人的天堂狠狠| 亚洲自偷自拍图片 自拍| 国产精品久久久久久亚洲av鲁大| 精品久久久久久久毛片微露脸| 18美女黄网站色大片免费观看| 欧美成狂野欧美在线观看| 97超级碰碰碰精品色视频在线观看| 午夜久久久久精精品| 日韩欧美在线乱码| 不卡一级毛片| www日本黄色视频网| 亚洲精品美女久久久久99蜜臀| 精品国产乱码久久久久久男人| 18禁裸乳无遮挡免费网站照片| 首页视频小说图片口味搜索| 成人永久免费在线观看视频| av在线蜜桃| 日本精品一区二区三区蜜桃| 麻豆av在线久日| tocl精华| 久久人妻av系列| 久久精品亚洲精品国产色婷小说| 亚洲自拍偷在线| 色吧在线观看| 久久亚洲真实| 最新在线观看一区二区三区| 日日夜夜操网爽| 国内少妇人妻偷人精品xxx网站 | 成人三级做爰电影| 久久久久国产精品人妻aⅴ院| 12—13女人毛片做爰片一| 精品不卡国产一区二区三区| 村上凉子中文字幕在线| 亚洲真实伦在线观看| 午夜免费观看网址| 久久久久久久精品吃奶| 亚洲欧美精品综合久久99| 特级一级黄色大片| 午夜福利18| 国产伦人伦偷精品视频| 美女高潮的动态| 看免费av毛片| 在线观看66精品国产| 他把我摸到了高潮在线观看| 精品久久久久久久末码| 人妻久久中文字幕网| 亚洲av中文字字幕乱码综合| 老司机在亚洲福利影院| 精品国产三级普通话版| 欧美zozozo另类| 国产成人一区二区三区免费视频网站| 天堂网av新在线| 麻豆久久精品国产亚洲av| 无遮挡黄片免费观看| 日韩三级视频一区二区三区| 精品欧美国产一区二区三| 欧美中文综合在线视频| h日本视频在线播放| 欧美日本亚洲视频在线播放| 男人舔女人的私密视频| 日韩欧美国产一区二区入口| 国产91精品成人一区二区三区| 俄罗斯特黄特色一大片| 精品久久久久久久末码| 在线观看免费午夜福利视频| 51午夜福利影视在线观看| 欧美精品啪啪一区二区三区| 国产激情偷乱视频一区二区| 精品久久久久久久久久免费视频| 国产欧美日韩精品一区二区| 丁香六月欧美| 亚洲av成人一区二区三| 亚洲欧美日韩高清在线视频| 国产97色在线日韩免费| 久久草成人影院| 久久午夜综合久久蜜桃| 啦啦啦韩国在线观看视频| 国产午夜精品论理片| 国内精品一区二区在线观看| 网址你懂的国产日韩在线| 欧美日韩瑟瑟在线播放| 亚洲 欧美一区二区三区| 国产成人精品久久二区二区免费| 亚洲中文日韩欧美视频| 国产高清videossex| 99久国产av精品| 精品久久久久久久久久久久久| 国产成人欧美在线观看| 日本 欧美在线| 国产伦精品一区二区三区视频9 | 国产一区二区在线观看日韩 | 亚洲无线在线观看| 国产一区二区在线观看日韩 | 日韩欧美 国产精品| 国产成人av激情在线播放| 精品无人区乱码1区二区| 亚洲aⅴ乱码一区二区在线播放| 午夜a级毛片| 黄色 视频免费看| 亚洲va日本ⅴa欧美va伊人久久| 国产 一区 欧美 日韩| 操出白浆在线播放| 亚洲欧美日韩卡通动漫| 亚洲熟妇熟女久久| 国产成人系列免费观看| av在线蜜桃| 国产免费男女视频| 午夜福利成人在线免费观看| 性色avwww在线观看| 亚洲性夜色夜夜综合| or卡值多少钱| 黄色视频,在线免费观看| 啪啪无遮挡十八禁网站| 丁香六月欧美| www日本在线高清视频| 国产99白浆流出| www日本黄色视频网| 狠狠狠狠99中文字幕| 国产真实乱freesex| 国产久久久一区二区三区| 白带黄色成豆腐渣| 美女免费视频网站| 久久久国产精品麻豆| 男人和女人高潮做爰伦理| 天堂网av新在线| 12—13女人毛片做爰片一| 女生性感内裤真人,穿戴方法视频| 后天国语完整版免费观看| 一级毛片精品| 日本五十路高清| 51午夜福利影视在线观看| 特大巨黑吊av在线直播| 亚洲国产欧美一区二区综合| 欧美一区二区国产精品久久精品| 国产美女午夜福利| 女人被狂操c到高潮| 99久国产av精品| 一级毛片高清免费大全| 好男人电影高清在线观看| 亚洲av日韩精品久久久久久密| 天天添夜夜摸| 后天国语完整版免费观看| 亚洲av电影不卡..在线观看| 国产av麻豆久久久久久久| 又紧又爽又黄一区二区| 国产精品av视频在线免费观看| 色精品久久人妻99蜜桃| 国产在线精品亚洲第一网站| 日韩欧美 国产精品| 久久久久久大精品| 日韩成人在线观看一区二区三区| 18禁美女被吸乳视频| 日本三级黄在线观看| 免费在线观看日本一区| 久久久久久人人人人人| 久久久久久久久久黄片| 少妇裸体淫交视频免费看高清| 可以在线观看的亚洲视频| 99热这里只有是精品50| 不卡av一区二区三区| 亚洲国产看品久久| 色播亚洲综合网| 亚洲男人的天堂狠狠| 日韩中文字幕欧美一区二区| 国产精品亚洲一级av第二区| 欧美午夜高清在线| 18美女黄网站色大片免费观看| 国产又黄又爽又无遮挡在线| 精品久久久久久久久久久久久| 成人永久免费在线观看视频| 麻豆成人av在线观看| 美女大奶头视频| 俄罗斯特黄特色一大片| 亚洲精品456在线播放app | 国产久久久一区二区三区| 久久精品国产综合久久久| 亚洲色图av天堂| 美女午夜性视频免费| 久久久成人免费电影| 香蕉av资源在线| 亚洲av成人av| av视频在线观看入口| 18禁黄网站禁片午夜丰满| 国产av一区在线观看免费| 全区人妻精品视频| 在线永久观看黄色视频| 久久久久九九精品影院| 人妻久久中文字幕网| 中文字幕久久专区| 亚洲最大成人中文| 一二三四在线观看免费中文在| 色综合婷婷激情| 啦啦啦观看免费观看视频高清| 国产探花在线观看一区二区| 国产亚洲精品一区二区www| 国产精品一区二区三区四区免费观看 | 成人18禁在线播放| 中文字幕久久专区| 在线国产一区二区在线| 国产黄片美女视频| 婷婷精品国产亚洲av在线| 精品国产三级普通话版| 亚洲,欧美精品.| 波多野结衣高清作品| 久9热在线精品视频| 日日摸夜夜添夜夜添小说| 国产一级毛片七仙女欲春2| 男女床上黄色一级片免费看| 很黄的视频免费| 中文资源天堂在线| 丁香欧美五月| 又大又爽又粗| 久久久久国产精品人妻aⅴ院| 精品国产三级普通话版| 国产亚洲欧美98| 老熟妇乱子伦视频在线观看| 国内毛片毛片毛片毛片毛片| 亚洲欧洲精品一区二区精品久久久| 夜夜看夜夜爽夜夜摸| 亚洲狠狠婷婷综合久久图片| 在线播放国产精品三级| 大型黄色视频在线免费观看| 久久性视频一级片| netflix在线观看网站| 亚洲美女黄片视频| 非洲黑人性xxxx精品又粗又长| 91在线观看av| 久久婷婷人人爽人人干人人爱| 久99久视频精品免费| 亚洲自偷自拍图片 自拍| 亚洲美女视频黄频| 欧美xxxx黑人xx丫x性爽| 99久久99久久久精品蜜桃| 美女免费视频网站| 99在线视频只有这里精品首页| 国产精品一区二区精品视频观看| 久久久久久久午夜电影| 亚洲国产精品久久男人天堂| av中文乱码字幕在线| cao死你这个sao货| 窝窝影院91人妻| 999久久久国产精品视频| 99国产精品一区二区蜜桃av| 亚洲在线观看片| 日本免费一区二区三区高清不卡| 欧美黑人欧美精品刺激| 国产麻豆成人av免费视频| 久久久久免费精品人妻一区二区| 在线永久观看黄色视频| aaaaa片日本免费| 高潮久久久久久久久久久不卡| 亚洲专区中文字幕在线| 偷拍熟女少妇极品色| 在线观看午夜福利视频| 可以在线观看毛片的网站| 国产男靠女视频免费网站| 日日干狠狠操夜夜爽| 国内久久婷婷六月综合欲色啪| 午夜两性在线视频| or卡值多少钱| 欧美大码av| 国产爱豆传媒在线观看| 啪啪无遮挡十八禁网站| 日日夜夜操网爽| 黄片小视频在线播放| 午夜激情欧美在线| 手机成人av网站| 美女黄网站色视频| 精品无人区乱码1区二区| 天天一区二区日本电影三级| 婷婷六月久久综合丁香| 亚洲成人久久爱视频| 国产伦在线观看视频一区| 黄色女人牲交| 国产精品美女特级片免费视频播放器 | 久久久久国产精品人妻aⅴ院| 久久中文看片网| 免费人成视频x8x8入口观看| 无限看片的www在线观看| 亚洲精品粉嫩美女一区| 窝窝影院91人妻| 97超视频在线观看视频| 女人被狂操c到高潮| 国产精品乱码一区二三区的特点| 亚洲,欧美精品.| 欧美日韩一级在线毛片| 天堂av国产一区二区熟女人妻| 日本五十路高清| 夜夜爽天天搞| 波多野结衣巨乳人妻| 成人国产一区最新在线观看| 国产成人精品久久二区二区91| 亚洲国产日韩欧美精品在线观看 | 久久精品夜夜夜夜夜久久蜜豆| 免费观看的影片在线观看| 久久久久久久久中文| 亚洲国产精品sss在线观看| 亚洲成人久久爱视频| 欧美激情久久久久久爽电影| 国产一区二区激情短视频| 婷婷丁香在线五月| 国产99白浆流出| 俄罗斯特黄特色一大片| 天堂影院成人在线观看| 国产精品影院久久| 嫁个100分男人电影在线观看| 亚洲国产欧美网| 午夜a级毛片| 一级作爱视频免费观看| 国产毛片a区久久久久| 国产黄色小视频在线观看| 国产高清三级在线| 巨乳人妻的诱惑在线观看| 欧美色视频一区免费| 熟妇人妻久久中文字幕3abv| 99精品欧美一区二区三区四区| 亚洲人成伊人成综合网2020| 全区人妻精品视频| 精品免费久久久久久久清纯| 色尼玛亚洲综合影院| 啦啦啦免费观看视频1| 久9热在线精品视频| 免费看光身美女| 美女高潮喷水抽搐中文字幕| 欧美丝袜亚洲另类 | 99视频精品全部免费 在线 | 91老司机精品| www日本黄色视频网| 久久性视频一级片| 手机成人av网站| a级毛片a级免费在线| 亚洲欧洲精品一区二区精品久久久| av天堂在线播放| 一本久久中文字幕| 色老头精品视频在线观看| 欧美一区二区国产精品久久精品| 中出人妻视频一区二区| 男女下面进入的视频免费午夜| 熟女人妻精品中文字幕| 亚洲人成网站在线播放欧美日韩| 国产高清有码在线观看视频| 欧美高清成人免费视频www| 日韩成人在线观看一区二区三区| 色综合欧美亚洲国产小说| 国产乱人视频| 成人无遮挡网站| 国产亚洲精品一区二区www| 深夜精品福利| 成人18禁在线播放| 一级a爱片免费观看的视频| 熟女人妻精品中文字幕| 国产午夜福利久久久久久| 国产亚洲精品久久久com| 又黄又爽又免费观看的视频| 嫩草影院精品99| 在线国产一区二区在线| 亚洲,欧美精品.| 99热这里只有是精品50| 97人妻精品一区二区三区麻豆| 亚洲人成电影免费在线| 婷婷精品国产亚洲av| 夜夜爽天天搞| 亚洲一区二区三区不卡视频| 黄频高清免费视频| 黑人操中国人逼视频| 日本精品一区二区三区蜜桃| 日韩大尺度精品在线看网址| 91麻豆精品激情在线观看国产| 美女高潮的动态| 91av网一区二区| 午夜激情福利司机影院| 色哟哟哟哟哟哟| 一级作爱视频免费观看| 国产99白浆流出| 制服丝袜大香蕉在线| 国产成人精品久久二区二区免费| 九色成人免费人妻av| 国产97色在线日韩免费| 麻豆久久精品国产亚洲av| 性色avwww在线观看| 成人欧美大片| 国产成人福利小说| 久久热在线av| 成年女人永久免费观看视频| 天堂av国产一区二区熟女人妻| 久久中文看片网| 夜夜躁狠狠躁天天躁| 国产高清激情床上av| 免费在线观看亚洲国产| 免费看a级黄色片| 精品熟女少妇八av免费久了| 欧美日本亚洲视频在线播放| 热99re8久久精品国产| 日本黄大片高清| 国产美女午夜福利| 午夜视频精品福利| 真人一进一出gif抽搐免费| 身体一侧抽搐| 桃红色精品国产亚洲av| 亚洲专区中文字幕在线| 一级作爱视频免费观看| av中文乱码字幕在线| 欧美日韩一级在线毛片| 国产精品国产高清国产av| 18禁美女被吸乳视频| 黄色丝袜av网址大全| 88av欧美| 国产精品一区二区免费欧美| 五月玫瑰六月丁香| 两个人看的免费小视频| 在线永久观看黄色视频| 国内精品久久久久精免费| 免费电影在线观看免费观看| 美女扒开内裤让男人捅视频| 黄片小视频在线播放| 成人三级做爰电影| 亚洲国产精品sss在线观看| 欧美中文综合在线视频| 国内揄拍国产精品人妻在线| bbb黄色大片| 成人av一区二区三区在线看| 麻豆国产av国片精品| 日韩人妻高清精品专区| 女警被强在线播放| 三级国产精品欧美在线观看 | 亚洲欧美日韩东京热| 国产精品久久久久久人妻精品电影| 国内毛片毛片毛片毛片毛片| 少妇的逼水好多| 在线观看美女被高潮喷水网站 | 叶爱在线成人免费视频播放| 午夜精品一区二区三区免费看| 国产精品一区二区三区四区久久| 天堂网av新在线| or卡值多少钱| 91在线观看av| 亚洲激情在线av| 一进一出抽搐gif免费好疼| 亚洲国产欧美一区二区综合| 国产一级毛片七仙女欲春2| 久9热在线精品视频| 国产男靠女视频免费网站| x7x7x7水蜜桃| 我的老师免费观看完整版| 最近在线观看免费完整版| 国产99白浆流出| 少妇裸体淫交视频免费看高清| 精品福利观看| 精品久久久久久久末码| 国产精品1区2区在线观看.| 欧美成人性av电影在线观看| 超碰成人久久| www.www免费av| 黄色成人免费大全| 亚洲av熟女| 99国产精品一区二区三区| 91av网站免费观看| 啦啦啦韩国在线观看视频| 久久中文字幕一级| 国产亚洲精品综合一区在线观看| 听说在线观看完整版免费高清| 黑人操中国人逼视频| 女生性感内裤真人,穿戴方法视频| 精品国内亚洲2022精品成人| 久久这里只有精品中国| 国产精品亚洲av一区麻豆| 欧美av亚洲av综合av国产av| 成人鲁丝片一二三区免费| 日本与韩国留学比较| 国产69精品久久久久777片 | 欧美日韩亚洲国产一区二区在线观看| www日本黄色视频网| 久久热在线av| 国产欧美日韩精品亚洲av| 国产精品亚洲美女久久久| 制服丝袜大香蕉在线| 国产亚洲av高清不卡| 一边摸一边抽搐一进一小说| 精品一区二区三区视频在线观看免费| 久久天躁狠狠躁夜夜2o2o| 久久久久久九九精品二区国产| 麻豆成人午夜福利视频| 国产精品久久久av美女十八| aaaaa片日本免费| 精品人妻1区二区| 国产日本99.免费观看| 又黄又爽又免费观看的视频| a级毛片a级免费在线| 曰老女人黄片| 午夜福利在线观看吧| 又紧又爽又黄一区二区| 一本久久中文字幕| 国产精品久久久久久久电影 | 国产精品久久视频播放| 精品不卡国产一区二区三区| 真人一进一出gif抽搐免费| 午夜两性在线视频| 在线观看日韩欧美| 午夜福利高清视频| 可以在线观看的亚洲视频| 一进一出抽搐gif免费好疼| 久久久久国产一级毛片高清牌| 国产精品久久久av美女十八| 欧美另类亚洲清纯唯美| 99久久精品热视频| 色尼玛亚洲综合影院| 亚洲av日韩精品久久久久久密| 真人做人爱边吃奶动态| 久久中文字幕人妻熟女| 国产亚洲精品一区二区www| 桃红色精品国产亚洲av| 亚洲一区二区三区色噜噜| 亚洲国产色片| 国产精品影院久久| aaaaa片日本免费| 国产私拍福利视频在线观看| 久久久精品大字幕| 琪琪午夜伦伦电影理论片6080| 欧美最黄视频在线播放免费| 久久精品综合一区二区三区| 香蕉久久夜色| 亚洲欧美一区二区三区黑人| 欧美精品啪啪一区二区三区| 成人性生交大片免费视频hd| 日韩欧美国产一区二区入口| 国产精品亚洲av一区麻豆| 国产不卡一卡二| 婷婷丁香在线五月| 国产黄a三级三级三级人| 亚洲中文日韩欧美视频| 日韩欧美免费精品| 国内精品一区二区在线观看| 亚洲成人久久爱视频| 2021天堂中文幕一二区在线观| 国产成人欧美在线观看| 婷婷六月久久综合丁香| 黑人操中国人逼视频| 在线观看美女被高潮喷水网站 | 免费高清视频大片| 亚洲国产精品合色在线| 日本 欧美在线| 精品人妻1区二区| xxx96com| 久久香蕉国产精品| 久久国产精品影院| 99在线视频只有这里精品首页| 最新中文字幕久久久久 | 国产一区二区三区视频了| 国产单亲对白刺激| 亚洲国产看品久久| 久久伊人香网站| 俄罗斯特黄特色一大片| 手机成人av网站| 香蕉国产在线看| 日韩欧美国产一区二区入口| www国产在线视频色| 91麻豆精品激情在线观看国产| 国产探花在线观看一区二区| 嫩草影院入口| 午夜影院日韩av| 国产探花在线观看一区二区| 老鸭窝网址在线观看| 嫩草影视91久久| 久久国产精品人妻蜜桃| av视频在线观看入口| 又粗又爽又猛毛片免费看| 亚洲天堂国产精品一区在线| 亚洲九九香蕉| 999久久久精品免费观看国产| 熟女少妇亚洲综合色aaa.| 欧美一区二区国产精品久久精品| 午夜两性在线视频| 色综合欧美亚洲国产小说| 久久久久久大精品| 精品一区二区三区视频在线观看免费| 午夜福利成人在线免费观看| 综合色av麻豆| 国产99白浆流出| 国产成人啪精品午夜网站| 久久精品国产亚洲av香蕉五月| 成人欧美大片| 18禁美女被吸乳视频| 搡老熟女国产l中国老女人| 国产淫片久久久久久久久 | 免费人成视频x8x8入口观看| 不卡av一区二区三区| 狂野欧美白嫩少妇大欣赏| 啦啦啦免费观看视频1| 中文在线观看免费www的网站| 国产高清有码在线观看视频| 一个人看的www免费观看视频| 亚洲五月天丁香|