?

基于無人機遙感的多特征組礦區草本植物地上生物量反演

2024-03-14 08:17張世文芮婷婷李唯佳蔡慧珍
草業科學 2024年1期
關鍵詞:植被指數波段反演

俞 靜 2, ,張世文,芮婷婷,李唯佳 2, ,蔡慧珍

(1.安徽理工大學空間信息與測繪工程學院, 安徽 淮南 232001;2.礦山采動災害空天地協同監測與預警安徽普通高校重點實驗室, 安徽 淮南 232001;3.礦山環境與災害協同監測煤炭行業工程研究中心,安徽 淮南 232001;4.安徽理工大學地球與環境學院, 安徽 淮南 232001)

植被生長狀況是對廢棄礦區生態修復情況的直接反映,礦區地表植被生物量是反映植被生長狀況的直接且重要的指標,可為礦區植被恢復的效果評價提供基礎數據,因此實現礦區植被生物量快速準確的預測對礦區修復效果評價具有重要意義[1-4]。傳統植被生物量估測采用野外人工實地采集方法,不僅費時費力,而且會對當地植被生長產生不同程度的不良影響[5]。隨著遙感技術的快速發展,無人機遙感憑借其高靈活性、高空間分辨率和能夠進行較小面積的植被監測等優勢,已成為植被長勢監測的重要手段,并廣泛應用于資源調查、環境監測等領域[6]。

近些年,眾多專家學者使用各種監測手段對不同植被的各種長勢指標進行了一系列研究,并獲得了諸多成果。如苗春麗等[7]基于無人機高光譜影像光譜數據構建隨機森林模型對天然草地地上生物量進行了反演研究,實現了高寒草甸地上生物量高精度監測;黃家興等[8]基于Sentinel-2 和Landsat 8數據計算5 種植被指數建立草地地上生物量反演模型,結果表明使用10 m 分辨率的Sentinel-2 構建的生物量反演模型精度高于30 m 的Landsat 8 數據;肖武等[9]運用無人機多光譜影像,在傳統植被指數基礎上引入紅邊波段擴展植被指數用于構建高潛水位采煤沉陷區玉米生物量反演模型,結果表明紅邊波段能提高玉米生物量反演精度;鄧尚奇等[10]利用圖像閾值分割去除無人機影像背景像元后提取小麥像元平均反射率,構建植被指數利用偏最小二乘回歸構建反演模型,結果顯示剔除背景像元后并使用植被指數構建的反演模型精度有較明顯的提 升; Stratoulias 等[11]采 集 同 期 的Sentinel-1 和Sentinel-2 衛星影像不同光譜信息,并比較了基于衛星影像不同光譜信息估算生物量的方法,證明了相比于傳統植被指數,所有Sentinel-1 和Sentinel-2 的波段的簡單光譜比值組合與生物量具有更高的相關性;Wan 等[12]基于無人機利用輻射傳輸模型建立水稻全生育期生物量與葉片葉綠素含量、葉面積指數以及冠層葉綠素含量之間的關系,驗證了模型性能。以上研究對于植被生物量的遙感監測,多集中于提取遙感影像的單一光譜特征作為特征變量實現地上生物量的反演,然而當植被過于稀疏或密集時,光譜特征與地上生物量相關性會降低。并且,衛星遙感時空分辨率低且受天氣影響較大,難以有效提取礦區的精細植被長勢信息,而新興且輕小型無人機遙感監測技術可以彌補衛星遙感在小尺度區域監測方面的不足,其可以搭載多光譜傳感器實時準確獲取修復礦區高分辨率遙感影像,從而實現修復礦區植被精細化監測。

因此本研究以無人機遙感為監測手段,以修復礦區草本植物地上生物量為研究對象,在單波段光譜和植被指數兩光譜特征的基礎上引入紋理特征和地形特征,將3 種特征分為光譜特征組和多特征組分別作為模型輸入變量構建反向傳播神經網絡(back propagation neural network, BPNN)、卷積神經網絡(convolutional neural network, CNN)和Elman 神經網絡草本植物生物量反演模型,并將該模型應用于全區生物量反演,對全區草本植物地上生物量進行分級評估,研究結果可為無人機遙感礦區修復地草本植物地上生物量反演研究提供參考。

1 研究區與數據

1.1 研究區概況

本研究以銅官山(117.81° E,30.91° N)為研究區,區域面積約為0.11 km2。其位于安徽省銅陵市市中心東南方向2 km。銅官山全年平均氣溫為18~19 ℃,年降水量為1 384.7 mm,平均海拔495.7 m,呈東南低、西北高的地勢,山頂標高493.1 m,屬于北亞熱帶濕潤季風氣候,特點是季風明顯,四季分明,雨量豐沛,日照充足,雨熱同季,無霜期長。春季雨水偏多,常出現低溫連綿陰雨天氣;夏季天氣炎熱;秋季晴朗少雨,常出現秋季干旱;冬季天氣晴朗、寒冷干燥。實地調查結果顯示研究區主要以灌草叢為地表覆蓋植被類型,其中自然生長的草本植物為主要植物類型,代表性草本植物主要有小蓬草(Conyza canadensis)、狗尾巴草(Setaria viridis)、雞眼草(Kummerowia)、狗牙根(Cynodon dactylon)、狼杷草(Bidens tripartita)和通奶草(Euphorbia hypericifolia)等,研究區域地理位置和采樣點分布如圖1 所示。

圖1 研究區位置及采樣點分布圖Figure 1 Location of the study area and distribution of sampling points

1.2 數據獲取與預處理

1.2.1 無人機數據獲取與預處理

無人機平臺采用大疆Phantom 4 多光譜版植保無人機搭載一體式多光譜成像系統,集成1 個可見光傳感器和5 個多光譜傳感器(藍光、綠光、紅光、紅邊和近紅外),其中藍光波段中心波長為450 nm,帶寬為32 nm;綠光波段中心波長為560 nm,帶寬為32 nm;紅光波段中心波長為650 nm,帶寬為32 nm;紅邊波段中心波長為730 nm,帶寬為32 nm;近紅外波段中心波長為840 nm,帶寬為52 nm,同時無人機內置GPS 和IMU 系統(圖2)。

圖2 無人機和多光譜相機Figure 2 UAV and multispectral camera

進行航拍作業前,首先需要考慮到天氣條件,避免雨雪、大風等天氣,飛行時應選擇天氣晴朗無云,風力較小的時候。同時需要對研究區域的地理情況及外部設施進行調查,確保在無人機飛行過程中不會受到高壓電線、高大樹木等物體的干擾。綜上考慮,無人機多光譜影像的采集時間于2021 年9 月26 日,飛行時間為11:00-14:00。起飛前手動控制無人機飛行至校準白板的正上方約2.5 m 處,采用相機單攝模式來拍攝標準白板。無人機的飛行模式按照提前規劃的S 型航線飛行,航向重疊度70%,旁向重疊度60%,無人機飛行高度50 m,傳感器鏡頭垂直向下,拍照模式為等時間間隔,共獲取了約3 000 張影像。

獲取影像后利用Pix4D-mapper 拼接軟件進行影像拼接得到正射影像,并以正射影像為參考影像,使用ArcGIS 10.6 軟件分別在5 個單波段影像上均勻選取30 個參考點進行幾何校正,且幾何校正誤差在0.5 個像元之內。獲取各個波段中的白板的DN(digital number)值并計算5 個波段的DN 均值,使用ENVI 5.3 軟件中的Band Math 工具對多光譜影像進行輻射校正。計算公式為:

式中:f目標為目標地物的反射率,DN目標為目標地物的DN值,DN白板為白板DN均值,f白板為標準白板反射率。

經上述預處理后即獲取到無人機多光譜影像,空間分辨率為0.03 m。在ENVI 5.3 軟件中導入預處理后的多光譜影像以及采樣點GPS 點位信息,以采樣點為中心,并對應于樣方大小以及影像分辨率,在圖像上裁剪出33 × 33 (像素)的影像作為感興趣區,提取每個樣方內的各個波段的平均反射率。

1.2.2 生物量數據獲取

在進行無人機數據采集前,需進行樣方布設,待無人機數據采集結束后立即進行樣方采集。根據研究區域范圍,利用均勻布點原則布設采樣點,研究區內共布設采樣點43 個,如圖1 所示,設置1 m ×1 m 的樣方,然后通過GPS 記錄每一個樣方點的位置信息,將樣方內的植被齊地面剪下并對每個樣方進行編號。實驗室處理時,將植被放入烘箱內殺青30 min,溫度設置為105 ℃;考慮到樣方內植株水分含量較大,為了避免烘干溫度不適而致使植被沒有充分烘干[13],將其再置于80 ℃的烘箱中24~48 h直到得到恒定的生物量,最后稱取每個樣方的生物量干質量[14]。

2 研究方法

2.1 變量提取

2.1.1 光譜特征

植被指數是不同光譜波段之間的線性或非線性組合,可以反映綠色植被的相對豐富度。植被的光譜反射率與非植被地區的光譜反射率之間具有明顯的差異,很多學者根據植被的光譜響應特征,利用不同波段反射率構造植被指數來表征植被狀況[15]。根據本研究參考傳感器波段設置和前人的生物量反演研究,共選擇了18 種植被指數以及無人機5 個單波段光譜反射率用于生物量反演,本研究所用到的植被指數計算公式如表1 所列。

表1 選用于反演的植被指數計算公式Table 1 Calculation formulae of vegetation indexes selected for inversion

2.1.2 紋理特征

紋理特征是物體固有的特性且穩定性較好,很難受到外界因素的影響,灰度共生矩陣法具有旋轉不變性和多尺度特性,常常用于提取紋理信息。本研究采用灰度共生矩陣法共提取5 個波段的8 個紋理特征,包括每個波段的均值(mean)、方差(variance)、同質性(homogeneity)、對比度(contrast)、非相似性(dissimilarity)、熵(entropy)、角二階矩(angular second moment)、相關性(correlation)[26-27]。提取出的紋理特征及其計算公式如表2 所列。

表2 紋理特征及其計算公式Table 2 Texture features and their calculation formulae

2.1.3 地形特征

自然地形影響降水和陽光輻射的空間分布,從而導致土壤的質地、養分和含水量分布差異,進而決定植被生物量分布以及生長更新[28]。由此本研究引入坡度、坡向等地形特征建立生物量反演模型,實現礦區植被生物量準確監測。在ArcGIS 10.6 中利用實測高程點數據生成不規則三角網(triangulatedirregular network, TIN)后轉成高精度DEM 柵格數據,并生成坡度、坡向、剖面曲率和平面曲率4 種地形數據,如圖3 所示。

圖3 研究區地形特征Figure 3 Topographic characteristics of the study area

2.2 變量篩選

經以上數據處理后,共獲取67 個特征變量,其中光譜特征變量23 個,紋理特征40 個,地形特征4 個。本研究先后采用灰度關聯法和熵權法對光譜特征和紋理特征分別進行篩選,由于地形特征較少,則全部作為輸入變量。

灰度關聯法是判斷變量是否相關并確定其相關程度的一種分析方法,其根據各變量間的相對變化趨勢表征變量間的關聯程度,如果變量間相對變化趨勢越一致,則該變量間的關聯度越大,其基本思想是通過計算關聯度找出系統各因素之間的主次關系,從而找出影響最大的因素[29]。

熵權法是通過信息熵計算變量包含的信息量大小并且對變量進行賦權[30],根據權重大小對變量進行篩選。問題涉及多個變量,需確定每個變量的相對權重,一個變量的信息熵值越小其權重越大,即越能反映真實信息,該方法可以消除人為的干擾,為一種較為客觀的變量篩選方式。在計算完熵權值后,可通過設定閾值的方式對變量進行篩選,即當熵權值大于閾值時,保留該變量,反之則刪除。本研究設定的閾值為0.1。

2.3 反演模型構建

實地共采集43 個生物量樣方數據,在43 個實測生物量樣方數據中隨機選取33 個用于反演模型的訓練,剩余10 個用于模型的驗證,以此作為模型的輸入樣本,構建神經網絡反演模型。將最終篩選出所有特征變量分為光譜特征組和多特征組。選用反向傳播神經網絡(BPNN)、卷積神經網絡(CNN)和Elman 神經網絡分別構建兩組輸入變量的草本植物生物量反演模型,3 種算法均于MATLAB 2022a軟件中實現。

反向傳播神經網絡是一種多層前饋神經網絡,該機器學習方法的特點是信號向前傳遞,誤差反向傳播,因其具備非線性映射能力且適用于具有復雜內部機制的數據問題,并具有自學習和自適應能力以及泛化能力[31-32]。本研究依據輸入變量個數設置神經網絡輸入層個數,輸出層為1,對模型進行多次訓練調試后最終確定隱含層為10,激活函數為“tansig”函數,學習率為0.01,迭代次數為1000。

卷積神經網絡是具有卷積結構的深度神經網絡,減少了網絡參數和過擬合問題,卷積神經網絡具有學習性可以構建高精度作物估產模型[33]。卷積神經網絡由輸入層、卷積層、池化層、全連接層和輸出層組成。其具有的局部連接、參數共享和下采樣降維3 個核心思想,使其可以加快網絡收斂和計算速度,去除冗余特征。本研究采用了2 個卷積層、2 個池化層和一個全連接層構建CNN 網絡結構,其中卷積核個數為32,大小設置為2 × 2,步長設為1;池化核尺寸為3 × 3,步長為1,池化方式為最大池化;卷積層和全連接層由ReLU 激活函數激活。

Elman 神經網絡是一種典型局部回歸神經網絡,可以被看作是具有局部記憶單元和反饋連接的遞歸神經網絡,其較好地適應動態變化的數據[34];Elman 神經網絡主要由輸入層、隱含層、上下文層和輸出層組成,上下文節點的自連接使其對歷史輸入數據較為敏感。本研究選擇sigmoid 函數作為隱含層的傳遞函數,輸出層采用purelin 傳遞函數,隱含層節點數設為10。網絡訓練函數為traingdm,其具有較快的收斂速度。

2.4 模型精度評價方法

本研究采用決定系數(R2)、均方根誤差(RMSE)兩個指標來評價模型的精度。其中R2值越接近1 表明模型預測精度越高,RMSE 值越小表明模型精度越高,預測值與實測值越接近。計算公式如下:

式中:yi為實測值,y?i為預測值,yˉ為實測值的均值。

3 結果與分析

3.1 變量篩選

本研究首先利用灰色關聯度方法分別對光譜特征和紋理特征進行處理,并根據關聯度大小進行從大到小排序,且保留排序前10 的光譜特征和紋理特征;然后使用熵權法計算所保留變量的熵權值,以設定的閾值0.1 對變量進行二次篩選,最終篩選出5 個光譜特征和4 個紋理特征,計算結果如表3 所列。將篩選結果與地形特征相結合,從而構建出兩個生物量反演模型的輸入變量組:將篩選出的光譜特征變量作為模型輸入變量,即為光譜特征組;在光譜特征變量基礎上同時引入紋理特征和地形特征變量,即為多特征組。

表3 本研究輸入模型變量Table 3 Model input variables used in this study

3.2 基于光譜特征組的礦區草本植物生物量反演

利用無人機多光譜影像提取采樣點的單波段光譜反射率,并構建相關植被指數,將5 個波段光譜反射率和構建的18 種植被指數作為光譜特征與實測地上生物量進行相關程度分析,先后使用灰色關聯法和熵權法選出相關程度較高且信息熵較大的光譜特征變量建立草本植物生物量反演模型?;诠庾V特征分別構建BPNN、CNN 和Elman 神經網絡的生物量反演模型,其反演結果如圖4 所示,可以看出,基于BPNN 構建的生物量反演模型的決定系數(R2)為0.762,均方根誤差(RMSE)為14.459 g·m-2;CNN 構建的反演模型R2為0.540,RMSE 為20.119 g·m-2;Elman 神經網絡反演模型R2為0.763,RMSE為14.442 g·m-2。結果表明,輸入相同的光譜特征,3 個模型表現出不同的反演精度,其中Elman 神經網絡模型反演效果最佳,其次為BPNN 模型,CNN模型反演精度最低。

圖4 基于光譜特征組反演模型精度驗證Figure 4 Accuracy verification of inversion models based on spectral features

3.3 基于多特征組的礦區草本植物生物量反演

根據表3 所篩選出的5 個光譜特征,4 個紋理特征以及4 個地形特征共13 個變量,分別構建BPNN、CNN 和Elman 神經網絡的生物量反演模型,反演結果如圖5 所示,可以看出,BPNN反演模型R2為0.841,RMSE 為11.813 g·m-2,CNN 反 演 模 型R2為0.683,RMSE 為16.709 g·m-2,Elman 神經網絡模型R2為0.814,RMSE 為12.784 g·m-2。其中,所建立的BPNN 反演模型的預測精度最高。

圖5 基于多特征組反演模型精度驗證Figure 5 Accuracy verification of inversion models based on multiple features

3.4 模型的綜合評價

3.4.1 模型精度對比

以單波段光譜和植被指數作為光譜特征組,以光譜特征、紋理特征和地形特征作為多特征組,分別用作輸入變量構建基于BPNN、CNN、Elman 神經網絡的草地地上生物量反演模型,對構建的生物量反演模型進行精度驗證,各特征組的模型反演結果如表4 所列。其中,在以光譜特征組作為模型輸入變量時,BPNN 和Elman 算法模型精度較為接近,CNN 算法精度最低;以多特征組作為輸入變量時,BPNN 算法精度最高,R2為0.841,RMSE 為11.813 g·m-2。根據表4 可得,構建的6 個生物量反演模型R2均在0.5 及以上,均表現出良好的反演精度。對比光譜特征組,基于多特征組構建的反演模型呈現出較優的性能,紋理特征變量和地形特征變量的加入使3 種機器學習算法模型都有不同程度的提高,相比于光譜特征組,基于多特征組BPNN 算法構建的生物量反演模型R2提高了10.37%,CNN 算法構建的反演模型R2提高了26.48%,Elman 神經網絡構建的模型的R2提升幅度為6.68%。對比6 個生物量反演模型結果精度可得,基于多特征組的BPNN 草地地上生物量反演模型表現出最優性能,模型預測精度最高,因此本研究采用該模型對研究區進行全區生物量反演。

表4 基于不同特征組的算法模型Table 4 Algorithm models based on different feature groups

3.4.2 模型穩定性分析

單獨分組構建的草本植物地上生物量模型的驗證集結果顯示BPNN 算法性能較好,且基于多特征組構建的模型精度較單特征有明顯的提升。為進一步評價模型的穩定性,分別以兩組特征作為模型輸入變量,將樣本數據分為6 組,對各模型進行交叉驗證,結果如表5 所列,可以看出,以均值分析,對于相同的輸入變量組,BPNN 構建的模型R2、RMSE 均為最優,其次為Elman 模型,CNN 模型精度最差;從標準差結果來看,以光譜特征組作為輸入變量組時,CNN 模型精度交叉驗證結果較為集中,而BPNN 結果較為分散,而多特征作為輸入變量時,BPNN 穩定性顯著提升,且穩定性最優。結合最優模型精度綜合分析可得,BPNN 擬合精度最優,且相較于光譜特征組,多特征組提高了3 種模型的擬合精度,且提升了BPNN 模型預測結果的穩定性。

表5 基于不同特征組的模型交叉驗證Table 5 Model cross-validation based on different feature groups

3.5 反演模型應用

為獲取全研究區草本植物生物量分情況,使用訓練好的多特征BPNN 神經網絡模型對全區草本植物生物量進行反演,反演所得的生物量分布情況如圖6 所示,其中區域為紅色表示該區域生物量高,區域為綠色表示該區域生物量低。由圖6 可以看出研究區生物量主要分布于0~100 g·m-2。將區域生物量分為5 個等級,如表6 所示,全區生物量多集中于第2 等級,即20~40 g·m-2,生物量低于60 g·m-2的區域的面積占全區87.39%,說明研究區大部分區域的生物量處于較低等級。反演結果與實地采樣情況大致相同。

表6 研究區生物量分級Table 6 Biomass classification in the study area

圖6 研究區生物量分布圖Figure 6 Biomass distribution map of the study area

4 討論

修復礦區地表的植被生物量可以反映礦區生態恢復情況,故實現對礦區植被生物量精確預測有利于修復礦區生態監測工作。本研究在選用光譜特征作為生物量反演模型輸入變量的基礎上,引入紋理特征變量和地形特征變量,并采用BPNN、CNN 和Elman 神經網絡3 種算法構建草本植物生物量反演模型,比較模型性能后選用基于多特征的BPNN 反演模型對全區生物量進行了反演。研究結果可為無人機遙感礦區修復地草本植物生物量反演研究提供參考。

本研究首先分別對光譜特征和紋理特征使用灰色關聯度和熵權法進行篩選,光譜特征中單光譜波段篩選結果為紅邊波段和近紅波段,說明紅邊波段和近紅波段對估算生物量的貢獻較大,這與李淑貞等[35]的研究篩選結果相同;植被指數篩選結果為綠度指數、比值指數和非線性植被指數,這些植被指數主要為綠波段、紅波段以及近紅波段的組合,同時發現篩選出的紋理特征也是來源于這些波段,這與劉暢等[5]紋理指標篩選結果一致,由此本研究篩選結果可為植被生物量反演相關研究提供參考。

從所構建的生物量反演模型精度可以看出,引入紋理信息和地形信息后,3 種機器學習模型的反演精度都有相應程度的提升,這主要的原因是紋理信息增加了地上植被的空間特性,而地形差異直接影響著植被分布和生長狀態,引入紋理特征和地形特征后同時考慮了3 種特征對生物量的貢獻性,解決了光譜特征的單一光譜信息不足以及光譜飽和問題[36-37]。其中基于CNN 生物量反演模型精度提升幅度最大,但是精度依然是最低,其原因是CNN訓練需要大量的樣本數據,而本研究樣本量過少,使其沒有得到充分的訓練。

本研究基于兩組輸入變量使用3 種機器學習算法構建生物量反演模型并進行反演精度比較,其中以多特征組為輸入變量構建的BPNN 反演模型精度最高,相同的輸入變量基于BPNN 構建的反演模型也表現出較優的預測性能。本研究結果與以往的研究結果一致,Lyu 等[38]以內蒙古草原為例,分別構建單因素參數模型、多因素參數模型和多因素非參數模型反演草地地上生物量,結果表明BPNN 反演模型的準確性和穩定性均高于其他模型;Yang 等[39]基于最大貢獻變量建立BPNN 反演模型,驗證了BPNN 反演模型比傳統多變量回歸模型有更好的性能。其主要原因是BPNN 較其他機器學習算法更具有魯棒性和容錯性且適用于解決復雜的非線性關系。

修復礦區草本植被生長狀態可以由多種長勢指標反映,如葉片葉綠素含量、植株含水率等;本研究僅以生物量作為單一指標反映礦區植被長勢,未來可以綜合多個指標實現對修復礦區植被監測;此外本研究僅使用了3 種機器學習算法進行反演模型的構建,未來可以選擇更多的算法并對算法進行優化以實現針對特定區域構建更高精度的反演模型。

5 結論

1) 對光譜特征和紋理特征變量分別使用灰度關聯法和熵權法進行初次篩選和二次篩選,最終確定光譜特征變量為紅邊、近紅外波段、綠度指數、比值植被指數和非線性植被指數,紋理特征變量為var-840、con-840、con-730 和var-560,地形特征變量包括坡度、坡向、平面曲率和剖面曲率。

2) 在單波段光譜和植被指數兩光譜特征的光譜特征基礎上引入紋理特征和地形特征,并劃分為光譜特征組和多特征組,將其作為模型輸入變量構建BPNN、CNN 和Elman 神經網絡反演模型,其中基于多特征組構建的3 種反演模型精度均高于基于光譜特征組構建的反演模型,且BP 神經網絡模型精度最高,該模型R2為0.841,RMSE 為11.813 g·m-2,并對3 種反演模型進行了交叉驗證,進一步表明了基于多特征構建的BPNN 草本植物地上生物量反演模型更加穩定,反演精度最優。

3) 基于多特征構建的BPNN 生物量反演模型對研究區生物量分布進行預測分析可得,全區草本植物生物量估測值范圍為0~100 g·m-2、20~40 g·m-2的所占比重較大,占全區面積的52.9%,繪制結果與實地調查情況較為相符,研究結果可實現對草地地上生物量快速、準確的預測。

猜你喜歡
植被指數波段反演
反演對稱變換在解決平面幾何問題中的應用
AMSR_2微波植被指數在黃河流域的適用性對比與分析
河南省冬小麥產量遙感監測精度比較研究
基于低頻軟約束的疊前AVA稀疏層反演
基于自適應遺傳算法的CSAMT一維反演
M87的多波段輻射過程及其能譜擬合
日常維護對L 波段雷達的重要性
主要植被指數在生態環評中的作用
基于MODIS數據的植被指數與植被覆蓋度關系研究
基于SPOT影像的最佳波段組合選取研究
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合