?

基于全站儀和快鳥影像協同的山地樣地定位與匹配研究

2017-12-11 01:49王書涵張曉麗劉會玲
中南林業科技大學學報 2017年10期
關鍵詞:單木全站儀樣地

王書涵,張曉麗,楊 銘,劉會玲

(北京林業大學 林學院,精準林業北京市重點實驗室、省部共建森林培育與保護重點實驗室,北京 100083)

基于全站儀和快鳥影像協同的山地樣地定位與匹配研究

王書涵,張曉麗,楊 銘,劉會玲

(北京林業大學 林學院,精準林業北京市重點實驗室、省部共建森林培育與保護重點實驗室,北京 100083)

山區高低起伏的地形容易影響樣地三維坐標的采集精度,從而限制了高分遙感在山區林業上的廣泛應用。以鷲峰林場為研究區,提出了采用全站儀和手持 GPS 結合遙感圖像協同對樣地進行定位的方法,并提出了坐標修正方案:即首先運用手持 GPS 結合全站儀采集樣地及單木的絕對大地坐標;然后利用大比例尺地形圖采用基于控制點、 RPC 的兩種正射校正方法對覆蓋研究區的 QuickBird 影像進行正射校正, RMSE 多光譜控制在1個像元之內,全色圖像控制在 10 個像元之內;并以樣地為例分析了樣地及單木的匹配精度。研究結果表明:本方案能快速準確地建立樣地中單木與高分影像的聯系,其中基于控制點的圖像正射校正法的匹配精度最高,相對提高幅度為 92.6%,而基于 RPC 的提高幅度70.05%,即 EControlPoint> ERPC。研究結果表明,對應用于地形復雜的山區的高空間分辨率的遙感影像正射校正是不可缺少的步驟,而本文提出的方案能實際應用于山區的高分遙感樣地定位,具有推廣利用的價值,并為山區的高分遙感進一步的應用打下基礎。

快鳥影像;高空間分辨率遙感圖像;全站儀;手持GPS;正射校正;

中低等空間分辨率影像,如Landsat TM、MODIS等已經被廣泛運用到大尺度的森林資源、生態環境變化監測、土地利用分類等研究中[1-3]。高分辨率遙感影像由于空間分辨率高、地物表現能力良好的優勢被許多學者應用于小范圍的森林資源調查[4],并取得一定成果[5-7]。然而林業遙感的研究區域經常是在地形起伏多變的山區,如何精確而準確的建立實際地物特征和影像信息間的聯系一直是林業研究者討論和關心的問題[8-9]。手持GPS采集目標點的精度誤差是主要的問題,容易受到山區地形和樹冠的遮擋無法及時收到差分信號,導致定位精度不高,往往單點解只能達到10 m左右,造成樣地坐標和遙感圖像無法很好的匹配。這種定位的不確定性誤差在中低等分辨率影像遙感反演植被參數的過程中往往可以通過樣地的設置將其平差或者忽略,不會對結果造成很大的影響(TM空間分辨率為30 m,MODIS空間分辨率最高為250 m)。但是這樣的精度誤差無法滿足高分遙感影像在林場級的林業應用,10 m級的誤差會導致偏差嚴重,無法精確定位到相應的樹種,因此也限制了高分遙感影像在林業上的應用。

在海拔變化明顯的山區,針對少有相關文獻討論山區樣地的定位精度問題,本文以快鳥影像為例,以鷲峰林場為研究區,運用全站儀與快鳥影像的協同,及后期的處理,精確定位樣地四個角及單木坐標,并在影像上準確還原出樣地的地理位置。并探討了正射校正對對位精度的影響,即采用影像數據自帶的RPC參數的正射校正方法和基于控制點的校正方法;采用全站儀和手持GPS相結合的外業調查方法,對樣地及單木進行快速而準確的定位。本文提出的方法可為高分遙感在林業上的應用提供借鑒。采集得到精確的GPS的坐標點還可作為高分影像進行分類和精度評價的訓練樣本和檢驗樣本,為遙感的定量分析、森林資源監測等提供依據。

1 試驗區與儀器設備

1.1 試驗區概況

研究區在位于北京市海淀區的鷲峰國家森林公園(圖 1 所示)。位于 39°54′ N、116°28′ E,最高海拔1 153 m,最低海拔60 m,總面積約為832.04 hm2,整體坡度起伏較大。山上林分較為破碎,地下灌木較厚,山下植被較為齊整,樹種結構單一[10]。針葉林面積為474.85 hm2,占整個林區面積的57.1%。主要是種植于20世紀50年代的人工林,主要的林分類型有:側柏林、油松林、落葉松林、刺槐林、栓皮櫟林、栓皮櫟與槲櫟混交林等[11],屬于華北大陸季風氣候,具有冬春季干燥多風、夏季涼爽多雨的特征[12]。

1.2 儀器設備

地面樣地調查采用Garmin手持式GPS接收機,因其小巧攜帶方便常被科研人員采用。地面單木的位置測量采用賓得免棱鏡全站儀R-422N,其免棱鏡測量模式下超過550 m的測量距離保證了樣地內單木坐標采集的準確和迅速,必要情況也可以轉換成相應的棱鏡測量模式進行地面實地調查。

1.3 信息源

采用的遙感圖像為覆蓋研究區的快鳥圖像LV2A數據,16Bit位深度。根據Digital Global公司的分級標準,將快鳥圖像分為基本圖像,標準圖像和正射校正圖像,其中標準影像數據本身已經做了幾何粗校正,但未做正射校正,RMSE通常為3 ~10 m。本研究獲取的QBLV2A影像屬于標準影像,采用了WGS84坐標系統,可以直接用來進行正射校正[13],成像時間為2009-10-16T 07:17:48。多光譜空間分辨率為2.4 m,全色分辨率為0.6 m。太陽高度角為37.3°,太陽天頂角是166.6°,平均衛星天頂角為74.6°,平均衛星高度角為 67.3°。

2 試驗方法

2.1 遙感圖像的正射校正

2.1.1 地形圖的選取與校正

正射校正需要研究區域的高程信息,全色高分影像的分辨率達到0.6 m,一般需要使用具有較高比例尺精度的地形圖進行幾何精校正。例如1∶5 000比例尺的比例尺精度為0.5 m,與QuickBird全色影像的空間分辨率較為一致。本研究采用由北京市測繪設計研究院提供的1∶2 000比例尺的地形圖分別對全色和多光譜影像進行正射校正。首先將1∶2 000比例尺的地形圖等高線數字化,將等高線轉換為TIN,再將TIN通過ArcGIS的地形分析模塊將TIN轉換為DEM。地形圖本身為柵格化的圖像已經進行過幾何精校正,采用Beijing_1954坐標系統,中央經線117° E,比例因子1。圖2為線狀等高線提取TIN的結果,并將TIN轉換成了1 m精度的DEM提供正射校正所需的高程信息圖。

圖1 研究區位置圖Fig.1 Geography position of the study area

2.1.2 控制點的選取

遙感圖像的正射校正采用兩種方式,分別為基于圖像自帶的RPC文件進行正射校正和基于控制點的校正。從已校正的1∶2 000比例尺地形圖中選取明顯地面特征點作為控制點,如房屋的角點,道路的交叉點,水塘,寺廟等。盡量保證圖像的四周和中央都有點分布以控制整體精度。分別在ENVI5.0的正射校正模塊中進行。

2.2 圖像融合

圖像融合技術既能保留影像的空間分辨率又能不丟掉光譜信息,有利于解譯和定量分析[14-15],將經過兩種方法校正得到的多光譜和全色圖像分別采用PCA(Principle Component Analysis)、HPF(高通濾波)和Wavelet(小波)進行融合。通過比較圖像均值和標準差評價融合圖像的效果,選擇了均值與原多光譜最為接近的,而標準差最大的一種方法進行圖像的融合。

2.3 地面單木絕對大地坐標的采集

手持GPS信號因易受到樹冠地形的干擾,直接量測濃密林區的單木坐標精度極差,且目標相互之間的相對位置精度無法保證。而全站儀和手持GPS的配合使用保證了絕對大地坐標的準確獲取。首先建立全站儀測量系統:選擇直角坐標測量,(1)建立基站點:于樣地外沒有樹木遮擋且較為平坦處設立基站點,然后利用手持GPS 采集基站點的坐標,利用采集得到的大地平面坐標作為基站點坐標,設(XB-A,YB-A)。(2)建立后視點:從基站點往南或者往北用皮尺量測直線D(用羅盤儀定出南北方向),將此點作為后視點。后視點的坐標為(XB-A,YB-A±D)其中后視點縱坐標加減的確定取決于基站點相對于后視點的位置。全站儀測量系統建立之后接著便可測量單木坐標。將全站儀測量得到的點導出即可得到單木的絕對位置。

圖2 研究區等高線DEM圖(單位:m)Fig.2 DEM in study area converted from contour line

此時采集的單木絕對坐標仍然有偏差,無法與快鳥圖像真實的樣地位置匹配。因此,在前述測量步驟的基礎上增加一個特征點輔助樣地定位,基本思路如下:首先,樣地周邊需要確定一個特征點,其需具備明顯的形狀特征保證其能從快鳥影像中快速目視定位,比如道路的交叉處、鐵路的交叉點、房屋、水池等均可以作為特征點。特征點的選取不宜離樣地太遠,若特征點超過了全站儀的觀測范圍,可以采用搬站的方式引點測量。

在外業過程中,由于調查人員的知識層次不一,不可避免會出現一些差錯。比如,將手持GPS 采集的坐標點(L,B)輸入全站儀中時,橫、縱坐標輸入顛倒,導致測量結果出現偏差,如圖3所示。

可以將(X0,Y0)做關于直線y=x對稱的點得到(x′,y′)(結果見表 1),接著,作(x′,y′)關于直線X=XB-A的對稱點,按照公式1計算。即可得到正確的坐標點(X,Y)(結果見表2,圖4)。其中,(X0,Y0)為直接測量得到的單木坐標,(x′,y′)為對稱點。

圖3 未修正的GPS 采集點的偏差示意圖Fig.3 Deviations schematic diagram with the uncorrected GPS points

3 結果與分析

3.1 全站儀測量點的修正結果

示例樣地中全站儀測量的初步結果記做(X0,Y0),由于橫、縱坐標顛倒導致全站儀測量值出現偏差。表1(X′,Y′)為處理的中間過程,(X,Y)為最終正確的坐標點(見表2)。B-A為基站點的坐標。

表1 將(X0, Y0)作關于直線y=x的對稱點Table 1 The symmetry points of (x0,y0) about line X=XB-A

3.2 QB圖像正射校正結果

控制點的選擇原則為盡量分散,表3顯示多光譜圖像正射校正控制點的均方根誤差(RMSE),表4顯示了全色圖像的正射校正結果。多光譜圖像的總校正精度控制在1個像元之內,全色圖像由于分辨率較高,而地形圖表達受限,不能很精細地選擇控制點,因此控制在10個像素之內即可滿足要求。

表2 將(X′, Y′)作關于直線X=XB-A的對稱點Table 2 The symmetry points of (x′,y′) about line X=XB-A

圖4 修正的GPS采集點示意圖Fig.4 Schematic diagram with the corrected GPS points

3.3 融合結果的對比

不同融合方法的融合效果有差異,本文通過目視和統計方法衡量融合的效果。從目視效果來看,HPF融合方法效果最好,PCA對比度稍低效果其次,而小波融合出現了較為模糊現像;比較3種方法的均值和標準差,從圖5可以看出HPF方法各個波段的均值與原始多光譜的均值最為接近,而HPF的標準差相比小波和PCA方法偏高,圖像層次更豐富(圖6)。目視效果和定量分析的結果一致。因此最終采用HPF方法進行融合。

表3 多光譜圖像控制點誤差Table 3 Multispectral image orthorectification control points error

表4 全色圖像正射校正控制點誤差Table 4 Panchromatic image orthorectification control points error

圖5 三種融合方法的均值比較Fig.5 Mean value with three kinds of fusion methods

3.4 定位結果的評定

本文以樣地B1、B2兩個樣地為例評定兩種正射校正的定位結果,所選取的這兩個樣地毗鄰,分別位于道路的兩側,并且樹冠高大容易辨認,因此在QuickBird影像中很容易識別,圖7右顯示了樣地的位置。分別采用目視效果和距離統計方法評價定位的精度。定位效果如圖7所示,從目視效果來看,經過控制點正射校正的定位精度明顯優于RPC系統自帶的校正精度,而未經過正射校正的圖像的定位偏移距離最大(如圖7左所示),林分單木的定位結果幾乎不可靠。定量分析定位的結果,采用直線距離衡量地面實測單木GPS點到圖像實際單木位置的匹配精度,結果見表5。其中圖像實際單木位置從圖像目視解譯出。從表5可以看出,無論是B1還是B2,采用控制點校正后的圖像定位精度最高,相對提高幅度最大,誤差距離也最小。不作正射校正的高分影像單木定位的結果最差,偏離距離最大。偏移距離采用公式2計算。

圖6 三種融合方法的標準差比較Fig.6 Standard deviation with three kinds of fusion methods

其中(x0,y0)為GPS采集點,(x1,y1)為從圖像中提取的單木位置點。

采用本文提出的樣地匹配方法,對研究區內另兩個樣地進行樣地GPS單木位置點同快鳥影像的匹配分析,對定位的方案進行檢驗,得到樣地精確的匹配效果圖。樣地定位結果如圖8所示,從目視檢驗來看匹配效果較好,一定程度上說明了本文提出的樣地匹配方法的可行性。

圖7 正射校正后定位精度評定(局部放大)Fig.7 Position accuracy assessment after orthorectification (Local Zoom)

表5 GPS測量單木點與實際圖像上單木點的誤差距離?Table 5 The error between individual trees points with GPS and individual trees points in image

4 結論與討論

本文提出了應用于山區的高分影像的單木坐標定位方法:利用全站儀和手持GPS輔以測定特征點的方法較為精確地測定了單木的坐標并還原了影像中對應的單木及樣地位置。該方法能直接獲取樣地及單木絕對大地坐標,同時還保留了目標點的相對位置的準確性。一定程度上解決了手持GPS無法準確獲取樣地內單木間相對位置的弊端,又解決了山地林區樣地定位不精準的矛盾。針對測量過程中容易出現的錯誤,本文還提出了相應的坐標修正方法。定位評價的結果顯示,應用于山區的高空間分辨率影像單木準確定位的實現,正射校正是不可缺少的步驟,基于控制點的正射校正方法定位精度最高,平均提高幅度達到了92.6%,而RPC正射校正的平均提高幅度為70.05%。研究結果表明:手持式GPS雖然定位精度有限,但通過合適的坐標投影設置,圖像經正射校正后,能達到較好的定位結果。從而在一定程度上解決了地面同名地物點與遙感影像匹配難的問題,拓展了高分影像在林業上的應用。

圖8 其他樣地匹配結果Fig.8 Matching results in other plots

隨著遙感圖像空間分辨率的不斷提高,單木參數因子如樹冠輪廓、樹高、胸徑等因子的準確提取成為了可能;而GPS定位技術的不斷改進又更進一步推動了高分辨率遙感影像在林業上的廣泛應用??梢哉f遙感影像的空間分辨率和GPS定位精度的不斷提高,為森林結構參數的提取提供了有利的條件[16-17]。但另一方面,隨著高分辨率遙感影像種類及獲取手段層出不窮[18],很多技術上的問題亟待得到解決?;诟叻直媛蔬b感的單木定位與匹配就是值得探討的問題之一。影響高分辨率遙感影像的單木定位和匹配困難的原因主要有兩個因素,一是遙感影像本身的幾何畸變,其主要來源于傳感器與衛星平臺的安裝,衛星位置、姿態觀測誤差和幾何處理模型造成的誤差等[19],地形起伏也是造成影像畸變的一個重要因素[20];二是地面精度的驗證問題,存在的主要困難就是很難獲取高精度的目標物的絕對大地坐標,即坐標真值。這方面的誤差通常來源于地面樣地或單木調查應用的手持GPS定位精度受到樹冠遮擋、地勢、天氣等因素影響造成結果的不確定性。這兩方面的因素影響了高分遙感的單木定位匹配精度。

本文所提出的方案也存在一些弊端,比如由于GPS采集的坐標不可能百分之百的準確,對于定位的結果仍然需要人工調整并移動位置點,以保證定位的準確性;當樣地的郁閉度較高并位于復雜的林分內,四周又沒有合適的特征點可供參考和選取的時候,在高分遙感影像上準確定位并找到單木的位置仍然存在很大的困難??傊?,本文結合遙感、測繪、地理信息系統等技術手段提高了樣地及單木的定位精度,特別是準確獲取了山區樣地的大地絕對坐標,對科研、生產均有指導意義,例如利用提取的單木位置信息可以進一步分析樹種的空間格局[21]、提取單木競爭指標[22]等。

致謝:本論文受到國家“863”計劃課題“數字化森林資源監測技術”(2012AA102001)、教育部北京市森林培育與保護省部共建重點實驗室項目“北京生態公益林節水增效撫育技術及碳儲量、三維綠量遙感測定技術研究”(2009GJSYS02)資金資助,野外基礎數據來源于實驗室學生野外采集,在此一并表示感謝。

[1]Muukkonen P, Heiskanen J. Biomass estimation over a large area based on standwise forest inventory data and ASTER and MODIS satellite data: A possibility to verify carbon inventories[J]. Remote Sensing of Environment, 2007, 107(4): 617-24.

[2]Foody G M. Status of land cover classification accuracy assessment[J]. Remote Sensing of Environment, 2002, 80(1):185-201.

[3]徐涵秋,唐 菲. 新一代Landsat系列衛星:Landsat 8遙感影像新增特征及其生態環境意義[J]. 生態學報, 2013, 33(11):3249-3257.

[4]吳 見,彭道黎. 基于面向對象的QuickBird影像退耕地樹冠信息提取[J]. 光譜學與光譜分析, 2010, 30(9): 2533-2536.

[5]林 輝,寧曉波,呂 勇. 基于高分辨率衛星圖像的立木材積表的編制[J]. 林業科學, 2004, 40(4): 33-39.

[6]熊軼群,吳健平. 基于高分辨率遙感圖像的樹冠面積提取方法[J]. 地理與地理信息科學. 2007, 23(6): 30-33.

[7]張慧春,鄭加強,周宏平,等. 農藥精確施用系統信息流集成關鍵技術研究[J]. 農業工程學報. 2007, 23(5): 130-136.

[8]羅 旭,李春干,李崇貴,等. 星基差分GPS用于林區高分辨率遙感圖像幾何精校正[J]. 北京林業大學學報,2005,27(S2):48-51.

[9]陳爾學,李增元,車學儉,等. DGPS在資源環境遙感中的應用方法研究[J]. 遙感學報, 2001, 5(4): 282-288.

[10] 劉 維,張曉麗,馬 菁. 鷲峰國家森林公園主要喬木樹種含碳率分析[J]. 西北林學院學報, 2011, 26(5): 214-218.

[11] 趙永泉,彭道黎. 北京鷲峰公園主要人工林群落多樣性研究[J]. 西南林學院學報, 2008, 28(1): 17-22.

[12] 王 昆,張曉麗,王 珊,等. 鷲峰地區QuickBird影像紋理特征與生物量估測關系初探[J]. 地理與地理信息科學, 2013,29(3): 52-55.

[13] Digital Globe, Inc. QuickBird Imagery Products-product guide[EB/OL]. http://glcf.umd.edu/library/guide/QuickBird_Product_Guide.pdf.

[14] Lu D, Li G, Moran E,et al. A comparison of multisensor integration methods for land cover classification in the Brazilian Amazon[J]. GIScience & Remote Sensing. 2011, 48(3): 345-370.

[15] Lu D S, Hetrick S, Moran E,et al.Detection of urban expansion in an urban-rural landscape with multitemporal QuickBird images[J]. Journal of applied remote sensing. 2010, 4(1).

[16] Ozdemir I, Karnieli A. Predicting forest structural parameters using the image texture derived from WorldView-2 multispectral imagery in a dryland forest, Israel[J]. International Journal of Applied Earth Observation and Geoinformation,2011,13(5):701-710.

[17] Meng J, Li S, Wang W,et al. Estimation of Forest Structural Diversity Using the Spectral and Textural Information Derived from SPOT-5 Satellite Images[J]. Remote Sensing, 2016, 8(2): 125.

[18] 曹明蘭, 張力小, 王 強. 無人機遙感影像中行道樹信息快速提取[J]. 中南林業科技大學學報, 2016, 36(10): 89-93.

[19] 袁修孝, 曹金山. 高分辨率衛星遙感精確對地目標定位理論與方法[M]. 北京: 科學出版社, 2012.

[20] 孫 軍,黎 琪,李和睿. QuickBird遙感圖像校正模型[J]. 四川兵工學報, 2012, 33(4): 107-109.

[21] 李艷麗,楊 華,亢新剛,等. 長白山云冷杉針闊混交林天然更新空間分布格局及其異質性[J]. 應用生態學報, 2013,25(2): 311-317.

[22] 孜來比·買木提名,楊 華,趙廣亮,等. 單木競爭指標的研究進展[J]. 西北林學院學報, 2012, 27(6): 152-158.

[本文編校:吳 彬]

Field location and matching method research in mountain areas based on the synergy of QuickBird image and total station

WANG Shuhan; ZHANG Xiaoli, YANG Ming, LIU Huiling
(Key Laboratory or Silviculture and Conservation, Ministry of Education, Beijing Forestry University, Beijing 100083, China)

Extensively applying of high resolution image was restricted because dramatical affects were caused by the rugged topography of mountain area to the geometric rectification of high resolution remote sensing imagery and coordinate collection with high precision.The forestry research region in remote sensing is often rugged topography of mountains. How to link precisely and accurately the actual land features and image information has long been the concerned and discussed issues by researchers. Jiufeng national forestry park was selected as the study area. One Chinese pine (Pinus tabulaeformis) field sample plot and Chinese arborvitae (Platycladus oriental) field sample plot B2 and B3 were selected respectively. The synergy method combining with remote sensing and total station and Handheld GPS was proposed as well as the revised scenario of coordinate. First the true value of ground objects was collected using total-station and handheld GPS, including the trees in plots and angular points. The original QuickBird multispectral and panchromatic image were respectively orthorectified using two types of orthorectification methods: method based on control points identified from topographic map and RPC method. The RMSE of multispectral image and panchromatic image is within 1 pixel and 10 pixels respectively. Road intersections, houses and other obvious features in the 1:2000 topographic map were selected as the features where the GCPs are. The topographic map surveyed by Beijing institute of surveying and mapping in Oct 2004 has been geometrically corrected. Both methods of process above were conducted in ENVI 5.0. First, feature points were selected near the sample plots whose coordinates were obtained using differential GPS system. Matching accuracy between 3D coordinates of feature points and feature positions in the QuickBird image was analyzed. The results shows that matching accuracy of orthorectification based on control points reached the highest level with 92.6% increased range and RPC orthorectification is better than non-corrected image with 70.05% increased ranged: Econtrolpoint>Erpc.The research indicates that the scenario proposed this paper can realize the position of plots in high resolution remote sensing imagery in the mountain area. And orthorectification is required in the application of high resolution remote sensing imagery especially in mountain area. The proposed method can provide reference for the application of remote sensing in the forestry. The GPS coordinates of the plots and trees can be used as training and testing samples on high image classification and accuracy evaluation, or quantitative analysis of remote sensing, and forest resources monitoring.

QuickBird; high spatial resolution remote sensing imagery; total station; handheld GPS; orthorectification

S715.3 文獻標志碼:A 文章編號:1673-923X(2017)10-0022-08

10.14067/j.cnki.1673-923x.2017.10.005

http: //qks.csuft.edu.cn

2016-05-29

國家“863”計劃課題“數字化森林資源監測技術”(2012AA102001);教育部北京市森林培育與保護省部共建重點實驗室項目“北京生態公益林節水增效撫育技術及碳儲量、三維綠量遙感測定技術研究”(2009GJSYS02)

王書涵,博士研究生 通訊作者:張曉麗,博士后,教授;E-mail:Zhang-xl@263.net

王書涵,張曉麗,楊 銘,等. 基于全站儀和快鳥影像協同的山地樣地定位與匹配研究[J].中南林業科技大學學報,2017,37(10): 22-29.

猜你喜歡
單木全站儀樣地
地基與無人機激光雷達結合提取單木參數
仁懷市二茬紅纓子高粱的生物量及載畜量調查
工程測量中智能化全站儀的應用解析
基于無人機激光雷達點云數據的單木分割研究
額爾古納市興安落葉松中齡林植被碳儲量研究
基于角尺度模型的林業樣地空間結構分析
15 年生鵝掌楸林分生長差異性研究
全站儀中間法在礦山高程測量中的應用——以河南鎮平縣地形測量為例
全站儀極坐標法監測點穩定性分析方法研究
基于雙尺度體元覆蓋密度的TLS點云數據單木識別算法
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合