亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種確定聯(lián)合反演中相對權(quán)比的兩步法

        2012-01-04 08:00:38王樂洋許才軍張朝玉
        測繪學報 2012年1期
        關(guān)鍵詞:方差反演觀測

        王樂洋,許才軍,張朝玉

        1.東華理工大學測繪工程學院,江西南昌330013;2.武漢大學測繪學院,湖北武漢430079

        一種確定聯(lián)合反演中相對權(quán)比的兩步法

        王樂洋1,許才軍2,張朝玉2

        1.東華理工大學測繪工程學院,江西南昌330013;2.武漢大學測繪學院,湖北武漢430079

        提出確定聯(lián)合反演中相對權(quán)比的兩步法。第1步是不考慮相對權(quán)比,建立兩類或多類數(shù)據(jù)的聯(lián)合反演模型,通過赫爾墨特方差分量估計使得兩類或多類觀測數(shù)據(jù)的單位權(quán)方差相等(或單位權(quán)方差的比值接近1),從而理論上確定觀測數(shù)據(jù)合理的權(quán)陣,同時對不同種類數(shù)據(jù)進行歸一化;第2步是在獲得合理的觀測數(shù)據(jù)權(quán)陣之后,建立顧及相對權(quán)比的聯(lián)合反演模型,以目標函數(shù)值為最小來確定權(quán)比例因子。通過模擬數(shù)據(jù),設(shè)計多個反演方案,證明兩步法的有效性和可行性。

        聯(lián)合反演;相對權(quán)比;兩步法;赫爾墨特方差分量估計;單位權(quán)方差

        1 引 言

        大地測量地球物理聯(lián)合反演可以融合多種數(shù)據(jù)(大地測量數(shù)據(jù)、地震波數(shù)據(jù)、地質(zhì)數(shù)據(jù)以及地球物理數(shù)據(jù)等),研究地球動力學,反演活動斷層、塊體的運動情況,探討地殼運動與地震的關(guān)系,是大地測量學深入地學研究領(lǐng)域的一個重要的手段。聯(lián)合反演應(yīng)用的一個關(guān)鍵問題是相對權(quán)比的確定,它體現(xiàn)了各類資料在聯(lián)合反演中貢獻的大小,一個成功的聯(lián)合反演必須有一個合理的相對權(quán)比,因此相對權(quán)比的研究是當今聯(lián)合反演研究的一個熱點問題,得到了廣泛關(guān)注[1-14]。目前相對權(quán)比確定的方法主要有4種:將相對權(quán)比當做未知數(shù)與模型參數(shù)一同反演求解;根據(jù)先驗信息(驗前方差)確定相對權(quán)比;相對權(quán)比人為定為0.5;赫爾墨特方差分量估計法[2-4,7,9,13]。本文針對這些方法的優(yōu)缺點,提出確定線性聯(lián)合反演中相對權(quán)比的兩步法,最后通過模擬計算顯示方法的有效性和可行性。

        2 確定聯(lián)合反演中相對權(quán)比的兩步法

        對于線性反演模型,假設(shè)有兩種不同類型的觀測數(shù)據(jù),聯(lián)合反演的目標函數(shù)為

        將式(1)對^X求導(dǎo)有

        轉(zhuǎn)置后得

        將V1和V2代入上式,整理得

        兩步法確定線性聯(lián)合反演模型中相對權(quán)比的具體步驟如下。

        ① 先不考慮相對權(quán)比λ,聯(lián)合反演的目標函數(shù)為

        法方程為

        所以,式(5)變?yōu)?/p>

        確定兩類觀測資料權(quán)的初值P1和P2,選取時根據(jù)前面計算得到的結(jié)果作為初值,如利用GPS和地震矩張量聯(lián)合反演地殼運動速度場時,P1可選為GPS站速度場的先驗權(quán)陣,P2可選為平均應(yīng)變率的先驗權(quán)陣。一般來說兩類觀測數(shù)據(jù)的權(quán)P1和P2是可能不準確的,也就是說P1和P2所對應(yīng)的單位權(quán)方差不相等聯(lián)合反演中利用赫爾墨特方差分量估計的前提是:兩類或多類觀測值之間是相互獨立的,不獨立時必須利用方差-協(xié)方差分量估計來確定單位權(quán)方差與單位權(quán)協(xié)方差,然后獲得合理的觀測值權(quán)陣。

        ③ 進行第一次方差分量估計[15]

        式中

        第1步所用的方程為

        目標函數(shù)是

        (2)第2步,根據(jù)第1步得到的兩類觀測資料權(quán)陣的估值^P1和^P2,在式(10)下求取相對權(quán)比和模型參數(shù),即

        具體步驟是:選定一個步長,讓權(quán)比例因子τ遍歷整個取值區(qū)間0<τ<1,每個τ對應(yīng)于一個參數(shù)解和,把其中的最小值對應(yīng)的參數(shù)向量和此時的權(quán)比例因子作為最終結(jié)果。

        由式(9)和式(10)得λ與τ的關(guān)系為

        式中,j=1,2,3,…為兩步法第1步的迭代次數(shù)。

        可以將兩步法簡單概括如下:第1步是不考慮相對權(quán)比,建立兩類或多類數(shù)據(jù)的聯(lián)合反演模型,通過赫爾墨特方差分量估計使得兩類或多類觀測數(shù)據(jù)的單位權(quán)方差相等(或單位權(quán)方差的比值接近1),從理論上確定觀測數(shù)據(jù)合理的權(quán)陣,同時對不同種類數(shù)據(jù)進行歸一化;第2步是在獲得合理的觀測數(shù)據(jù)權(quán)陣之后,建立顧及相對權(quán)比的聯(lián)合反演模型,將相對權(quán)比和模型參數(shù)當做未知數(shù)同時進行反演,以目標函數(shù)值為最小來確定權(quán)比例因子。兩步法是一種概括反演方法,赫爾墨特方差分量估計法、將相對權(quán)比當做未知數(shù)同模型參數(shù)一同反演的方法以及單一反演方法都是兩步法的特例;其優(yōu)點是可以獲得觀測數(shù)據(jù)合理權(quán)陣,同時使得顧及相對權(quán)比的聯(lián)合反演目標函數(shù)取得最小值,并對各類觀測數(shù)據(jù)進行歸一化。

        當不同源數(shù)據(jù)出現(xiàn)量綱不一致情況時,首先需要對多源數(shù)據(jù)進行無量綱化處理,常用的方法為或其中為觀測數(shù)據(jù),F(xiàn)i(m)為各自模型擬合的觀測值,為觀測數(shù)據(jù)的標準差,下標i表示第i類數(shù)據(jù)[10]。

        3 算例及分析

        3.1 算 例

        在假定介質(zhì)均勻且鄰近測點間應(yīng)變均勻的情況下,可以建立鄰近點間相對形變量與地殼應(yīng)變張量的線性關(guān)系[16]。設(shè)監(jiān)測網(wǎng)第j點兩期縱坐標位移為uj,橫坐標位移為vj,網(wǎng)中共有m點,則由坐標位移反演應(yīng)變參數(shù)的公式(二維)為[17]

        式中,xj、yj(j=1,2,…,m)為第j點近似坐標;ω為其轉(zhuǎn)動量;εx、εy、γxy為應(yīng)變狀態(tài)參數(shù),且γxy=2εxy。

        利用該應(yīng)變參數(shù)反演模型模擬生成了反演應(yīng)變參數(shù)的GPS坐標位移數(shù)據(jù)和邊角網(wǎng)平差后的坐標位移數(shù)據(jù)(u1、v1)和(u2、v2),并對兩類觀測資料分別施加了σ2=0.5mm2和ˉσ2=1.0mm2的隨機噪聲后的數(shù)據(jù)為(u′1、v′1)和(u′2、v′2)。具體見表1和表2,相應(yīng)的坐標位移協(xié)因數(shù)陣數(shù)據(jù)分別見表3和表4。

        表1 GPS坐標數(shù)據(jù)及坐標位移觀測值數(shù)據(jù)Tab.1 Coordinate and displacement observation of GPS stations

        表2 邊角網(wǎng)的坐標數(shù)據(jù)及坐標位移觀測值數(shù)據(jù)Tab.2 Coordinate and displacement observation of side-angle network

        為了比較本文提出的確定線性聯(lián)合反演模型相對權(quán)比的兩步法與其他方法之間的優(yōu)缺點,利用模擬數(shù)據(jù)進行了計算,分以下10個方案進行,λ表示GPS坐標位移在聯(lián)合反演中所占的比重。方案1:GPS坐標位移的單一反演;方案2:邊角網(wǎng)平差坐標位移的單一反演;方案3:λ=0.5;方案4:方案5:以把λ當做未知數(shù)同模型參數(shù)一同進行反演;方案6:赫爾墨特方差分量估計法[4],迭代終止的條件為兩類數(shù)據(jù)單位權(quán)方差估值之差小于等于0.000 1;為了比較兩步法中第2步的目標函數(shù)ˉΦ(V1,V2)(式(10)所示)的不同對確定相對權(quán)比的影響,分別計算了以以及為目標函數(shù)確定相對權(quán)比的兩步法,分別稱為方案7、方案8、方案9和方案10,迭代計算時τ步長取為0.000 1;為便于閱讀將各方案的所用方法列于表5中。

        各方案應(yīng)變參數(shù)反演結(jié)果見表6。

        表3 GPS坐標位移協(xié)因數(shù)陣數(shù)據(jù)Tab.3 Cofactor matrix of GPS coordinate displacement

        表4 邊角網(wǎng)平差后坐標位移的協(xié)因數(shù)陣數(shù)據(jù)Tab.4 Cofactor matrix of side-angle network coordinate displacement after adjustment

        表5 各方案所用方法列Tab.5 The method of each project

        表6 應(yīng)變參數(shù)反演結(jié)果Tab.6 Inversion results of model parameters

        由于方案9與方案2等價,因此,表6中沒有列出方案9的計算結(jié)果,由方案7、方案8、方案9和方案10的兩步法計算出來的相對權(quán)比τ分別為0.500 0、0.212 7、0(不加0<τ<1的限制,允許取值0和1)和0.308 2。根據(jù)式(11)將對應(yīng)于的相對權(quán)比τ轉(zhuǎn)換成對應(yīng)于式(1)中Pi(i=1,2)的相對權(quán)比λ分別為0.955 1、0.851 7、0和0.904 5(如表6所示);‖ΔX‖表示模型參數(shù)反演結(jié)果與真值差值的范數(shù);R表示參數(shù)的真值。然后,本文繪制了方案5、方案7、方案8、方案9和方案10的目標函數(shù)與相對權(quán)比之間的關(guān)系圖,如圖1所示。

        3.2 結(jié)果分析

        (1)對于反演結(jié)果的評價一般要從內(nèi)符合精度和外符合精度兩個方面進行綜合考慮。本文參數(shù)估值的中誤差是反演結(jié)果的內(nèi)符合精度指標,反演結(jié)果與真值的差值范數(shù)‖ΔX‖作為外符合精度指標。整體上從內(nèi)符合精度來看,方案6和方案7的結(jié)果最優(yōu),而且二者是等價的。從外符合精度來看,方案4的差值范數(shù)最小,這與理論相符,因為給定的權(quán)比即為實際權(quán)比。然而這種方法不具操作性,因為在反演之前兩類數(shù)據(jù)的驗前單位權(quán)方差是很難準確知道的;方案3雖然外符合精度較高,但是內(nèi)符合精度稍差;方案5的結(jié)果與方案1幾乎相同,若不加0<λ<1的限制允許取值0和1,則方案5與方案1是等價的;方案8具有較高的外符合和內(nèi)符合精度;方案10具有較高的外符合精度,但內(nèi)符合精度較差;方案9與方案2等價,為邊角網(wǎng)坐標位移的單一反演。綜合起來看,當觀測值不含粗差時,兩步法(除方案10外)可以得到較好的反演結(jié)果。

        (2)驗前單位權(quán)方差法(方案4)的主要缺點是反演結(jié)果與給予的驗前(初始)單位權(quán)方差有著很大的關(guān)系。驗前單位權(quán)方差給的準確,則可以得到較好的結(jié)果(如本文的算例),若初始單位權(quán)方差不準確,就會得到不好的結(jié)果[2,4]。赫爾墨特方差分量估計法(方案6)的主要優(yōu)點是得到的結(jié)果具有最優(yōu)無偏的特性,不會受到驗前單位權(quán)方差不準確的限制。兩步法是一種概括反演方法,赫爾墨特方差分量估計法(方案6)以及單一反演方法(方案2)等都是兩步法的特例。兩步法通過目標函數(shù)的選擇統(tǒng)一了這些常用的方法。兩步法同樣不會受到驗前單位權(quán)方差不準確的限制,可以獲得觀測數(shù)據(jù)合理權(quán)陣,同時使得顧及相對權(quán)比的聯(lián)合反演目標函數(shù)取得最小值,并對各類觀測數(shù)據(jù)進行了歸一化。由圖1可以看出,兩步法的第2步采用不同的目標函數(shù)確定相對權(quán)比的結(jié)果是不同的。兩步法目標函數(shù)的選取方法為,觀測值僅含偶然誤差時,若要獲得最優(yōu)無偏的反演結(jié)果,則選為目標函數(shù);相對權(quán)比和模型參數(shù)都當做未知數(shù)同時計算時,可選為目標函數(shù);相對權(quán)比和模型參數(shù)都當做未知數(shù)同時計算時,不可選作為目標函數(shù),否則變成單一反演。因此,目標函數(shù)的選擇可以根據(jù)實際情況確定。

        (3)由殘差圖1(a)(方案5)和圖1(d)(方案4)可以看出與相對權(quán)比λ以及與相對權(quán)比τ之間是成遞減或遞增的關(guān)系,因此當以此類目標函數(shù)作為確定相對權(quán)比的標準,且將其當做未知數(shù)與模型參數(shù)一同進行求解時,得到相對權(quán)比的值為1或者0,如文獻[4]中的λ=1。

        (4)聯(lián)合反演中相對權(quán)比的大小雖然是由所用的方法確定的,但并不能以得到的相對權(quán)比的大小作為評價確定相對權(quán)比方法好壞的標準。無論采用何種方法確定相對權(quán)比,一般精度較高、數(shù)目較多的一類觀測值應(yīng)該在聯(lián)合反演中起主要作用,該觀測值的相對權(quán)比也較大,如本文算例中的GPS觀測數(shù)據(jù)。

        相對權(quán)比的確定是聯(lián)合反演中的關(guān)鍵問題。本文提出的確定聯(lián)合反演中相對權(quán)比的兩步法是一種概括反演方法,其優(yōu)點是不受驗前單位權(quán)方差的限制,可以獲得觀測數(shù)據(jù)的合理權(quán)陣,同時使得顧及相對權(quán)比的聯(lián)合反演目標函數(shù)取得最小值,并對各類觀測數(shù)據(jù)進行了歸一化。模擬數(shù)據(jù)的計算,顯示了在處理線性聯(lián)合反演問題時該方法是非常有效的。但是,兩步法還需要在實際應(yīng)用中進一步驗證和改進。

        [1] XU Caijun,SHEN Wenbin,CHAO Dingbo.Geophysical Geodesy Principles and Methods[M].Wuhan:Wuhan University Press,2006.(許才軍,申文斌,晁定波.地球物理大地測量學原理與方法[M].武漢:武漢大學出版社,2006.)

        [2] DING Kaihua.Research on Stochastic Modeling in Geodetic Joint Inversion[D].Wuhan:Wuhan University,2006.(丁開華.大地測量聯(lián)合反演中的隨機模型研究[D].武漢:武漢大學,2006.)

        [3] FU Yuning.Joint Inversion of Movement Parameters of Main Faults and Blocks in Sichuan and Yunnan Area Using GPS and Gravity Data[D].Wuhan:Wuhan University,2007.(富宇寧.利用GPS和重力數(shù)據(jù)聯(lián)合反演川滇地區(qū)主要斷層、塊體的運動參數(shù)[D].武漢:武漢大學,2007.)

        [4] XU Caijun,DING Kaihua,CAI Jianqing,et al.Methods of Determining Weight Scaling Factors for Geodeticgeophysical Joint Inversion[J].Journal of Geodynamics,2009,47(1):39-46.

        [5] DU Zhixing,OU Jikun,JIN Fengxiang,et al.Optimization Inversion of the Relative Weight Ratio in the Joint Inversion Models[J].Acta Geodaetica et Cartographica Sinica.2003,32(1):15-19.(獨知行,歐吉坤,靳奉祥,等.聯(lián)合反演模型中相對權(quán)比的優(yōu)化反演[J].測繪學報,2003,32(1):15-19.)

        [6] XU Caijun.Progress of Joint Inversion on Geodesy and Geophysics[J].Geomatics and Information Science of Wuhan University,2001,26(6):555-561.(許才軍.大地測量聯(lián)合反演理論和方法研究進展[J].武漢大學學報:信息科學版,2001,26(6):555-561.)

        [7] FU Yuning,XU Caijun.On Relative Weight Problem in Joint Inversion Using Levelling and Gravity Data[J].Journal of Geodesy and Geodynamics,2007,27(2):68-74.(富宇寧,許才軍.水平和重力資料聯(lián)合反演中權(quán)問題的研究[J].大地測量與地球動力學,2007,27(2):68-74.)

        [8] LI Shuang.Research on Models and Algorithms of Geodesy and Geophysics Joint Inversion[D].Wuhan:Wuhan University,2005.(李爽.大地測量聯(lián)合反演的模式及算法研究[D].武漢:武漢大學,2005.)

        [9] DING Kaihua,XU Cajjun.Current Crustal Strain Field in the Sichuan-Yunnan Area by Joint Inversion of GPS and Seismic Moment Tensor[J].Geomatics and Information Science of Wuhan University,2009,34(3):265-268.(丁開華,許才軍.川滇地區(qū)地殼應(yīng)變場的GPS與地震矩張量聯(lián)合反演研究[J].武漢大學學報:信息科學版,2009,34(3):265-268.)

        [10] ZHAO Shaorong.Joint Inversion of Observed Gravity and GPS Baseline Changes for the Detection of the Active Fault Segment at the Red River Fault Zone[J].Geophysical Journal International,1995,122(1):70-88.

        [11] XU Caijun,WANG Leyang.Progress of Joint Inversion of Geodetic and Seismological Data for Seismic Source Rupture Process[J].Geomatics and Information Science of Wuhan University,2010,35(4):457-462.(許才軍,王樂洋.大地測量和地震資料聯(lián)合反演地震震源破裂過程研究進展[J].武漢大學學報:信息科學版,2010,35(4):457-462.)

        [12] WANG Leyang,XU Caijun.Comparative Research on Equality Constraint Inversion and Joint Inversion[J].Journal of Geodesy and Geodynamics,2009,29(1):74-78.(王樂洋,許才軍.等式約束反演與聯(lián)合反演的對比研究[J].大地測量與地球動力學,2009,29(1):74-78.)

        [13] XU Caijun,LIU Yang,WEN Yangmao,et al.Coseismic Slip Distribution of the 2008Mw 7.9Wenchuan Earthquake from Joint Inversion of GPS and InSAR Data[J].Bulletin of the Seismological Society of America,2010,100(5B):2736-2749.

        [14] LIU Chongbing,NING Jinsheng,ZHANG Yushen.Method for Joint Inversion of Seismic Surface Waves and Gravity Data for Three Dimensional Density Structure[J].Acta Geodaetica et Cartographica Sinica.1999,28(2):103-109.(劉崇兵,寧津生,張禹慎.利用地震面波和重力資料聯(lián)合反演地殼-上地幔三維密度結(jié)構(gòu)的方法探討[J].測繪學報,1999,28(2):103-109.)

        [15] CUI Xizhang,YU Zongchou,TAO Benzao,et al.Generalized Surveying Adjustment[M].New ed.Wuhan:Wuhan Technical University of Surveying and Mapping Press,2001.(崔希璋,于宗儔,陶本藻,等.廣義測量平差[M].新版.武漢:武漢測繪科技大學出版社,2001.)

        [16] TAO Benzao.Robust Homogeneous Analysis of Displacement-Strain Models[J].Acta Geodaetica et Cartographica Sinica,2000,29(4):293-296.(陶本藻.位移-應(yīng)變模型的穩(wěn)健均勻分析[J].測繪學報,2000,29(4):293-296.)

        [17] CHEN Jian,TAO Benzao.Geodetic Deformation Surveying[M].Beijing:Earthquake Publishing House,1987.(陳健,陶本藻.大地形變測量學[M].北京:地震出版社,1987.)

        A Two-step Method to Determine Relative Weight Ratio Factors in Joint Inversion

        WANG Leyang1,XU Caijun2,ZHANG Chaoyu2
        1.Faculty of Geomatics,East China Institute of Technology,Nanchang 330013,China;2.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China

        The two-step method to determine relative weight ratio factors in joint inversion is derived.In the first step,without considering relative weight ratio factors,the joint inversion model of two or more kinds of observation data is established.Through Helmert method of variance components estimation,unit weight variances of observation data are equal(or the ratio of unit weight variances is approximate to1),and the reasonable and theoretical weight matrix is determined.At the sametime,different kinds of observation data are normalized.In the second step,on the base of reasonable weight matrix,the joint inversion model considering relative weight ratio factors is established.The relative weight ratio factors are determined when objective function is minimum.Through simulated data and several solutions,the effectiveness and feasibility of two-step method are proved.

        joint inversion;relative weight ratio factors;two-step method;Helmert method of variance components estimation;unit weight variance

        WANG Leyang(1983—),male,PhD,lecturer,majors in geodetic inversion and geodetic data processing.

        WANG Leyang,XU Caijun,ZHANG Chaoyu.A Two-step Method to Determine Relative Weight Ratio Factors in Joint Inversion[J].Acta Geodaetica et Cartographica Sinica,2012,41(1):19-24.(王樂洋,許才軍,張朝玉.一種確定聯(lián)合反演中相對權(quán)比的兩步法[J].測繪學報,2012,41(1):19-24.)

        P207

        A

        1001-1595(2012)01-0019-06

        國家自然科學基金(40874003;40974017;41074007;41021061);國家863計劃(2009AA12Z317);國家公益地震行業(yè)科研專項(200808080);教育部博士點基金(20090141110055);東華理工大學地理空間信息采集處理及應(yīng)用科技創(chuàng)新團隊項目

        叢樹平)

        2011-01-10

        2011-07-28

        王樂洋(1983—),男,博士,講師,主要研究方向為大地測量反演及大地測量數(shù)據(jù)處理。

        E-mail:wleyang@163.com

        猜你喜歡
        方差反演觀測
        觀測到恒星死亡瞬間
        軍事文摘(2023年18期)2023-11-03 09:45:42
        方差怎么算
        反演對稱變換在解決平面幾何問題中的應(yīng)用
        概率與統(tǒng)計(2)——離散型隨機變量的期望與方差
        計算方差用哪個公式
        基于低頻軟約束的疊前AVA稀疏層反演
        基于自適應(yīng)遺傳算法的CSAMT一維反演
        方差生活秀
        天測與測地VLBI 測地站周圍地形觀測遮掩的討論
        可觀測宇宙
        太空探索(2016年7期)2016-07-10 12:10:15
        一本色道久久综合亚洲精品蜜臀| 国产av丝袜旗袍无码网站| 亚洲av无码男人的天堂在线| 一区二区精品| 少妇人妻字幕一区二区| 国产一区二区三区久久悠悠色av| 国产成人av一区二区三区| 国产精品久久久久乳精品爆| 亚洲国产理论片在线播放| 亚洲最稳定资源在线观看| 久久国产精品色av免费看| 久久精品国产亚洲av无码偷窥| 亚洲va国产va天堂va久久| 一性一交一口添一摸视频| 中文字幕一区二区三区乱码不卡 | 麻豆夫妻在线视频观看| 不卡日韩av在线播放| 亚洲人成精品久久久久| 久久99精品久久久久久hb无码| 91久久精品国产91久久| 老肥熟女老女人野外免费区| 精品人妻69一区二区三区蜜桃| 天天摸夜夜摸夜夜狠狠摸| 亚洲欧美激情精品一区二区| 亚州AV无码乱码精品国产| 亚洲精品国产福利在线观看 | 人人做人人爽人人爱| 国产欧美日韩a片免费软件| 女人被躁到高潮嗷嗷叫| 在线观看麻豆精品视频| 美女张开腿让男人桶爽| 97精品伊人久久大香线蕉| 日日噜噜噜夜夜爽爽狠狠视频| 羞羞色院99精品全部免| 国产成人精品免费视频大全软件| 国产高潮视频在线观看| 日韩精品国产自在久久现线拍| 国产精品农村妇女一区二区三区 | 91精品国自产拍老熟女露脸| 成人片黄网站a毛片免费| 久久久久久人妻一区精品 |