999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于InSAR技術(shù)的地表三維形變獲取方法綜述

2017-01-26 04:03:06高明亮宮輝力陳蓓蓓周超凡
測繪通報(bào) 2017年1期
關(guān)鍵詞:融合信息方法

高明亮,宮輝力,陳蓓蓓,周超凡,司 遠(yuǎn)

(1. 首都師范大學(xué)資源環(huán)境與旅游學(xué)院,北京 100048; 2. 城市環(huán)境過程 與模擬國家重點(diǎn)實(shí)驗(yàn)室培育基地,北京 100048; 3. 水資源安全北京實(shí)驗(yàn)室,北京 100048)

基于InSAR技術(shù)的地表三維形變獲取方法綜述

高明亮1,2,3,宮輝力1,2,3,陳蓓蓓1,2,3,周超凡1,2,3,司 遠(yuǎn)1,2,3

(1. 首都師范大學(xué)資源環(huán)境與旅游學(xué)院,北京 100048; 2. 城市環(huán)境過程 與模擬國家重點(diǎn)實(shí)驗(yàn)室培育基地,北京 100048; 3. 水資源安全北京實(shí)驗(yàn)室,北京 100048)

對近年來提出的基于InSAR技術(shù)的地表三維形變獲取方法進(jìn)行了綜述,并對該技術(shù)的發(fā)展進(jìn)行了展望。對一些重要的知識背景和理論方法進(jìn)行了簡單介紹,并針對每種方法列舉典型案例,對其適用條件進(jìn)行了簡單總結(jié),為相關(guān)領(lǐng)域的研究者提供了最新研究進(jìn)展和技術(shù)參考。

InSAR;視線向形變;方位向形變;MAI;偏移量追蹤

20世紀(jì)70年代發(fā)展起來的InSAR技術(shù)是一種極具潛力的主動(dòng)式對地觀測技術(shù)。InSAR技術(shù)除了具有分辨率高、覆蓋范圍廣、觀測周期短等星載平臺(tái)普遍的優(yōu)勢外,對形變敏感度高、全天候全天時(shí)觀測等優(yōu)勢使其成為獨(dú)一無二的基于面觀測的形變監(jiān)測手段[1]。如今InSAR技術(shù)已廣泛應(yīng)用于地形測繪[2]、DEM重建[3]、地面沉降[4]、地震形變[5]、礦山形變[6]、火山活動(dòng)[7]、冰川漂移[8]、山體滑坡[9]及大型線性工程形變監(jiān)測[10]等領(lǐng)域。

SAR系統(tǒng)的側(cè)視成像幾何和低觀測角使得沿其視線(line of sight,LOS)方向的觀測對地形的隆起和沉降非常敏感,但是對于垂直LOS向的形變無法感知[11]。實(shí)際情況中垂直向形變通常伴隨著一定程度的水平向形變,尤其是在研究由地質(zhì)構(gòu)造(如地裂縫)控制下的地表形變時(shí),通常是不可忽略的。通過地表三維形變信息的獲取,常常可以揭示常規(guī)InSAR所觀測不到的形變特征。

本文擬總結(jié)利用InSAR技術(shù)獲取地表三維形變的近期進(jìn)展,介紹目前主流的地表形變?nèi)S解算方法和應(yīng)用案例,并對不同方法的適用條件進(jìn)行總結(jié),以便相關(guān)研究者進(jìn)行方法選擇和參考借鑒。

1 傳統(tǒng)InSAR觀測的視線向模糊問題

根據(jù)雷達(dá)側(cè)視觀測幾何可知,InSAR觀測到的是地表在正東、正北和垂直向形變量在雷達(dá)視線方向的投影的和(即矢量和)。因此,由InSAR測量得到的形變結(jié)果,均可以按照成像幾何的空間關(guān)系分解為上述3個(gè)方向的分量。

約定形變遠(yuǎn)離雷達(dá)(沉降)時(shí)LOS向觀測結(jié)果為負(fù),反之為正。分別將LOS向、正東、正北和垂直向形變用dLOS、dE、dN、dU表示,則InSAR觀測的地表形變可以表示為[20]

dLOS=dUcosθ-dEsinθsin (α-3π/2)-

dNsinθcos (α-3π/2)

式中,θ為雷達(dá)觀測方向(LOS向)與垂直向的夾角;α為衛(wèi)星航向與正北方向沿順時(shí)針方向的夾角。

2 利用InSAR技術(shù)獲取地表三維形變

由于星載雷達(dá)觀測系統(tǒng)均運(yùn)行在近極地軌道,且是側(cè)視系統(tǒng),造成方位向信息的缺失,對南北向形變并不敏感。為了獲取可靠的南北向形變信息,國內(nèi)外學(xué)者轉(zhuǎn)而關(guān)注與其密切聯(lián)系的方位向形變信息獲取技術(shù),如偏移量追蹤(offset-tracking)技術(shù)和多孔徑SAR干涉(multi-aperture InSAR,MAI)技術(shù)。

2.1 偏移量追蹤(offset-tracking)技術(shù)

偏移量追蹤法一開始用于彌補(bǔ)InSAR方位向信息缺失。早在20世紀(jì)90年代末,Michel等[12]利用兩景ERS-1圖像配準(zhǔn)后的偏移量,成功獲取了由于地震引起的LOS向(即SAR距離向)和垂直LOS向(即SAR方位向)的二維同震形變信息,結(jié)果滿足地震形變測量的精度需要。此外,距離向和方位向形變可以通過互相關(guān)(cross-correlation)技術(shù)和相干追蹤法(coherent tracking)獲取[13],但是其精度遠(yuǎn)不及干涉處理的結(jié)果,約為SAR圖形像元尺寸的1/30~1/10[14]。互相關(guān)技術(shù)又稱強(qiáng)度追蹤法(intensity tracking),相干追蹤法又稱條紋法(the fringe visibility based algorithm),通過分塊配準(zhǔn)確定全局相干性最高時(shí)給定點(diǎn)的方位向和距離向偏移量,據(jù)此估算方位向和距離向的位移量,具體過程可參考文獻(xiàn)[15]。

2.2 多孔徑SAR干涉技術(shù)

MAI技術(shù)首先由Bechor和Zebker在2006年提出[16],成為繼offset-tracking之后基于InSAR技術(shù)獲取方位向形變信息的另一個(gè)選擇。MAI的核心思想是通過方位向?yàn)V波處理,將構(gòu)成干涉對的兩景SAR復(fù)數(shù)影像分別通過子孔徑處理技術(shù)生成前、后視復(fù)影像;然后分別將前、后視像對進(jìn)行差分干涉處理;最后將前、后視干涉相位共軛相乘,得到多孔徑差分干涉結(jié)果。王曉文等[17]于2014年提出了一種基于MAI和InSAR聯(lián)合獲取地表三維形變信息的方法,并以2003年伊朗Bam地震為案例進(jìn)行了試驗(yàn),結(jié)果表明其精度和效率優(yōu)于前人研究成果。

MAI技術(shù)基于差分干涉相位,其精度取決于干涉相干性和視數(shù),因此MAI技術(shù)得到的方位向形變精度通常高于offset-tracking技術(shù),而且效率比后者要高。然而,由于MAI技術(shù)對相位失相干非常敏感,在獲取短時(shí)間劇烈地表形變方面(如地震、冰川漂移)與offset-tracking相比略遜一籌。

2.3 融合方位向觀測技術(shù)

基于融合方位向形變信息觀測技術(shù)的核心思想是:通過不同途徑得到方位向信息,然后聯(lián)合距離向形變信息進(jìn)行三維信息的解算。從上文可以得知,差分干涉處理或偏移量追蹤的方法均可以獲取LOS向形變信息,而方位向形變信息可以通過offset-tracking或MAI技術(shù)得到。在獲取到至少兩個(gè)觀測軌道(視角)的距離向和方位向形變之后,采用一種加權(quán)最小二乘算法得到三維形變信息。這里需要注意的是,采用升降軌數(shù)據(jù)可以保證其觀測的相對獨(dú)立性,以便得到可信的解算結(jié)果。具體的算法可參考文獻(xiàn)[18]。

2.4 融合GPS觀測技術(shù)

眾所周知,在連續(xù)工作模式下,GPS觀測精度(水平和垂直方向)可以達(dá)到毫米級。其全天時(shí)、全天候、大范圍觀測優(yōu)勢可以與雷達(dá)衛(wèi)星媲美。然而,其過于稀疏的站點(diǎn)分布(gt;10 km)及昂貴的安裝和維護(hù)費(fèi)用使得其普適性大大降低。

事實(shí)上,GPS除了可以提供地面觀測點(diǎn)的形變信息,也可為InSAR提供點(diǎn)位上的氣象觀測數(shù)據(jù),用于消除大氣延遲的影響[19],是InSAR技術(shù)的絕佳數(shù)據(jù)輔助和補(bǔ)充數(shù)據(jù)源。融合InSAR和GPS觀測技術(shù)最早由Gudmundsson等[20]在2002年提出,用于獲取雷克珍半島的三維形變速率場信息。作者將GPS點(diǎn)位信息插值到與InSAR觀測結(jié)果同樣的分辨率,然后基于馬爾可夫隨機(jī)場和模擬退火算法的全局最優(yōu)化方案到三維形變速率場信息。目前,融合InSAR和GPS觀測技術(shù)的發(fā)展方向主要為改進(jìn)GPS插值算法、消除大氣延遲影響[21]及改進(jìn)融合算法[22]。

2.5 多平臺(tái)InSAR技術(shù)

InSAR技術(shù)已成為較成熟的地表形變信息獲取技術(shù),因此基于本文的觀測幾何模型,聯(lián)合多個(gè)視角觀測結(jié)果可以實(shí)現(xiàn)三維形變場的恢復(fù)與重建。通過對至少3個(gè)衛(wèi)星平臺(tái)SAR數(shù)據(jù)的干涉處理,得到LOS向地表形變信息,基于最小二乘的方法解算三維位移分量,典型的處理案例可參考文獻(xiàn)[23]。

2.6 基于先驗(yàn)知識的技術(shù)

在對研究對象形變或位移特征(如地裂縫的發(fā)育、滑坡或冰川的水平移動(dòng))有一定認(rèn)識的基礎(chǔ)上,可以作出相應(yīng)假設(shè),從而簡化地表三維形變信息的求解過程。Kumar等[24]在冰川貼近地表地形起伏運(yùn)動(dòng)的假設(shè)下,通過升降軌干涉獲取的LOS向位移信息,解算出格陵蘭島冰川移動(dòng)的三維速度場。研究表明,基于此假設(shè)解算得到的三維速率場信息與融合方位向觀測技術(shù)(MAI、offset-tracking)獲取的三維解算結(jié)果一致。

正如上文所述,目前在軌的雷達(dá)觀測衛(wèi)星運(yùn)行在近極地軌道,在非極地地區(qū),南北向地表形變信息對雷達(dá)側(cè)視系統(tǒng)獲取的LOS觀測方向形變貢獻(xiàn)不大。因此,在研究對象形變方向主要沿東西向發(fā)展的情況下,可以忽略南北方向的形變信息。在此前提下,只需兩組來自不同視角的LOS向觀測結(jié)果,便可通過簡化后的觀測幾何模型得到垂直向和東西向的形變信息[25]。

3 需要解決的關(guān)鍵問題

在過去幾十年,地表形變?nèi)S獲取技術(shù)有了長足的發(fā)展,然而,并沒有一種技術(shù)可解決任何應(yīng)用問題。在研究不同問題時(shí),需要根據(jù)已有的數(shù)據(jù)條件和精度需要選取相應(yīng)的技術(shù)方法。

常規(guī)InSAR技術(shù)易受到時(shí)、空失相干和大氣延遲等因素的影響,因此基于InSAR結(jié)果的地表三維形變信息獲取技術(shù)勢必也受到其制約。時(shí)序InSAR技術(shù)可以有效解決時(shí)空失相干的問題,大氣延遲也可以通過模擬或結(jié)合外部數(shù)據(jù)的方法借以弱化。如何將時(shí)序InSAR技術(shù)與三維解算方法進(jìn)行有機(jī)融合,得到時(shí)間序列的三維形變信息,是有待解決的難題。

此外,以地面沉降為典型代表的地表長時(shí)間緩慢形變,由于其普遍性和危害性受到社會(huì)各界包括學(xué)術(shù)界的高度重視[26]。然而,對于此類長時(shí)間緩慢地表形變的三維形變速率場獲取問題,融合方位向觀測技術(shù)(offset-tracking和MAI)由于其過低的觀測精度(厘米至分米級別),無法獲取可靠的結(jié)果。目前最有效的方法是結(jié)合GPS觀測結(jié)果作為補(bǔ)償,然而由于后者的空間密度過低,遠(yuǎn)不能滿足監(jiān)測地表形變這種面觀測的要求,對于區(qū)域復(fù)雜地表形變信息的獲取方面依然有所欠缺。

4 InSAR技術(shù)獲取三維形變技術(shù)展望

近年來,一批高分辨率雷達(dá)衛(wèi)星的發(fā)射為InSAR技術(shù)提供了豐富的數(shù)據(jù)源,也為進(jìn)一步發(fā)現(xiàn)區(qū)域形變特征提供了可能。通過獲取高分辨率的地形數(shù)據(jù)(如DEM或DSM),可以進(jìn)一步削弱地形相位誤差,提高影響配準(zhǔn)的精度,從而獲得高精度的形變觀測結(jié)果。此外,新型的高級InSAR觀測技術(shù)如SqueezeSARTM[27]和SAR層析技術(shù)[28](SAR tomography)的發(fā)展,極大拓展了InSAR技術(shù)的研究和應(yīng)用領(lǐng)域,解決了高密度復(fù)雜建筑區(qū)目標(biāo)點(diǎn)的識別,并提高了數(shù)據(jù)利用率。此外,GRACE(gravity recovery and climate experiment)計(jì)劃實(shí)施以來,利用重力衛(wèi)星研究全球質(zhì)量變化受到廣泛關(guān)注[29],InSAR技術(shù)與GRACE技術(shù)的融合互補(bǔ)也成為國內(nèi)外地學(xué)界一個(gè)研究熱點(diǎn)。此外,正如前文所提到的,發(fā)展融合時(shí)序InSAR觀測結(jié)果的三維形變獲取技術(shù),監(jiān)測地表長時(shí)間序列緩慢形變特征,也是InSAR對地觀測技術(shù)的重要需求。

5 結(jié)束語

本文總結(jié)了基于InSAR觀測結(jié)果獲取地表三維形變信息的主要技術(shù)。其中,融合方位向觀測技術(shù)通過聯(lián)合距離向與方位向觀測結(jié)果得到地表的三維形變信息,多用于獲取地表短時(shí)間劇烈形變,如地震同震形變場、冰川漂移速率場、火山運(yùn)動(dòng)引起的形變場等。融合InSAR和GPS觀測結(jié)果的方法已廣泛用于地震監(jiān)測、地裂縫監(jiān)測及滑坡監(jiān)測等應(yīng)用領(lǐng)域。基于先驗(yàn)知識的技術(shù)方法,常用于已掌握其形變特征或規(guī)律的研究對象,如冰川漂移、滑坡和地裂縫發(fā)育等方向性明顯的地表形變問題。由于每種方法均各有利弊,在進(jìn)行方法選取時(shí)要根據(jù)研究問題的具體情況進(jìn)行最優(yōu)選擇。

有理由相信,在不久的將來,一定會(huì)有更多更先進(jìn)的基于InSAR觀測獲取地表三維形變信息的技術(shù)方法涌現(xiàn)出來,將進(jìn)一步推動(dòng)InSAR及其相關(guān)技術(shù)在更多領(lǐng)域的應(yīng)用和發(fā)展。

[1] 廖明生, 田馨, 趙卿. TerraSAR-X/TanDEM-X雷達(dá)遙感計(jì)劃及其應(yīng)用[J]. 測繪信息與工程, 2007, 32(2): 44-46.

[2] ZEBKER H A, GOLDSTEIN R M. Topographic Mapping from Interferometric Synthetic Aperture Radar Observations[J].Journal of Geophysical Research, 1986, 91(B5): 4993-4999.

[3] 蔣厚軍, 廖明生, 張路, 等. 高分辨率雷達(dá)衛(wèi)星COSMO-SkyMed干涉測量生成DEM的實(shí)驗(yàn)研究[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2011, 36(9): 1055-1058.

[4] 宮輝力, 張有全, 李小娟, 等. 基于永久散射體雷達(dá)干涉測量技術(shù)的北京市地面沉降研究[J]. 自然科學(xué)進(jìn)展, 2009, 19(11): 1261-1266.

[5] ALBANO M, BARBA S, SAROLI M, et al. Gravity-driven Postseismic Deformation Following the Mw 6.3 2009 L’Aquila (Italy) Earthquake[J]. Scientific Reports, 2015(5): 16558.

[6] 王志勇, 張繼賢, 黃國滿. 基于InSAR的濟(jì)寧礦區(qū)沉降精細(xì)化監(jiān)測與分析[J]. 中國礦業(yè)大學(xué)學(xué)報(bào), 2014, 43(1):169-174.

[7] WALTER T R, SUBANDRIYO J, KIRBANI S, et al. Volcano-tectonic Control of Merapi’s Lava Dome Splitting: The November 2013 Fracture Observed from High Resolution TerraSAR-X Data[J].Tectonophysics, 2015, 639: 23-33.

[8] NEELMEIJER J, MOTAGH M, WETZEL H. Estimating Spatial and Temporal Variability in Surface Kinematics of the Inylchek Glacier, Central Asia, Using TerraSAR-X Data[J]. Remote Sensing, 2014, 6(10): 9239-9259.

[9] 王志勇, 張金芝. 基于InSAR技術(shù)的滑坡災(zāi)害監(jiān)測[J]. 大地測量與地球動(dòng)力學(xué), 2013, 33(3):87-91.

[10] 師紅云, 劉廣, 楊松林. 基于時(shí)序InSAR技術(shù)的京津高鐵區(qū)域沉降穩(wěn)定性評估[J]. 北京交通大學(xué)學(xué)報(bào)(自然科學(xué)版), 2014, 38(6): 78-81.

[11] WRIGHT T J, PARSONS B E, ZHONG L. Toward Mapping Surface Deformation in Three Dimensions Using InSAR[J]. Geophysical Research Letters, 2004, 31(1): 169-178.

[12] MICHEL R, AVOUAC J P, TABOURY J. Measuring Near Field Coseismic Displacements from SAR Images: Application to the Landers Earthquake[J]. Geophysical Research Letters, 1999, 26(19): 3017-3020.

[13] STROZZI T, LUCKMAN A, MURRAY T, et al. Glacier Motion Estimation Using SAR Offset-tracking Procedures[J].IEEE Transactions on Geoscience amp; Remote Sensing, 2002, 40(11): 2384-2391.

[14] SIMONS M, ROSEN P A. Interferometric Synthetic Aperture Radar Geodesy[J]. Treatise Geophys, 2007, 3: 391-446.

[15] DERAUW D. DInSAR and Coherence Tracking Applied to Glaciology: the Example of Shirase Glacier[C]∥Fringe ’99: Advancing ERS, SAR Interferometry from Applications towards Operations. Noordwijk:ESA Publications Division, 1999.

[16] BECHOR N B D, ZEBKER H A. Measuring Two-dimensional Movements Using a Single InSAR Pair[J]. Geophysical Research Letters, 2006, 331(16):275-303.

[17] 王曉文, 劉國祥, 張瑞,等. 聯(lián)合多孔徑雷達(dá)干涉與常規(guī)合成孔徑雷達(dá)干涉提取三維形變場[J]. 大地測量與地球動(dòng)力學(xué), 2014, 34(4): 127-134.

[18] FIALKO Y, SIMONS M, AGNEW D. The Complete (3-D) Surface Displacement Field in the Epicentral Area of the 1999 M(w) 7.1 Hector Mine Earthquake, California, from Space Geodetic Observations[J]. Geophysical Research Letters, 2011, 28 (16): 3063-3066.

[19] CATALAO J, NICO G, HANSSEN R, et al. Merging GPS and Atmospherically Corrected InSAR Data to Map 3-D Terrain Displacement Velocity[J].IEEE Transactions on Geoscience amp; Remote Sensing, 2011, 49(6): 2354-2360.

[20] 楊成生, 張勤, 趙超英,等. 基于地形和GPS觀測的SAR干涉圖大氣延遲估計(jì)[J]. 大地測量與地球動(dòng)力學(xué), 2011, 31(1): 142-146.

[21] GUDMUNDSSON S, TUMI GUDMUNDSSON M, BJ?RNSSON H, et al. Three-dimensional Glacier Surface Motion Maps at the Gjálp Eruption Site, Iceland, Inferred from Combining InSAR and other Ice-displacement Data[J]. Annals of Glaciology, 2001, 34(1): 315-322.

[22] 胡俊, 李志偉, 朱建軍,等. 基于BFGS法融合InSAR和GPS技術(shù)監(jiān)測地表三維形變[J]. 地球物理學(xué)報(bào), 2013, 56(1): 117-126.

[23] 劉國祥, 張瑞, 李陶,等. 基于多衛(wèi)星平臺(tái)永久散射體雷達(dá)干涉提取三維地表形變速度場[J]. 地球物理學(xué)報(bào), 2012, 55(8): 2598-2610.

[24] KUMAR V, VENKATARAMANA G, H?GDA K A. Glacier Surface Velocity Estimation Using SAR Interferometry Technique Applying Ascending and Descending Passes in Himalayas[J]. International Journal of Applied Earth Observation amp; Geoinformation, 2011, 13(4): 545-551.

[25] GERNHARDT S, CONG X, EINEDER M, et al. Geometrical Fusion of Multitrack PS Point Clouds.[J].IEEE Transactions on Geoscience amp; Remote Sensing Letters, 2012, 50(1): 38-42.

[26] 陳蓓蓓, 宮輝力, 李小娟,等. 基于InSAR技術(shù)北京地區(qū)地面沉降監(jiān)測與風(fēng)險(xiǎn)分析[J]. 地理與地理信息科學(xué), 2011, 27(2): 16-20.

[27] FERRETTI A, FUMAGALLI A, NOVALI F, et al. A New Algorithm for Processing Interferometric Data-Stacks: SqueeSAR[J]. IEEE Transactions on Geoscience amp; Remote Sensing, 2011, 49(9): 3460-3470.

[28] FORNARO G, REALE D, SERAFINO F. Four-Dimensional SAR Imaging for Height Estimation and Monitoring of Single and Double Scatterers[J].IEEE Transactions on Geoscience amp; Remote Sensing, 2009, 47(1): 224-237.

[29] 羅志才, 李瓊, 鐘波. 利用GRACE時(shí)變重力場反演黑河流域水儲(chǔ)量變化[J]. 測繪學(xué)報(bào), 2012, 41(5):676-681.

ReviewofThree-dimensionalSurfaceDeformationAcquisitionfromInSARMeasurementsandItsApplication

GAO Mingliang1,2,3,GONG Huili1,2,3,CHEN Beibei1,2,3,ZHOU Chaofan1,2,3,SI Yuan1,2,3

(1. College of Resource Environment and Tourism, Capital Normal University, Beijing 100048,China; 2. Base of the State Key Laboratory of Urban Environmental Process and Digital Modeling, Capital Normal University, Beijing 100048,China; 3. Beijing Laboratory of Water Resources Security, Capital Normal University, Beijing 100048, China)

In this paper, methods for obtaining three-dimensional deformation based on InSAR techniques in recent years are reviewed, and the development of these techniques is discussed. Important knowledge and theory are introduced simply while typical cases for each method are enumerated. Then applicable condition is summerized to provide relevant researchers of the latest research progress and technical reference.

InSAR; deformation in LOS; deformation in azimuth; MAI; offset-tracking

P23

A

0494-0911(2017)01-0001-04

高明亮,宮輝力,陳蓓蓓,等.基于InSAR技術(shù)的地表三維形變獲取方法綜述[J].測繪通報(bào),2017(1):1-4.

10.13474/j.cnki.11-2246.2017.0001

2016-03-25;

2016-06-24

國家自然科學(xué)基金重點(diǎn)項(xiàng)目(41130744/D0107);國家自然科學(xué)基金(41401492/D010702)

高明亮(1989—),男,博士生,主要從事地面沉降方面的研究。E-mail:b-19890320@163.com

宮輝力。E-mail:gonghl1956@126.com

猜你喜歡
融合信息方法
村企黨建聯(lián)建融合共贏
融合菜
從創(chuàng)新出發(fā),與高考數(shù)列相遇、融合
《融合》
訂閱信息
中華手工(2017年2期)2017-06-06 23:00:31
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
展會(huì)信息
健康信息
祝您健康(1987年3期)1987-12-30 09:52:32
主站蜘蛛池模板: 亚洲成人在线网| 国产精品毛片一区| 国产亚洲欧美另类一区二区| 91网址在线播放| 日韩欧美亚洲国产成人综合| 91福利在线看| 人妻一区二区三区无码精品一区| 天堂av综合网| 久久精品嫩草研究院| 国产又色又刺激高潮免费看| 精品乱码久久久久久久| 最新日本中文字幕| 无码福利日韩神码福利片| 色噜噜在线观看| 亚洲美女久久| 国产久操视频| 午夜天堂视频| 婷婷成人综合| 国产鲁鲁视频在线观看| 亚洲国产午夜精华无码福利| 亚洲成a人片| 亚洲欧美日韩动漫| 一级毛片视频免费| 国产18在线播放| 欧美第二区| 国产精品爆乳99久久| 伊人激情久久综合中文字幕| 人妻夜夜爽天天爽| 无码国产伊人| 九色最新网址| 538国产视频| 一区二区三区四区精品视频| 99久久国产综合精品2020| 玩两个丰满老熟女久久网| 黄色网在线| 久久亚洲国产视频| 日本高清免费不卡视频| 国产视频只有无码精品| 日韩精品一区二区三区中文无码| 精品久久久久成人码免费动漫| 国产成人久久777777| 无码电影在线观看| 亚洲国产高清精品线久久| 国产特级毛片aaaaaaa高清| 国产真实二区一区在线亚洲| 日本成人精品视频| 在线综合亚洲欧美网站| 伊人无码视屏| 狠狠色狠狠综合久久| 毛片一区二区在线看| 国产在线自在拍91精品黑人| 亚洲第一中文字幕| 亚洲三级成人| 激情国产精品一区| 91免费国产在线观看尤物| 久久久久久国产精品mv| 久久成人国产精品免费软件| 无码专区国产精品一区| 99久久精品视香蕉蕉| 美女国内精品自产拍在线播放| 国产91久久久久久| 亚洲黄网在线| 九月婷婷亚洲综合在线| 亚洲国产成人精品无码区性色| 一级不卡毛片| 国产69精品久久| 呦女精品网站| 亚洲欧美精品一中文字幕| a级毛片免费网站| 欧美一区福利| 亚洲午夜福利在线| 在线播放精品一区二区啪视频| 一本一道波多野结衣一区二区 | 中文字幕在线免费看| 综合色区亚洲熟妇在线| 色欲国产一区二区日韩欧美| 欧美国产精品不卡在线观看| 青青热久麻豆精品视频在线观看| 国产嫖妓91东北老熟女久久一| 亚洲国产成人自拍| 亚洲国产精品VA在线看黑人| 色悠久久综合|