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

基于集合經驗模態分解算法的土壤肥力脅迫甄別與監測

2022-02-15 01:21:22李旭青王小丹趙辰雨張文龍王春暖
農業工程學報 2022年21期
關鍵詞:研究

李旭青,劉 帝,王小丹,趙辰雨,張文龍,王春暖

基于集合經驗模態分解算法的土壤肥力脅迫甄別與監測

李旭青1,2,劉 帝3,王小丹1,2,趙辰雨1,2,張文龍1,2,王春暖1,2

(1. 北華航天工業學院遙感信息工程學院,廊坊 065000;2.河北省航天遙感信息處理與應用協同創新中心,廊坊 065000;3.航天宏圖信息技術股份有限公司,北京 100195)

土壤肥力是農作物生長所需的關鍵要素之一,其水平直接影響農作物的長勢甚至產量,然而作物生長過程中受多種脅迫因素的綜合影響。避免土壤肥力監測結果受到其他脅迫因素的干擾是土壤肥力精準監測的關鍵問題之一。該研究旨在分析利用集合經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)算法進行土壤肥力脅迫甄別和監測的可行性。以河北省廊坊市大廠縣冬小麥耕地為研究區域,以EEMD算法為基礎,將分解后的本征模函數(Intrinsic Mode Function,IMF)分量按照年際、年間和年內的尺度進行合成,結合不同的時間尺度及脅迫特征剔除土壤水分脅迫、病蟲害脅迫及重金屬脅迫等因素,實現對土壤肥力脅迫的有效甄別提取。結合主成分分析相關方法,將有機質、全氮、有效磷以及速效鉀4項養分指標轉換為3個主成分,初步得到土壤肥力綜合評價簡易模型,后與分解結果進行擬合,構建并實現了土壤肥力綜合水平的定量評價模型。結果表明:1)在年際、年間以及年內3組波動組分中,年際波動組分可以較好地反映研究區內土壤肥力脅迫作用對農作物長勢的影響;2)利用最小二乘法對土壤養分指標測定數據進行線性擬合,擬合結果與原始數據的決定系數達到了0.857,能夠較好地反映出原始數據的變化水平;3)最終土壤肥力水平評價模型評價結果與實測結果間的平均誤差為11.82%,表明模型預測結果與實際情況契合程度高,模型反演結果能夠較好地反映研究區的土壤肥力水平。該研究結合EEMD算法與統計學分析實現了土壤肥力脅迫的有效甄別及土壤肥力定量評價模型的構建,為遙感技術在土壤肥力研究領域的應用提供了參考。

遙感;土壤;肥力;集合經驗模態分解;脅迫甄別;定量分析

0 引 言

對于農作物而言,其生長過程受諸如水分脅迫、病蟲害脅迫、重金屬脅迫以及土壤肥力脅迫等多種脅迫因素的影響[1-2]。其中,土壤肥力在農作物生長過程中發揮著尤為重要的作用[3-4],它為農作物生長提供良好的環境,使農作物達到高產優質的效果,最終提高農作物的經濟效益。土壤肥力脅迫的甄別與評價研究,對于提高糧食產量、保護耕地土壤質量以及改善環境條件都有十分重要的意義。

傳統土壤肥力水平監測方法主要通過實地采樣測量的方式開展,即對采樣點進行農作物的定期長勢監測和土壤樣本的采集化驗[5]。這類方法的精確度能夠得到保證,但其研究過程花費的人力物力成本較高,難以滿足大面積實際生產應用中的需求。也有研究通過利用傳統的地統計分析法[6]對土壤有機質(Soil Organic Matter,SOM)的空間分布進行預測來間接監測土壤肥力水平,但是土壤空間變異顯著的特點使得地統計分析法需要大量的樣本點作為數據支撐才能保證樣點的代表性足夠強,在土壤空間分布的復雜程度更高的區域,地統計分析法的誤差甚至會進一步加大。

近年來,通過遙感技術開展的土壤肥力監測手段[7]具備實時、快速獲取海量信息的能力,在很大程度上改善了傳統方法存在的缺陷[8]。劉煥軍等[9]通過對比多種降噪方法處理后的高分五號衛星影像的SOM反演結果,證實了降噪方法應用高光譜衛星數據處理中的可行性,為遙感影像數據反演SOM提供了參考;焦彩霞等[10]采用600 nm“弓曲差”以及最小二乘法2種方法構建了一種基于高光譜數據的土壤有機質含量估測模型,為衛星遙感數據在SOM研究中的應用提供了可靠思路。遙感監測方法實現了對土壤肥力的實時、快速監測,但在應用過程中依然存在精確度不高等方面的不足,為此,國內外研究從監測指標和監測模型方法2個角度開展了不同程度的研究。一方面,構建新的模型描述農作物在不同土壤條件下的長勢差異是一條可行的方法。Xu等[11-12]基于WorldView-2影像構建了一種2 m分辨率的土壤總氮含量預測模型以及一種基于貝葉斯克里格模型的土壤可交換鉀含量的預測模型;方慧婷等[13]依據WOFOST作物模型對土壤速效養分進行反演,為速效養分時空變異的預測提供了參考。另一方面,部分研究著眼于尋找對土壤養分指標變換足夠敏感的觀測參數以提高監測精度[14]。唐海濤等[15]對土壤類型特征波段進行了提取,利用篩選的預測因子采用隨機森林模型可實現SOM含量的快速預測;Ibrahim等[16]使用歸一化植被指數(Normalized Difference Vegetation Index,NDVI)、裸土指數以及土壤調整型植被指數對土壤有機質含量進行反演,具有較好的反演效果;吳亞鵬等[17]篩選出包括土壤調整型植被指數、修正型紅邊比率在內的多種植被指數,實現了土壤中氮元素含量變化與敏感植被指數間的有效擬合。

綜上,以往研究證實了遙感技術在土壤肥力監測上的可行性,但農田生態系統多源脅迫下耕地土壤肥力監測的關鍵問題在于如何避免監測結果受到其他脅迫因素的干擾,然而卻少有研究進行土壤肥力脅迫的甄別及提取,導致未能充分剝離其他脅迫因素的影響。而集合經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)算法可以將原始信號分解為幾個本征模函數(Intrinsic Mode Function,IMF),通過分析不同IMF的時域和頻域統計參數,形成不同維度的特征向量[18]。相關學者利用EEMD算法分解提取腦電信號特征[19]、閃電及應力波信號特征[20],進行復合故障檢測與識別[21],在長時間尺度上獲取重金屬應力信號特征[22]。還有學者利用EEMD算法和其他算法結合實現共同降噪,如和小波變換[23]、小波閾值等方法[24]聯和去噪,與形態分量分析、廣義S變換、獨立分量分析方法相結合[25-26]實現混合降噪。EEMD算法可以簡化去噪過程,減少模態混疊現象[27-28],為非線性多尺度時序信號的快速處理和不同時間尺度的特征甄別提供了可靠的解決方案[29]。

因此,本研究旨在分析利用EEMD算法進行土壤肥力脅迫甄別和監測的可行性。以河北省廊坊市大廠縣冬小麥耕地為研究區域,采用高分一號WFV影像為數據源獲取冬小麥NDVI長時間序列,利用EEMD算法將序列進行分解后按照年際、年間和年內3種尺度進行合成,剔除水分脅迫、病蟲害脅迫及重金屬脅迫等因素的干擾,降低其他脅迫因素對土壤肥力監測結果的影響。以土壤有機質、全氮、有效磷以及速效鉀含量4項指標為評價因子,結合主成分分析相關方法構建土壤肥力綜合水平的定量評價模型,實現對土壤肥力脅迫的準確甄別提取,以期提高土壤肥力的監測精度。

1 材料與方法

1.1 研究區概況

本研究選擇廊坊市大廠縣為研究區域,冬小麥為目標作物。大廠縣位于廊坊北部,中心經緯度116.98°E、39.88°N,總面積為176.29 km2,其中耕地為11 581.54 hm2,地貌以平原為主,土壤多為潮褐土或潮土。大廠縣冬小麥耕地主要分布在東部地區,即陳府鄉以及大廠鎮地區。此外,由于研究選取的影像數據源分辨率為16 m,為盡可能避免混合像元導致的監測結果誤差,選取16 m×16 m以上的耕地地塊作為采樣點。結合研究區影像與實地考察結果,最終確定位于大廠縣東南部的大廠鎮及陳府鎮的6處采樣區,共23個采樣點,研究區及采樣區位置如圖1所示。

圖1 研究區及采樣區位置圖

1.2 遙感影像數據與處理

高分一號衛星搭載的WFV傳感器影像幅寬800 km,空間分辨率為16 m,覆蓋周期達到41 d,具備高空間和高時間分辨率的優勢且有足夠的寬幅,能夠滿足本研究的數據需求。綜合發射時間等多方面因素考慮,本次研究選取大廠縣自2015年1月—2020年11月,共計102景影像作為研究的影像數據源。高分一號WFV影像的原始影像為1A級影像,只帶有基本的地理坐標信息,經過輻射定標、大氣校正、正射校正、幾何校正以及圖像融合的預處理流程后備用。

1.3 野外實測數據采集

野外實測數據用于表征試驗區域土壤肥力情況,并與EEMD分解提取結果進行擬合。在進行土壤采樣時,額外選取若干相鄰冬小麥耕地進行采樣測定備用數據,以用作研究結果精度驗證。

1.3.1 采樣方法

采集土壤樣品過程的科學性是影響野外實測數據獲取結果準確程度的關鍵,本研究土壤采樣遵循隨機以及等量2個原則。一方面,土壤樣品的采集地點隨機分布在整個采樣區;另一方面,每個采樣地區內土壤樣品的采樣數量應該保持一致。一般情況下,采樣點數量的選擇取決于研究總體范圍的大小、研究對象的復雜程度以及研究要求的精確度等條件。理想情況下,應滿足采樣點盡可能少的前提下,樣品的代表性最大。為滿足以上條件,在本次采樣工作中,采用“X”形或者“S”形2種土壤采樣方法,如圖2所示。

一個地區采樣點的數量在8~10個點之間,土壤樣品的質量在1 kg以上。采樣時,使用管型土鉆、打孔取土器、土鏟等工具,采集深度距離地表0~20 cm之間。此外,采樣時避開田邊、路邊、樹邊、溝邊、堆放過肥料的地點,以確保樣品的代表性。

圖2 土壤樣本采集方式及過程示意圖

1.3.2 采樣指標及處理流程

選取4種具有較強代表性的土壤養分指標作為土壤肥力水平的評價因子,即土壤有機質、全氮、有效磷以及速效鉀含量4項指標。在冬小麥進入越冬期前對選取的耕地樣本點進行土壤樣品采集,并封裝在密封袋中帶回實驗室進行處理。

1)風干:風干過程在通風情況良好的室內進行,將土壤樣品置于鐵盤中,避免曝曬,并避免酸、堿性氣體以及灰塵的污染。當土壤樣品達到半干狀態時,將樣品中體積較大的土塊碾碎,以避免完全風干后板結成塊。

2)磨細過篩:將經風干處理后的土壤樣品放置在平整的平臺上,磨碎成塊的土壤之后用1 mm孔徑的篩子過篩。對于未過篩部分的土粒重新進行壓碎過篩處理,石子、植物殘根等雜質予以剔除。

3)標注送檢:將過篩處理后的土壤樣品重新混合均勻,裝入專用樣品袋中封存。在外封裝上注明土壤樣品的名稱、采集時間、采集地點等條目,并依照編號對所有土壤樣品進行登記,避免日光、高溫、潮濕等環境,送往專業檢測部門進行養分指標含量的測定。

經專業監測部門測定完成后,在各采樣地區取各項土壤養分指標的均值作為最終的測定結果,如表1所示。

表1 冬小麥種植區土壤養分指標測定結果

除對土壤樣品進行采集化驗外,測量土壤中的重金屬元素含量是否正常是保證土壤肥力水平監測精確度的重要環節之一。手持式土壤重金屬分析儀能夠對土壤中的鐵、鋅、銅、汞等多種重金屬元素的含量進行定量測定,測定時保證探頭上的保護膜完整,一旦出現破損情況,立即對保護膜進行更換,以避免造成探頭內部的損壞。采用原位測量法,即不進行額外的土壤采樣,直接使用儀器對耕地土壤進行測量,如圖2c所示。原位測量時使測量探頭緊貼地表土壤,避免外部光線對儀器自身發出的光束造成影響。

1.4 集合經驗模態分解算法

EEMD算法是一種參考數據本身時間尺度特征的時域信號分析方法,適用于分析處理非線性和非平穩信號。通過多次向原始序列中添加不同振幅的噪聲信號最終結果取均值的方法在很大程度上抵消“模態混疊”造成的誤差[22],對于研究單一脅迫作用對農作物長勢的影響具有很大的優勢。

在EEMD算法的實現過程中,加入的噪聲信號的振幅以及計入噪聲信號的次數是影響最終分解結果精度最重要的2個因素。選擇合適的噪聲振幅以及加入白噪聲次數是利用EEMD算法進行時間序列分解的關鍵。理論上,加入噪聲信號的次數越多,EEMD分解的結果越準確,但加入白噪聲次數設置過大會明顯降低效率,同時計算成本也隨之增加。同樣,當噪聲信號的振幅太大時,EEMD分解過程將產生多項缺少有效信息的IMF分量,導致分解信號失去原始信號的特征;當噪聲信號的振幅過小時,其對原始信號產生的影響不足以消除模態混疊現象,同樣無法保證分解結果的精度。通常在EEMD算法分解中,噪聲信號振幅的取值范圍在0.01~0.4之間,加入白噪聲次數的取值為50或100。根據實際情況,當農作物的長勢不受環境脅迫因素影響時,其每年的生長趨勢應基本保持不變,因此可將作物的年間波動組分在不同年間的相關性作為確定參數的基準。經過多次試驗,當加入噪聲信號的振幅相同時,加入白噪聲次數的次數對分解效果的影響很小,為使結果真實性更強,將加入白噪聲次數設置為100;當噪聲信號的振幅取0.3時,作物每年的固有生長趨勢最為相似,因此將所加入噪聲信號的振幅的值設置為0.3。

1.5 模型構建與驗證

耕地土壤的肥力水平是多種土壤養分指標共同作用的結果,在不采用試驗田方法的情況下研究單一某種養分含量變化對農作物長勢的影響是較為困難的。通過合適的方法構建一個能夠代表土壤肥力整體水平的參數能夠有效降低分析過程的難度。在統計學中,主成分分析是考察多個變量間相關性的一種多元統計手段[30],它能夠盡可能地反映原始變量信息的少數幾個主成分,并對幾個原始成分進行線性組合,形成一個綜合指標,從而實現對原始序列的表達。通過對上文提到的4項典型土壤養分指標的主成分分析確定土壤肥力綜合水平與各項指標間的關系。在SPSS環境中進行指標間的相關性檢驗以確定各項指標之間的相關性,相關系數計算方法采用Pearson法。根據相關性檢驗的結果,對土壤肥力水平進行主成分分析,將3個主成分線性組合得到總體關于3個主成分的線性表達式,最終得到土壤肥力整體水平關于4項養分指標的線性表達式。將其結果與EEMD算法的分解結果進行擬合進一步構建土壤肥力綜合水平的定量評價模型,在進行精度驗證時采用野外實測數據采集時選取的采樣點備用實測結果作為驗證數據,根據擬合誤差高低選出最優擬合方式。

2 結果與分析

2.1 土壤肥力綜合評價簡易模型構建

采用Pearson法進行相關分析以確定各項指標之間的相關性關系,根據相關性檢驗的結果,有機質與有效磷的相關性較高、與速效鉀的相關性較低;全氮和有機質的相關性相比與其他2種元素的相關性更高一些,其與速效鉀的相關性最低;有效磷與全氮的相關性最高。對土壤肥力水平進行主成分分析,結果如表2所示,前3項主成分對總體的累計貢獻率已經達到了99.16%,基本能夠反映總體的全部信息。因此,借助前3項主成分對整體進行表達的準確度能夠滿足需求。

表2 主成分分析結果

根據表2可以得出3個主成分關于4項土壤養分指標的線性表達式:

1=0.561+0.322+0.033+0.544(1)

2=?0.141?0.212?0.873+0.424(2)

3=0.791?0.312?0.273?0.454(3)

式中1為有機質含量,g/kg;2為全氮含量,mg/kg;3為有效磷含量,mg/kg;4為速效鉀含量,mg/kg;為第個主成分。

通過3個主成分的線性組合能夠得到總體關于3個主成分的線性表達式,如式(4)所示:

=0.5761+0.3172+0.0993+Δ() (4)

式中表示土壤養分整體水平,Δ()表示土壤養分剩余成分。進一步地,代入每個主成分對總體的貢獻率,即可得到土壤肥力整體水平關于4項養分指標的線性表達式:

=0.3591+0.2882+0.1493+0.2024(5)

結合式(5)及表1數據可以得出大廠縣各樣區的土壤肥力綜合水平,如表3所示。

表3 土壤肥力綜合水平計算結果

2.2 集合經驗模態分解結果的合成

土壤中的重金屬不易遷移、抗降解,各種養分的含量通過不同途徑的生物界循環在整體上保持一定幅度內的平衡,土壤重金屬與土壤肥力脅迫因素具有相對持久性和穩定性。與此相反,水分脅迫以及病蟲害脅迫在不同時間和空間上的影響通常更短暫,可能只發生在一個或幾個生長周期,具有易揮發性。土壤肥力脅迫對作物的影響具有相對穩定的年際變化特征,通過對作物的長期連續監測,可以根據不同的時間尺度及頻率特征有效區分土壤肥力脅迫與其他脅迫。

NDVI時間序列經過EEMD分解為了6個IMF分量,結合上述脅迫因素特征分析以及頻率差異將所有IMF分為2項周期大于年的分量、1項周期約等于年的分量、2項周期小于年的分量以及1個殘差項,將這些分量根據不同時間尺度特征劃分為年際波動組分、年間波動組分、年內波動組分以及殘差項。選取幾處典型冬小麥耕地采樣點為例,對EEMD分解結果的各組分進行甄別。

2.2.1 年內波動組分甄別

EEMD分解結果中的分量IMF1以及IMF2的頻率較高,二者的周期均小于1 a,可以判定IMF1以及IMF2為信號的年內波動組分,這一部分分量包含代表植物生長周期中受到的瞬時脅迫影響作用和信號處理過程中產生的誤差結果。其中,IMF1具有明顯的不規則性,屬于信號處理過程產生的誤差噪聲信號,這一組分不參與模型最終分解結果的合成。因此,IMF2代表了冬小麥長勢長時間序列分解結果中的年內波動組分,如圖3所示。

注:IMFj為第j個本征模函數,下同。

如圖3所示,其頻率分辨率較低,能夠對原始序列的高頻振蕩部分進行很好地表現,可以檢測信號夾帶的瞬時異?,F象。這些信息的特征表現為年內波動雜亂,振幅規律性不強。對相鄰年份同一月份的信號數值進行相關性分析能夠較好地對這一特征進行驗證,選擇3月、4月、5月、6月共4個冬小麥生長較為旺盛的月份進行相關性分析,結果如表4所示。從表中結果可以看出,相鄰年份間的同一月份數值相關性較小,表明曲線的規律性較弱。因此,可以認為冬小麥長勢長時間序列的年內波動組分代表了冬小麥受各種瞬時脅迫因素的作用結果。

表4 年內波動組分相關性分析結果

表5給出了廊坊市大廠縣2015年至2019年的年降水量。其中,2018年的年降水量遠高于其他年份,較多的降水和連續陰雨天氣會對冬小麥的長勢和最終產量造成影響。研究區降水主要集中于每年的7—9月,結合該地區冬小麥耕作“春灌1水”的灌溉方式,水分脅迫作用應在上半年和下半年各出現1個峰值。在圖3代表年內波動組分的曲線中,2018年的曲線波動幅度明顯小于其他年份,并且曲線在每年的4月以及9月左右出現了2個波峰。這些特性均符合水分脅迫對冬小麥長勢的影響,進一步證實了分解結果中瞬時脅迫因素的甄別有效性。因此,可以認為信號的高頻分量即年內波動組分能夠代表冬小麥受水分、病蟲害等瞬時脅迫作用影響的結果。

表5 大廠縣年降水量

2.2.2 年間波動組分甄別

為準確甄別出研究區域內冬小麥長勢的年間波動組分,在此引入季節平均模型方法[31]。對長時間序列所包含年份中的每年同月的數據求取平均值,并重復6 a以獲取平均季節數據。此平均季節數據代表冬小麥長時間序列中所有周期特征的平均波動結果,能夠反映不同年份間冬小麥生長所共有的周期特征。利用EEMD算法對平均季節數據進行分解,得到代表年間波動的IMF分量,如圖4所示。圖4a為基于研究中長時間序列每年同月數據平均值計算得到的平均季節數據,其曲線周期約為1 a;圖4b為EEMD分解結果中,曲線波動周期與平均季節數據結果最為接近的IMF分量,即季節項分量。根據季節平均模型方法,此分量即為代表冬小麥生長過程中年間波動組分脅迫作用的結果。

圖4 平均季節模型示意圖

對原始序列EEMD分解結果中所有與季節項IMF分量周期一致的分量進行篩選,得到冬小麥長勢的年間波動組分,如圖5所示。年間波動組分代表了冬小麥植株在不受外界條件影響下的生長趨勢,表現為一個近似年周期性的序列。在實際研究中,耕地作物類型的變更以及試驗過程中數據獲取的誤差會導致年間波動組分在不同年份存在一定的差異,這種差異主要體現在波峰與波谷的數值差異上,但在時間尺度上峰值對應的時間點在不同年份內大體是一致的。因此,年間波動組分能夠代表冬小麥固有生長趨勢的作用結果。

圖5 年間波動組分示意圖

2.2.3 年際波動組分甄別

年際波動組分又稱為長時間序列分解結果的趨勢項,它代表了信號在一個長時間尺度上的振動趨勢。對于農作物而言,趨勢項代表的意義是其受到時間尺度大于季節項的脅迫因素影響時所產生的響應。因此,將振動頻率小于年頻率的若干IMF分量進行合成,得到冬小麥長勢長時間序列分解結果中的年際波動組分,結果如圖6所示。從曲線可以看出,年際脅迫是長期而平緩的。土壤肥力脅迫的作用形式主體為農作物對土壤中養分的吸收,土壤中各種養分的含量通過不同途徑的生物界循環在整體上是保持一定幅度內的平衡的。土壤中的重金屬不易遷移、抗降解,重金屬脅迫對作物的脅迫具有相對的持久性和穩定性。這些特征符合年際波動組分曲線平穩且振幅較小的特點,因此,可以認為年際波動組分中包含著土壤重金屬脅迫和土壤肥力脅迫的作用結果。

圖6 年際波動組分示意圖

結合上述分析,已將基于EEMD算法分解所得到的6個IMF分量甄別劃分為年內、年間以及年際3種波動組分,排除年內以及年間波動組分可有效剔除研究區冬小麥所受水分脅迫、病蟲害脅迫以及固有生長趨勢的影響。

將野外實測工作中利用手持式土壤重金屬分析儀對采樣區耕地土壤中重金屬元素含量的測量數據與標準含量參考范圍進行對比,所有樣區耕地土壤中的重金屬元素含量均在正常水平范圍內,能夠確定幾個樣區不存在重金屬污染的情況,該地區未受到重金屬脅迫作用的影響。因此,可以認為本研究中的年際波動組分表征的是研究區域內土壤肥力脅迫作用對于農作物長勢的影響。將EEMD算法的分解結果與實測數據分析相結合,剔除了水分脅迫、病蟲害脅迫以及土壤重金屬脅迫,實現了土壤肥力脅迫的準確甄別與提取。

2.3 基于EEMD分解結果的土壤肥力綜合評價模型構建

為實現土壤肥力水平的宏觀監測,將土壤肥力綜合評價模型與EEMD方法的分解結果進行擬合。研究中將23個采樣點的采樣數據按照一定比例分為建模集與驗證集,其中,17個點的采樣數據用作建模集,另外6個點的采樣數據用作驗證集。表6列出了17個建模集對應采樣點的土壤肥力水平與EEMD分解的年際波動組分數據,其中,IMF10表示10月的土壤肥力脅迫量,IMF11表示11月的土壤肥力脅迫量。為避免擬合過程過于復雜,首先利用最小二乘法分別對數據進行一次、二次、三次線性擬合,結果如表7所示。可以看出,雖然IMF11-IMF10比IMF10更能反映冬小麥長勢的變化情況,但其數值變化幅度過小,導致擬合過程中的誤差被大幅度提高。因此,選擇單時相EEMD分解年際波動組分數據更適合對土壤肥力水平進行擬合,結果見式(6)。

0.1686.93×IMF6.31×10-2×IMF20.0002×IMF3(6)

利用最小二乘法對數據進行線性擬合的方法,在實現過程上較為簡單,且擬合結果的決定系數達到了0.857,能夠較好地反映出原始數據的變化水平,實現了土壤肥力水平基于遙感影像數據的定量評價。

表6 模型擬合所用數據

注:IMF10代表10月和11月IMF,此處分別代表10月的土壤肥力脅迫分量,IMF11代表11月的土壤肥力脅迫量,IMF11-IMF10表示11月與10月的土壤肥力脅迫量的差值。

Note: IMF10and IMF11are IMF in October and November, respectively. They represent soil fertilizer stress in in October and November, respectively. IMF11-IMF10is difference between IMF11and IMF10.

表7 不同擬合方式的擬合結果

注:為土壤肥力水平。

Note:is soil fertility level.

2.4 模型驗證

在上文中,利用主成分分析法實現了土壤肥力水平的定量評價,本節將利用前文研究中確定的驗證集對應的6個采樣點的采樣數據進行精度驗證。表8給出了該部分冬小麥耕地4項土壤養分指標測定結果、根據上述擬合模型計算得到的土壤肥力整體水平以及EEMD的擬合結果,并計算了擬合誤差。

結合測定結果以及計算結果來看,平均誤差為11.82%,最高誤差為14.363%,利用EEMD分解得到的土壤肥力水平與實際測量數據結果吻合情況較好。

表8 冬小麥耕地土壤肥力水平模型擬合結果及試驗結果對比

3 討 論

基于EEMD算法依據冬小麥長勢自身的時間尺度特征進行分析,對比已有土壤肥力分析方法具有更強的自適應性,其適合于分析非線性、非平穩信號序列,具有很高的信噪比[32]。利用統計學中相關性檢驗以及主成分分析的方法實現土壤肥力水平的定量反演,在完成對集合經驗模態分解結果的驗證工作的同時,提出了基于土壤有機質含量、全氮含量、有效磷含量以及速效鉀含量4項指標的土壤肥力水平定量評價模型,為使用遙感手段實現大范圍、長時間土壤肥力水平監測提供了一定的參考價值。

水分脅迫、病蟲害脅迫等脅迫對于農作物長勢的影響相較于土壤肥力、重金屬污染等環境因素的響應更快,而年內波動組分與其他組分相比,在曲線上表現為振動頻率更大且振幅更大。這一特點能夠反映出水分脅迫、病蟲害脅迫的脅迫信息。因此,信號的高頻分量可以代表冬小麥受水分、病蟲害等瞬時脅迫作用的響應結果。年間波動組分又可稱為長時間序列分解結果中的季節項,其典型特征是信號組分具有一個明顯的周期。根據農作物的實際物候參數不同,這個周期也會出現一定的變化。通常情況下,這個時間周期的長度是1 a。但由于實際生產過程中耕地的利用方式以及農作物的耕作方式的影響,這個年間波動組分會呈現出周期在1 a左右變動的趨勢。在剔除代表水分脅迫、病蟲害脅迫等瞬時脅迫的年內波動組分以及代表冬小麥固有生長趨勢的年間波動組分后,剩余的年際波動組分可代表持續的環境脅迫對農作物長勢的影響,包括土壤肥力脅迫、重金屬脅迫等因素。同時,根據野外實測工作中利用手持式重金屬分析儀對幾個典型樣區耕地土壤中重金屬元素含量的測定結果,所有樣區耕地土壤中的重金屬元素含量均在正常水平范圍內,說明該地區受到重金屬脅迫作用的影響很弱。因此,本研究中的年際波動組分能夠較好地表征研究區域內土壤肥力脅迫作用對于農作物長勢的影響。

從精度評價的結果可以看出,基于遙感影像提取出的土壤肥力水平總體數值略高于實際測量值,主要原因可能在于:1)受數據源空間分辨率限制導致的“混合像元”現象使得獲取的NDVI值與實際值間存在偏差;2)分解結果中存在除文中提到4項主要脅迫因素外的諸如溫度等其他長期脅迫因素作用的結果,后續研究可重點考慮如何排除此類與土壤肥力脅迫時間尺度差異較小的脅迫因素。

4 結 論

本研究基于高分一號遙感影像數據,以集合經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)算法為基礎,實現了對土壤肥力脅迫的有效甄別提取。結合統計學中主成分分析相關方法,構建并實現了土壤肥力綜合水平的定量評價模型,得到如下結論:

1)將EEMD算法引入研究中,利用EEMD在時間尺度上對信號進行分解的特性,對冬小麥NDVI長時間序列進行了分解,并結合相關脅迫機理特征及野外實測數據對分解結果進行了分析,確定了分別代表冬小麥長勢受年內(瞬時)脅迫、年間脅迫以及年際(長期)脅迫作用的本征模函數(Intrinsic Mode Function,IMF)分量曲線,驗證了EEMD算法在本研究中的可行性;

2)結合野外實測數據的結果構建土壤肥力綜合水平評價模型。利用統計學在中的主成分分析方法,以土壤中有機質、全氮、有效磷、速效鉀4項養分指標為出發點,將4項指標轉換為3個主成分,累計貢獻率達到99.16%,能夠反映出土壤肥力的整體水平,以此實現了土壤肥力水平的初步定量評價;

3)將模型與EEMD分解結果進行擬合,得出以年際波動組分IMF分量作為自變量的土壤肥力水平線性表達式,并與實測數據進行對比分析,最終結果的平均誤差為11.82%,能夠滿足實際需求,實現了基于遙感影像的土壤肥力水平定量評價。

[1] 王林,代永欣,張勁松,等. 水分和光照條件對核桃-黃豆農林復合系統中黃豆光合作用和生長的影響[J]. 林業科學,2020,56(4):188-196.

Wang Lin, Dai Yongxin, Zhang Jinsong, et al. Effects of water and light conditions on photosynthesis and growth of soybean in walnut-soybean agroforestry system[J]. Scientia Silvae Sinicae, 2020, 56(4): 188-196. (in Chinese with English abstract)

[2] 翁亞偉,張磊,張姍,等. 鹽旱復合脅迫對小麥幼苗生長和水分吸收的影響[J]. 生態學報,2017,37(7):2244-2252.

Weng Yawei, Zhang Lei, Zhang Shan, et al. Effects of salt with drought stress on growth and water uptake of wheat seedlings[J]. Acta Ecologica Sinica, 2017, 37(7): 2244-2252. (in Chinese with English abstract)

[3] 黃云鑫,李裕瑞,劉彥隨,等. 土層復配方案對治溝造地新增耕地土壤肥力的影響[J]. 農業工程學報,2021,37(12):64-72.

Huang Yunxin, Li Yurui, Liu Yansui, et al. Effects of soil-layer compounding schemes on the soil fertility of newly-constructed cultivated land[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(12): 64-72. (in Chinese with English abstract)

[4] 范拴喜,崔佳茜,李丹,等. 不同改良措施對設施蔬菜土壤肥力和番茄品質的影響[J]. 農業工程學報,2021,37(16):58-64.

Fan Shuanxi, Cui Jiaqian, Li Dan, et al. Effects of different improvement measures on soil fertility and the tomato quality of facilities vegetables[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(16): 58-64. (in Chinese with English abstract)

[5] 陳紅艷,趙庚星,張曉輝,等. 去除水分影響提高土壤有機質含量高光譜估測精度[J]. 農業工程學報,2014,30(8):91-100.

Chen Hongyan, Zhao Gengxing, Zhang Xiaohui, et al. Improving estimation precision of soil organic matter content by removing effect of soil moisture from hyperspectra[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(8): 91-100. (in Chinese with English abstract)

[6] 劉艷芳,宋玉玲,郭龍,等. 結合高光譜信息的土壤有機碳密度地統計模型[J]. 農業工程學報,2017,33(2):183-191.

Liu Yanfang, Song Yuling, Guo Long, et al. Geostatistical models of soil organic carbon density prediction based on soil hyperspectral reflectance[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(2): 183-191. (in Chinese with English abstract)

[7] 舒田,熊康寧,陳麗莎,等. 基于遙感的石漠化治理下土壤肥力變化特征分析[J]. 生態環境學報,2019,28(4):776-785.

Shu Tian, Xiong Kangning, Chen Lisha, et al. Analysis on variation characteristics of soil fertility under rocky desertification control based on remote sensing [J]. Ecology and Environmental Sciences, 2019, 28(4): 776-785. (in Chinese with English abstract)

[8] 劉煥軍,趙春江,王紀華,等. 黑土典型區土壤有機質遙感反演[J]. 農業工程學報,2011,27(8):211-215.

Liu Huanjun, Zhao Chunjiang, Wang Jihua, et al. Soil organic matter predicting with remote sensing image in typical black soil area of Northeast China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2011, 27(8): 211-215. (in Chinese with English abstract)

[9] 劉煥軍,鮑依臨,孟祥添,等. 不同降噪方式下基于高分五號影像的土壤有機質反演[J]. 農業工程學報,2020,36(12):90-98.

Liu Huanjun, Bao Yilin, Meng Xiangtian, et al. Inversion of soil organic matter based on GF-5 images under different noise reduction methods[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(12): 90-98. (in Chinese with English abstract)

[10] 焦彩霞,鄭光輝,解憲麗,等. 可見-短近紅外成像光譜數據的土壤有機質含量估算[J]. 光譜學與光譜分析,2020,40(10):3277-3281.

Jiao Caixia, Zheng Guanghui, Xie Xianli, et al. Prediction of soil organic matter using visible-short near-infrared imaging spectroscopy[J]. Spectroscopy and Spectral Analysis, 2020, 40(10): 3277-3281. (in Chinese with English abstract)

[11] Xu Y M, Smith S E, Grunwald S, et al. Estimating soil total nitrogen in smallholder farm settings using remote sensing spectral indices and regression kriging[J]. Catena, 2017, 163: 111-122.

[12] Xu Y M, Smith S E, Grunwald S, et al. Evaluating the effect of remote sensing image spatial resolution on soil exchangeable potassium prediction models in smallholder farm settings[J]. Journal of Environmental Management, 2017, 200: 423-433.

[13] 方慧婷,蒙繼華,程志強. 基于遙感與作物模型的土壤速效養分時空變異分析[J]. 中國農業科學,2019,52(3):478-490.

Fang Huiting, Meng Jihua, Cheng Zhiqiang. Spatio-temporal variability of soil available nutrients based on remote sensing and crop model[J]. Scientia Agricultura Sinica, 2019, 52(3): 478-490. (in Chinese with English abstract)

[14] 趙瑞,崔希民,劉超. GF-5高光譜遙感影像的土壤有機質含量反演估算研究[J]. 中國環境科學,2020,40(8):3539-3545.

Zhao Rui, Cui Ximin, Liu Chao. Inversion estimation of soil organic matter content based on GF-5 hyperspectral remote sensing image[J]. China Environmental Science, 2020, 40(8): 3539-3545. (in Chinese with English abstract)

[15] 唐海濤,孟祥添,蘇循新,等. 基于CARS算法的不同類型土壤有機質高光譜預測[J]. 農業工程學報,2021,37(2):105-113.

Tang Haitao, Meng Xiangtian, Su Xunxin, et al. Hyperspectral prediction on soil organic matter of different types using CARS algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(2): 105-113. (in Chinese with English abstract)

[16] Ibrahim M, Ghanem F, Al-Salameen A, et al. The estimation of soil organic matter variation in arid and semi-arid lands using remote sensing data[J]. International Journal of Geosciences, 2019, 10(5): 13.

[17] 吳亞鵬,賀利,王洋洋,等. 冬小麥生物量及氮積累量的植被指數動態模型研究[J]. 作物學報,2019,45(8):1238-1249.

Wu Yapeng, He Li, Wang Yangyang, et al. Dynamic model of vegetation indices for biomass and nitrogen accumulation in winter wheat[J]. Acta Agronomica Sinica, 2019, 45(8): 1238-1249. (in Chinese with English abstract)

[18] Nguyen V H, Cheng J S, Yu Y, et al. An architecture of deep learning network based on ensemble empirical mode decomposition in precise identification of bearing vibration signal[J]. Journal of Mechanical Science and Technology, 2019, 33(1): 41-50.

[19] 宋軻,杜世昌,奚立峰,等. 基于集合經驗模態分解和模糊集的異質信號融合方法研究[J]. 機械設計與研究,2017,33(6):15-20.

Song Ke, Du Shichang, Xi Lifeng, et al. Heterogeneous data fusion based on ensemble empirical mode decomposition and vague sets[J]. Machine Design & Research, 2017, 33(6): 15-20. (in Chinese with English abstract)

[20] Fang Y M, Feng H L, Li J, et al. Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[J]. Sensors, 2011, 11(8): 7554-7567.

[21] 施杰,伍星,劉韜. 采用HHT算法與卷積神經網絡診斷軸承復合故障[J]. 農業工程學報,2020,36(4):34-43.

Shi Jie, Wu Xing, Liu Tao. Bearing compound fault diagnosis based on HHT algorithm and convolution neural network[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(4): 34-43. (in Chinese with English abstract)

[22] Lingwen T, Xiangnan L, Biyao Z , et al. Extraction of rice heavy metal stress signal features based on long time series leaf area index data using ensemble empirical mode decomposition[J]. International Journal of Environmental Research and Public Health, 2017, 14(9): 1018-1035.

[23] 劉慶杰,黃輝,雷曉燕. 基于集合經驗模態分解和小波變換的輪軌力應變信號降噪[J]. 城市軌道交通研究,2016,19(11):26-29,37.

Liu Qingjie, Huang Hui, Lei Xiaoyan. Wheel/rail strain signals de-noising by using EEMD and wavelet transforill[J]. Urban Mass Transit, 2016, 19(11): 26-29, 37.

[24] 李一博,劉嘉瑋,芮小博,等. 基于集合經驗模態分解和小波閾值的真空泵振動信號降噪方法[J]. 航天器環境工程,2019,36(5):450-457.

Li Yibo, Liu Jiawei, Rui Xiaobo, et al. De-noising based on EEMD and wavelet threshold for vacuum pump vibration signals[J]. Spacecraft Environment Engineering, 2019, 36(5): 450-457. (in Chinese with English abstract)

[25] 趙偉,肖世校,張保燦,等. 基于形態分量分析和集合經驗模態分解的心電信號工頻干擾消除法[J]. 生物醫學工程學雜志,2015,32(6):1179-1184.

Zhao Wei, Xiao Shixiao, Zhang Baocan, et al. Removal algorithm of power line interference in electrocardiogram based on morphological component analysis and ensemble empirical mode decomposition[J]. Journal of Biomedical Engineering, 2015, 32(6): 1179-1184. (in Chinese with English abstract)

[26] Liu B, Feng Y Q, Yin Z H, et al. Radar signal emitter recognition based on combined ensemble empirical mode decomposition and the generalized S-Transform[J]. Mathematical Problems in Engineering, 2019, 12: 1-15.

[27] 李明,趙迎,崔飛鵬,等. 基于集合經驗模態分解的拉曼光譜信號特征研究[J]. 光譜學與光譜分析,2020,40(1):54-58.

Li Ming, Zhao Ying, Cui Feipeng, et al. Research on raman spectral signal characteristics based on ensemble empirical mode decomposition[J]. Spectroscopy and Spectral Analysis, 2020, 40(1): 54-58. (in Chinese with English abstract)

[28] 孔麗,劉美玲,劉湘南,等. 利用作物生長模型和時序信號甄別水稻鎘脅迫[J]. 農業工程學報,2021,37(4):249-256.

Kong Li, Liu Meiling, Liu Xiangnan, et al. Identifying heavy metal (Cd) stress in rice using time-series signals and crop growth model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(4): 249-256. (in Chinese with English abstract)

[29] 楊航,朱永利. 基于Storm的局部放電信號集合經驗模態分解[J]. 計算機工程與應用,2020,56(10):261-267.

Yang Hang, Zhu Yongli. Ensemble empirical mode decomposition of partial discharge signal based on storm[J]. Computer Engineering and Applications, 2020, 56(10): 261-267. (in Chinese with English abstract)

[30] 單薇,金曉斌,孟憲素,等. 基于多源遙感數據的土地整治生態環境質量動態監測[J]. 農業工程學報,2019,35(1):234-242.

Shan Wei, Jin Xiaobin, Meng Xiansu, et al. Dynamical monitoring of ecological environment quality of land consolidation based on multi-source remote sensing data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(1): 234-242. (in Chinese with English abstract)

[31] 孫亞楠,李仙岳,史海濱,等. 基于多源數據融合的鹽分遙感反演與季節差異性研究[J]. 農業機械學報,2020,51(6):169-180.

Sun Yanan, Li Xianyue, Shi Haibin, et al. Remote sensing inversion of soil salinity and seasonal difference analysis based on multi-source data fusion[J]. Transactions of the Chinese Society for Agricultural Machinery, 2020, 51(6): 169-180. (in Chinese with English abstract)

[32] 鄧凱,丁建麗,楊愛霞,等. EEMD在土壤剖面反射光譜消噪中的應用[J]. 光譜學與光譜分析,2015,35(1):162-166.

Deng Kai, Ding Jianli, Yang Aixia, et al. EEMD de-noising of reflecting spectrum in soil profiles[J]. Spectroscopy and Spectral Analysis, 2015, 35(1): 162-166. (in Chinese with English abstract)

Screening and monitoring of soil fertility stress using ensemble empirical mode decomposition

Li Xuqing1,2, Liu Di3, Wang Xiaodan1,2, Zhao Chenyu1,2, Zhang Wenlong1,2, Wang Chunnuan1,2

(1.,,065000,;2.,065000,; 3.100195,)

Soil fertility is one of the most important elements required for crop growth. The specific level of soil fertility directly affects the growth and even the yield of crops. However, the development of crop growth depends mainly on a variety of stress factors. Therefore, it is very necessary to avoid the interference of stress factors in the monitoring of soil fertility. The purpose of this study was to accurately identify and monitor the soil fertility stress using Ensemble Empirical Mode Decomposition (EEMD). The research area was taken as the winter wheat cultivated land in the central latitude and longitude of 116.98°E and 39.88°N of the Dachang Hui Autonomous County, Langfang City, Hebei Province, China. The total area was 176.29 km2, among which the cultivated land accounted for 11 581.54 hm2. The landform was mainly plain with mostly tidal brown or tidal soil. The decomposed Intrinsic Mode Function (IMF) components were synthesized using EEMD, according to the annual, inter-annual and intra-annual scales. Different time scales and stress characteristics were combined to eliminate the stress of soil water, pest, and heavy metal, in order to obtain the effective screening and extraction of soil fertility stress. The organic matter, total nitrogen, available phosphorus, and available potassium were also collected from the field soil samples in each planting area of winter wheat. A principal component analysis was performed on the four nutrient indexes to convert them into the three principal components. A simple model was preliminarily obtained for the comprehensive evaluation of soil fertility. Finally, a quantitative evaluation model was established for the comprehensive level of soil fertility after fitting with the decomposition. The results showed that: 1) The components were mainly divided into three groups of fluctuation components after EEMD. Among them, the inter-annual fluctuation component better reflects the effect of soil fertility stress on crop growth in the study area. 2) The principal component analysis of the four indicators showed that the accuracy of the overall expression with the first three principal components fully met the needs, where the cumulative contribution rate reached 99.16%. The Least Square Method (LSM) was used to fit the soil nutrient index data. The correlation coefficient between the fitting and the original data reached 0.857, indicating the better representative for the change level of the original data. 3) The average error between the evaluation and the measurement was 11.82%, indicating that the predicted performance of the model was in better agreement with the actual situation. Therefore, the inversion of the model can be expected to better reflect the soil fertility level in the study area. The effective screening of soil fertility stress was achieved to quantitatively evaluate the soil fertility level using remote sensing image data.

remote sensing; soils; fertility; EEMD; coercion screening; quantitative analysis

10.11975/j.issn.1002-6819.2022.21.017

S158.2

A

1002-6819(2022)-21-0137-10

李旭青,劉帝,王小丹,等. 基于集合經驗模態分解算法的土壤肥力脅迫甄別與監測[J]. 農業工程學報,2022,38(21):137-146.doi:10.11975/j.issn.1002-6819.2022.21.017 http://www.tcsae.org

Li Xuqing, Liu Di, Wang Xiaodan, et al. Screening and monitoring of soil fertility stress using ensemble empirical mode decomposition[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(21): 137-146. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2022.21.017 http://www.tcsae.org

2022-07-08

2022-09-10

河北省青年科學基金項目(D2018409029);河北省高等學??茖W技術研究青年拔尖人才項目(BJ2020056);廊坊市第一批“青年拔尖人才”支持項目(LFBJ202005);河北省高等學??茖W技術研究重點項目(ZD2016126);高分辨率對地觀測系統重大專項?。ㄗ灾螀^)域產業化應用項目(67-Y40G09-9002-15/18);高分共性應用技術規范和高分遙感數據云平臺處理應用共性關鍵技術項目(67-Y20A07-9002-16/17)

李旭青,博士,副教授,研究方向為農業遙感。Email:meililixuqing@163.com

猜你喜歡
研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關于遼朝“一國兩制”研究的回顧與思考
EMA伺服控制系統研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側面碰撞假人損傷研究
關于反傾銷會計研究的思考
焊接膜層脫落的攻關研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 亚洲成人播放| 国产精品亚洲一区二区在线观看| 91在线视频福利| 日本欧美午夜| 亚洲a级毛片| 在线中文字幕网| 日本高清有码人妻| 婷婷六月综合网| 国产精品九九视频| 谁有在线观看日韩亚洲最新视频| 日韩亚洲综合在线| 国产福利一区二区在线观看| 亚洲福利片无码最新在线播放| 国产精品不卡片视频免费观看| 在线高清亚洲精品二区| 中文字幕调教一区二区视频| 成人福利在线看| 国产精品手机视频一区二区| 粉嫩国产白浆在线观看| 国产亚洲精久久久久久久91| 久久99精品国产麻豆宅宅| 2021国产乱人伦在线播放| 亚洲色图欧美| 国产成人精品视频一区二区电影 | 国产成人亚洲欧美激情| 狠狠干综合| 99久久99这里只有免费的精品| 美女内射视频WWW网站午夜 | 92午夜福利影院一区二区三区| 一级毛片免费高清视频| 亚洲资源站av无码网址| 在线观看国产精品第一区免费| 久久国产av麻豆| 欧美激情网址| 最新日韩AV网址在线观看| 在线观看亚洲人成网站| 国产白丝av| 久青草国产高清在线视频| 亚洲人成网7777777国产| 综合久久久久久久综合网| 日韩毛片免费观看| 亚洲国产亚洲综合在线尤物| 暴力调教一区二区三区| 亚洲中文精品人人永久免费| 午夜性刺激在线观看免费| 成人在线视频一区| 日韩一区二区三免费高清| 91精品久久久无码中文字幕vr| 91欧美亚洲国产五月天| 1级黄色毛片| 精品久久蜜桃| 99久久99这里只有免费的精品 | 伊人无码视屏| 国产成人精品亚洲77美色| 国产日韩久久久久无码精品| 99热国产这里只有精品无卡顿"| 免费观看亚洲人成网站| 日韩高清欧美| 成人a免费α片在线视频网站| 一级毛片在线播放免费观看| 亚洲国产天堂在线观看| 中文字幕色在线| 国产导航在线| 中文字幕亚洲无线码一区女同| 国产在线欧美| 日韩精品无码免费专网站| 久久香蕉国产线看观看式| 国产视频只有无码精品| 亚洲资源在线视频| 97se亚洲综合在线韩国专区福利| 国产白浆在线| 久久黄色一级视频| 日韩精品无码免费一区二区三区| 国产成人精品视频一区二区电影 | 亚洲一区二区三区国产精品| 99久久精品免费看国产电影| 亚洲第一区在线| 欧美日韩精品一区二区在线线| 亚洲无码高清免费视频亚洲| 国产波多野结衣中文在线播放| 99视频全部免费| 亚洲av无码牛牛影视在线二区|