• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于分?jǐn)?shù)階微分的土壤重金屬高光譜遙感圖像反演

    2023-10-19 02:48:14丁松滔張霞尚坤李儒孫偉超
    遙感學(xué)報(bào) 2023年9期
    關(guān)鍵詞:模型

    丁松滔,張霞,尚坤,李儒,孫偉超

    1.中國(guó)科學(xué)院空天信息創(chuàng)新研究院,北京 100101;2.中國(guó)科學(xué)院大學(xué),北京 100049;3.自然資源部國(guó)土衛(wèi)星遙感應(yīng)用中心,北京 100048

    1 引言

    成像高光譜遙感能夠獲取圖像上連續(xù)且精細(xì)的土壤反射光譜,具有低成本大范圍快速監(jiān)測(cè)土壤重金屬含量的潛力(Khosravi 等,2018)。然而目前大多數(shù)研究采用實(shí)驗(yàn)室理想條件下采集的光譜,基于高光譜圖像的重金屬反演研究較少(Ding 等,2022),并且由于土壤樣本采集受制于地形、交通可達(dá)性等因素,且土壤樣本的重金屬實(shí)驗(yàn)室化驗(yàn)分析成本高,導(dǎo)致樣本不充足的問(wèn)題普遍存在,而對(duì)于高光譜圖像而言,圖像上的土壤像元數(shù)量與土壤樣本數(shù)相差懸殊,小樣本問(wèn)題更加突出。光譜特征選擇是解決小樣本問(wèn)題的有效途徑之一,而服務(wù)于光譜特征選擇的特征增強(qiáng)方法也是研究熱點(diǎn),研究顯示經(jīng)過(guò)特征增強(qiáng)的光譜曲線比原始光譜曲線可以更有效地提取出光譜的特征波段(沈強(qiáng) 等,2019)。

    微分處理是常用的增強(qiáng)光譜特征的方法,一階微分和二階微分是最常用的光譜微分方式,能夠突出土壤光譜中與重金屬相關(guān)的光譜信息,提高反演精度(郭學(xué)飛 等,2020)。然而,傳統(tǒng)的整數(shù)階微分(即一階、二階微分)缺乏對(duì)可能包含土壤重金屬有益信息的漸變傾斜度或曲率的敏感性(Hong等,2018)。分?jǐn)?shù)階微分作為整數(shù)階微分的擴(kuò)展,對(duì)光譜微分的處理更加細(xì)致,階數(shù)上更精細(xì)的量化意味著能更敏感地捕捉各波段的曲率和傾斜度變化,突出光譜特征。已有學(xué)者將FOD應(yīng)用于實(shí)驗(yàn)室光譜(Chen 等,2022),成功地反演了土壤重金屬Cr、Pb、Zn的含量;Meng等(2021)將FOD應(yīng)用于GF-5圖像,使用FOD處理后的像元光譜計(jì)算光譜指數(shù)作為自變量,構(gòu)建的有機(jī)質(zhì)反演模型精度R2達(dá)到0.84,但FOD 應(yīng)用于高光譜圖像上進(jìn)行土壤重金屬反演的有效性還有待驗(yàn)證。

    特征選擇算法通過(guò)選出特征性較強(qiáng)的波段,提高模型的反演性能?,F(xiàn)有很多研究都已經(jīng)證明,使用優(yōu)選出的特征波段比使用全波段可以獲得更好的反演結(jié)果(Tan 等,2020;Zhang 等,2021),這表明波段選擇的重要性。土壤重金屬反演中最為常見(jiàn)的篩選方法是皮爾森相關(guān)系數(shù)法,但該方法只考慮單個(gè)波段與理化性質(zhì)間的關(guān)系(周冰 等,2021),沒(méi)有考慮多個(gè)波段間的共線性。GA 是一種以“適者生存”為原則的優(yōu)化算法,以最大化精度為目標(biāo),在特征空間中進(jìn)行啟發(fā)式搜索,可以快速獲取近似最優(yōu)解(柏晗 等,2022),已有將其應(yīng)用于遙感圖像上反演重金屬含量的先例(Wang等,2022),然而GA的缺點(diǎn)是會(huì)過(guò)早收斂,導(dǎo)致結(jié)果難以跳出局部解(Pavez-Lazo 和Soto-Cartes,2011)。CARS 同樣是遵循“適者生存”的一種篩選算法(Li 等,2009),通過(guò)逐步去除冗余和不重要變量來(lái)選擇信息,使模型計(jì)算效率與穩(wěn)定性提升(Vohland 等,2017),能夠從全波段數(shù)據(jù)中選擇最優(yōu)的波段組合(Vohland 等,2016),Cheng 等(2021)在土壤全氮反演研究中指出CARS 算法在有效變量的選擇上優(yōu)于GA。綜上所述,CARS 算法在這3 種算法中有利于篩選出更有效的波段組合,達(dá)到全局最優(yōu)解。

    本文以新疆維吾爾自治區(qū)哈密市的黃山南礦區(qū)為研究區(qū),應(yīng)用航空高光譜圖像開(kāi)展對(duì)土壤重金屬Pb、Zn、Ni 含量的反演研究,提出了一種基于FOD 估算高光譜圖像上土壤重金屬濃度的反演方法,通過(guò)擴(kuò)充樣本、FOD 增強(qiáng)光譜特征、CARS篩選波段選取優(yōu)光譜特征,探索高光譜圖像反演中小樣本問(wèn)題的解決途徑,分析確定了研究區(qū)3種重金屬反演的最佳FOD 階數(shù),并將CARS 與CC、GA 的建模結(jié)果進(jìn)行對(duì)比;進(jìn)一步分析了CARS 的波段篩選結(jié)果,得出對(duì)3種重金屬反演貢獻(xiàn)最高的波段范圍及其物理意義,分析了重金屬含量分布圖的可靠性。

    2 研究區(qū)與研究數(shù)據(jù)

    2.1 研究區(qū)概況

    研究區(qū)位置如圖1所示,位于新疆維吾爾自治區(qū)哈密市黃山南銅鎳礦區(qū),是重要的銅、鎳、鐵、鉛、鋅等大型礦床集中區(qū)(王京彬 等,2006),該礦區(qū)于2007年開(kāi)始勘探,至今仍在開(kāi)采中。地處新疆維吾爾自治區(qū)東部,屬于典型的溫帶大陸性干旱氣候,干燥少雨,年降水量33.8 mm,年蒸發(fā)量3300 mm,夏季酷熱、蒸發(fā)強(qiáng)。研究區(qū)內(nèi)主要為空曠戈壁,無(wú)植被覆蓋,土壤表面多為裸露狀態(tài)。

    圖1 研究區(qū)位置示意圖Fig.1 The geographical location of the study area

    2.2 數(shù)據(jù)獲取及預(yù)處理

    航空高光譜圖像獲取于2021 年8 月23 號(hào),傳感器為中國(guó)科學(xué)院上海技術(shù)物理研究所研制的全譜段多模態(tài)成像光譜儀,本次航飛航高3 km,波長(zhǎng)范圍覆蓋350—2500 nm,350—1000 nm 光譜區(qū)間的空間分辨率為0.5 m,有251 個(gè)波段,1000—2500 nm 波長(zhǎng)范圍空間分辨率為1 m,有508 個(gè)波段。研究區(qū)由兩條航帶所覆蓋,圖像經(jīng)過(guò)輻射定標(biāo)、大氣校正后得到地表反射率數(shù)據(jù),對(duì)兩航帶圖像進(jìn)行配準(zhǔn)、拼接等處理后,裁剪出研究區(qū),經(jīng)緯度范圍為:42.197°N—42.216°N,94.625°E—94.684°E。

    成像光譜儀在獲取光譜數(shù)據(jù)時(shí),大氣中的水汽對(duì)1400 nm 和1900 nm 波長(zhǎng)附近的輻射能量存在強(qiáng)吸收,且1900—2500nm 的像元光譜存在較嚴(yán)重噪聲。為去除水汽等噪聲的影響,并保留盡可能多的波段,剔除1350—1450 nm 和1800 nm 之后的波段。對(duì)剩余的光譜進(jìn)行SG 濾波去噪處理,由于不同區(qū)間的光譜噪聲不同,對(duì)其采用不同的窗口大小進(jìn)行濾波,350—1000 nm 區(qū)間的噪聲較弱,采用窗口大小為9 的二次多項(xiàng)式,1000—1800 nm區(qū)間為中度噪聲,采用窗口大小為13 的二次多項(xiàng)式。

    與航空飛行準(zhǔn)同步,2021年7月22日—7月30日開(kāi)展了土壤樣本地面獲取實(shí)驗(yàn),由于研究區(qū)屬于干旱區(qū),且該時(shí)段內(nèi)土壤未受降雨影響,兩者獲取數(shù)據(jù)時(shí)土壤狀況可認(rèn)為一致。土壤樣本采集主要針對(duì)采礦區(qū)及周邊,沿道路兩旁選擇土壤的顏色和粒徑有明顯差異的樣點(diǎn),采集0—20 cm 表層土,共采集土壤樣本72 個(gè)。圖像及土壤樣本點(diǎn)位置見(jiàn)圖2,根據(jù)土壤樣本的經(jīng)緯度信息,提取對(duì)應(yīng)的圖像像元反射光譜數(shù)據(jù)。

    圖2 航空遙感圖像及采樣點(diǎn)位置Fig.2 Aerial remote sensing image and sampling locations

    將土樣在實(shí)驗(yàn)室風(fēng)干研磨后過(guò)100目篩制成標(biāo)準(zhǔn)樣,用火焰原子吸收分光光度法測(cè)定重金屬含量。圖3 所示為樣本中3 種重金屬含量的直方圖,可見(jiàn)重金屬Pb 和Zn 的含量分布都趨向于正態(tài)分布,而重金屬Ni 含量呈現(xiàn)出明顯的偏分布,本文采用對(duì)數(shù)變換對(duì)其進(jìn)行校正,轉(zhuǎn)換后Ni 含量分布如圖4所示,偏分布現(xiàn)象得到明顯改善。

    圖3 土壤樣本的重金屬含量直方圖Fig.3 Histogram of the heavy metal content of soil samples

    圖4 對(duì)數(shù)變化后Ni含量直方圖Fig.4 Histogram of the Ni content after log change

    3 實(shí)驗(yàn)方法

    圖5為本文研究技術(shù)路線圖,首先根據(jù)土壤樣本的經(jīng)緯度提取對(duì)應(yīng)像元和相似的鄰近像元光譜,對(duì)像元光譜數(shù)據(jù)進(jìn)行SG 濾波及分?jǐn)?shù)階微分處理,處理后的光譜采用CARS 算法進(jìn)行特征波段優(yōu)選,選出的波段組合用于偏最小二乘回歸PLSR(Partial Least Squares Regression)建立反演模型,最終將構(gòu)建的最優(yōu)模型應(yīng)用于航空高光譜圖像上反演重金屬含量,獲得重金屬含量分布圖。

    圖5 技術(shù)路線圖Fig.5 Flowchart of the proposed method

    3.1 樣本擴(kuò)充

    由于研究區(qū)的地形較平坦,地表一致性較好,可以認(rèn)為所采集土壤樣本的重金屬含量可以代表該樣本小范圍內(nèi)的土壤重金屬含量,又因?yàn)楸狙芯克@取的高光譜圖像空間分辨率為1 m,可認(rèn)為鄰近像元所對(duì)應(yīng)的土壤中重金屬含量差異極小。因此,本研究根據(jù)樣本經(jīng)緯度信息定位到圖像上的像元后,以該像元為中心像元,提取中心像元以及其八鄰域的像元光譜,依次計(jì)算八個(gè)鄰近像元光譜與中心像元光譜的歐氏距離,距離越近則認(rèn)為其與中心像元的相似度越高,選出歐氏距離最小的兩個(gè)像元,將3個(gè)像元(中心像元和兩個(gè)相似鄰近像元)的光譜共同作為該土壤樣本的反射率光譜信息,達(dá)到將土壤樣本的光譜數(shù)據(jù)擴(kuò)大3 倍的效果。本研究采集的72 個(gè)土壤樣本,用此方法進(jìn)行樣本擴(kuò)充后,獲得216條像元反射率光譜數(shù)據(jù)。

    3.2 分?jǐn)?shù)階微分(FOD)

    分?jǐn)?shù)階微分是數(shù)學(xué)中的重要分支,它將經(jīng)典的整數(shù)階微分推廣為任意階,能更敏感地捕捉光譜反射率細(xì)節(jié)的變化。到目前為止,還沒(méi)有一個(gè)統(tǒng)一的公式來(lái)定義分?jǐn)?shù)階微分。數(shù)學(xué)家們從不同的角度分析,得出了分?jǐn)?shù)階微積分的不同定義。目前,分?jǐn)?shù)階微分常見(jiàn)的表達(dá)形式主要包括Riemann-Liouville(R-L)、Grunwald-Letnikov(G-L)和Caputo。其中最常用的形式為G-L(Wang等,2017),本研究采用G-L形式進(jìn)行微分,其定義如下:

    式中,α為任意階數(shù);h為微分步長(zhǎng);t與α分別為微分的上、下限;Γ(α)為Gamma 函數(shù),Γ(α)滿足:

    令h=1,式(1)能夠推導(dǎo)出函數(shù)f(x)分?jǐn)?shù)階微分的表達(dá)式為

    3.3 競(jìng)爭(zhēng)自適應(yīng)重加權(quán)采樣法(CARS)

    CARS 是一種以“適者生存”為原則的篩選算法,基于與模型性能相關(guān)的統(tǒng)計(jì)數(shù)據(jù),通過(guò)逐步去除冗余和不重要的變量來(lái)選擇信息量大的變量,將其應(yīng)用于高光譜數(shù)據(jù)進(jìn)行篩選時(shí),其本質(zhì)是使用自適應(yīng)加權(quán)采樣技術(shù)來(lái)保留回歸中絕對(duì)系數(shù)較大的光譜波段(Tan等,2021)。

    圖6 為CARS 的方法流程圖,CARS 通過(guò)N次蒙特卡羅抽樣迭代地選擇N個(gè)波段子集,最終目的是選擇一個(gè)能使交叉驗(yàn)證的均方根誤差(RMSEcv)最低的最佳變量子集。在高光譜數(shù)據(jù)上使用CARS篩選的具體步驟如下:

    圖6 CARS算法流程圖Fig.6 Flowchart of CARS

    (1)使用蒙特卡羅抽樣方法從數(shù)據(jù)集中隨機(jī)選擇具有固定比例的樣本,然后用這些樣本建立PLSR 模型。PLSR 模型可用公式表達(dá)為y=bX+e,e為常數(shù)項(xiàng),b為回歸系數(shù)向量,b=[b1,b2,…,bn],b中第i個(gè)元素的絕對(duì)值|b|(1≤i≤n)表示第i個(gè)波段對(duì)y的貢獻(xiàn),波段對(duì)目標(biāo)變量的貢獻(xiàn)越大,該波段就越重要。為評(píng)價(jià)每個(gè)波長(zhǎng)的重要性,定義權(quán)重Wi:

    通過(guò)CARS 算法去掉的變量,其權(quán)重Wi均設(shè)為0;

    (2)利用指數(shù)遞減函數(shù)EDF(Exponentially Decreasing Function)和自適應(yīng)重加權(quán)采樣ARS(Adaptive Reweighted Sampling)分別強(qiáng)行和競(jìng)爭(zhēng)性地消除權(quán)重低的波段;

    (3)重復(fù)步驟(1)—(2),直到達(dá)到N次采樣運(yùn)行,最后選擇RMSEcv 最低的波段子集作為最優(yōu)波段組合。

    由于蒙特卡洛抽樣每次迭代抽取固定比例的樣本參與運(yùn)算,不使用所有樣本,因此選出的波段組合具有更好的適應(yīng)性;步驟(2)中使用指數(shù)遞減函數(shù)去除波段,能在迭代前期篩選掉大量重要性低的波段,使篩選過(guò)程的運(yùn)算量顯著下降。根據(jù)多次實(shí)驗(yàn)結(jié)果,將蒙特卡羅抽樣次數(shù)設(shè)置為50次,每次迭代抽取90%的樣本用于運(yùn)算。

    3.4 模型構(gòu)建與評(píng)價(jià)指標(biāo)

    PLSR 是土壤光譜分析中最為常用的反演模型,在處理共線性強(qiáng)、計(jì)算復(fù)雜度高的問(wèn)題上具有天然優(yōu)勢(shì)(Rossel 和Behrens,2010),能夠應(yīng)對(duì)自變量多于樣本個(gè)數(shù)的情況,因此其非常適合處理高光譜數(shù)據(jù),并且已被廣泛的應(yīng)用于土壤重金屬反演。本文采用CARS+PLSR 的模型構(gòu)建方法,先通過(guò)CARS篩選出特征波段,再將特征波段輸入PLSR 建立反演模型。在PLSR 建模中,通過(guò)留一交叉驗(yàn)證的最小均方根誤差確定PLSR 的潛變量個(gè)數(shù)。由于CARS的蒙特卡羅抽樣具有隨機(jī)性,為了避免實(shí)驗(yàn)的偶然性影響,本文對(duì)每組反演實(shí)驗(yàn)均進(jìn)行5 次,以5 次反演實(shí)驗(yàn)中的精度最優(yōu)值作為該模型精度將精度最優(yōu)的模型應(yīng)用于高光譜圖像,得到土壤重金屬的含量分布圖。

    精度評(píng)定采用預(yù)測(cè)均方根誤差RMSEP(Root Mean Square Error of Prediction)、相對(duì)分析誤差RPD(Ratio of Prediction to Deviation)和決定系數(shù)R2(Coefficient of Determination)3 個(gè)評(píng)價(jià)指標(biāo)。3 個(gè)評(píng)價(jià)指標(biāo),RMSEP 值越小,RPD 值越大,R2值越接近1,說(shuō)明反演模型的精度越高;反之,RMSEP值越大,RPD 值越小,R2值越小,說(shuō)明反演模型的精度越低,本文模型優(yōu)劣參考現(xiàn)有的土壤屬性含量高光譜反演的評(píng)價(jià)標(biāo)準(zhǔn)(Wang 等,2014):出色模型,R2>0.9;良好模型,0.9>R2>0.8;近似模型,0.8>R2>0.65;具有一定反演能力,0.65>R2>0.50;不具備反演能力,0.50>R2。

    4 實(shí)驗(yàn)結(jié)果與分析

    反演中將216個(gè)樣本按照2∶1的比例劃分為訓(xùn)練集和測(cè)試集。由于樣本擴(kuò)充出的3個(gè)樣本重金屬含量相同,如果將樣本重金屬含量排序后采用分層抽樣方法,按2∶1進(jìn)行劃分,會(huì)使訓(xùn)練集和測(cè)試集的樣本含量分布完全一樣,從而導(dǎo)致反演精度偏高。為了解決該問(wèn)題,選出具有代表性的訓(xùn)練集,采用Kennard-Stone(KS)算法(Zhang 等,2017)來(lái)劃分樣本集,具體的實(shí)施步驟如下:

    (1)首先計(jì)算兩兩樣本之間的光譜歐氏距離,選擇光譜距離最大的兩個(gè)樣本;

    (2)然后分別計(jì)算剩余樣本與已選兩樣本之間的光譜歐氏距離;

    (3)對(duì)于每個(gè)剩余樣本而言,計(jì)算其與已選各樣本之間的最短光譜距離,選擇這些最短光譜距離中相對(duì)最大的距離所對(duì)應(yīng)的樣本,作為新入選的樣本;

    (4)重復(fù)步驟(3),直至所選樣本的個(gè)數(shù)等于事先設(shè)定的數(shù)目為止。

    通過(guò)KS 算法選出144 個(gè)樣本作為訓(xùn)練集,剩下的72 個(gè)樣本作為測(cè)試集,樣本集劃分結(jié)果通過(guò)直方圖展示于圖7,3 種重金屬的訓(xùn)練集、測(cè)試集和總樣本集的含量直方圖分布一致且有一定差異性,表明該方法劃分出的訓(xùn)練集和測(cè)試集具有較好的代表性。

    圖7 各樣本集的重金屬含量直方圖Fig.7 Histogram of heavy metal content for each sample set

    4.1 分?jǐn)?shù)階微分的最佳階數(shù)確定

    由于分?jǐn)?shù)階微分需使用相鄰的波段進(jìn)行計(jì)算,理想條件下的土壤反射率光譜應(yīng)是一條連續(xù)的光滑曲線,相鄰波段的反射率不應(yīng)出現(xiàn)較大變化。本研究使用的航空高光譜圖像采集350—1000 nm和1000—2500 nm 波長(zhǎng)范圍的傳感器不同,導(dǎo)致光譜在1000 nm 波長(zhǎng)處的反射率曲線不夠平滑。此外,由于1350—1450 nm 的水汽吸收帶被剔除后,導(dǎo)致該光譜區(qū)間兩端波長(zhǎng)相差100 nm 的兩個(gè)波段變?yōu)橄噜彶ǘ?。為了確保微分時(shí)相鄰波段反射率的差異在合理范圍內(nèi),將光譜分為350—1000 nm、1000—1350 nm 和1450—1800 nm 這3 個(gè)區(qū)間來(lái)分別進(jìn)行分?jǐn)?shù)階微分處理。

    從圖8中可以看出,隨著階數(shù)的增加,光譜的特征變的越來(lái)越明顯,峰谷差異越發(fā)增大;0—1階的微分結(jié)果顯示出,階數(shù)增大的過(guò)程中,平緩區(qū)間的波段反射率值越發(fā)趨近于0,而波峰波谷區(qū)間的波段反射率值被逐漸放大;1—2 階的微分結(jié)果中,反射率曲線的波形已較為相似(由于1—2階內(nèi)的微分曲線較為相似,圖8中僅展示1階、1.5階和2 階的微分結(jié)果),均與原始光譜曲線差異較大,特征放大的效果明顯,隨著階數(shù)的增加,變化主要體現(xiàn)光譜曲線的極值不斷增大。

    圖8 不同階數(shù)微分后的像元光譜反射曲線Fig.8 Image spectral curves with different orders of differential

    在不同階數(shù)的微分處理后,反射率曲線中的光譜信息有明顯差異,為了選出適合于反演各重金屬的最佳微分階數(shù),在0—2階的范圍內(nèi),以0.1階為間隔,將不同階數(shù)微分后的光譜輸入CARS+PLSR構(gòu)建反演模型,以測(cè)試集精度為選擇依據(jù)。為了避免實(shí)驗(yàn)的偶然性影響,確保選出的階數(shù)能使反演精度達(dá)到最優(yōu),每組階數(shù)下均進(jìn)行了5次反演建模。

    圖9顯示了在不同階數(shù)微分條件下,CARS+PLSR反演3種重金屬的測(cè)試集R2最大值和平均值。對(duì)于Pb,階數(shù)為1.2 時(shí)測(cè)試集精度的最大值和平均值都達(dá)到了峰值(分別為0.7974 和0.7431);對(duì)于Zn,當(dāng)階數(shù)為0.8 時(shí),R2最大值曲線達(dá)到峰值0.8690,此時(shí)平均值曲線的值為0.8096,與平均值曲線的峰值0.8126 相差極小,因此將階數(shù)0.8 認(rèn)為是反演重金屬Zn 的最佳微分階數(shù);對(duì)于Ni,當(dāng)階數(shù)為0.3 時(shí)最大值曲線和平均值曲線均達(dá)到峰值(分別為0.8303 和0.7681)。同時(shí)考慮精度的最優(yōu)值和平均值,保證選出的最佳階數(shù)能實(shí)現(xiàn)反演精度最優(yōu)且模型穩(wěn)定性好。

    圖9 不同階數(shù)下反演3種重金屬的測(cè)試集精度Fig.9 Estimation accuracy of three heavy metals at different orders

    每種重金屬反演時(shí)光譜區(qū)間內(nèi)對(duì)不同種類重金屬敏感有效的波段不同,即各重金屬的特征波段不同,不同階數(shù)的微分處理后對(duì)各區(qū)間波段的特征化效果有所差異,因此各重金屬的最佳微分階數(shù)不盡相同。圖9 中當(dāng)階數(shù)為0 時(shí),代表不做微分處理使用原始像元光譜,當(dāng)階數(shù)為1 和2 時(shí)相當(dāng)于使用一階微分和二階微分,3種重金屬的反演結(jié)果都顯示出使用最佳階數(shù)分?jǐn)?shù)階微分的精度高于原始像元光譜、一階微分和二階微分光譜的反演精度,證明了使用分?jǐn)?shù)階微分能夠有效增強(qiáng)光譜特征,提高重金屬反演精度,并且分?jǐn)?shù)階微分比整數(shù)階微分能更敏感地突出對(duì)土壤重金屬反演有益的光譜信息。根據(jù)本節(jié)的實(shí)驗(yàn)結(jié)果,本文在反演重金屬Pb、Zn、Ni的各項(xiàng)實(shí)驗(yàn)中,分別選定1.2、0.8和0.3階的分?jǐn)?shù)階微分對(duì)像元光譜進(jìn)行處理。

    4.2 樣本擴(kuò)充有效性分析

    為驗(yàn)證樣本擴(kuò)充的有效性,每種重金屬均進(jìn)行兩組反演實(shí)驗(yàn),一組對(duì)樣本進(jìn)行擴(kuò)充,使用216個(gè)樣本參與建模,另一組不擴(kuò)充樣本,使用72 個(gè)樣本,反演每種重金屬時(shí),都以其最佳階數(shù)對(duì)像元光譜進(jìn)行微分處理后,建立CARS+PLSR模型。

    表1 中展示了各重金屬的兩組對(duì)照實(shí)驗(yàn)結(jié)果,在沒(méi)有進(jìn)行樣本擴(kuò)充前,小樣本問(wèn)題引起的過(guò)擬合現(xiàn)象明顯,3 種重金屬的反演模型訓(xùn)練集精度R2都大于0.98,然而只有Zn 的測(cè)試集精度為0.8178,Pb和Ni的測(cè)試集精度都較低,R2小于0.7;對(duì)樣本集進(jìn)行擴(kuò)充后建立的反演模型,過(guò)擬合現(xiàn)象得到了很好的緩解,3種重金屬的訓(xùn)練集精度和測(cè)試集精度的R2差距小于0.05。樣本擴(kuò)充后測(cè)試集的精度都得到了明顯的提升,Pb 的測(cè)試集R2從0.6128提升到0.7974,Zn 的精度從0.8178 提升到0.8690,Ni 的R2從0.6969 提升到0.8303。由此可見(jiàn),樣本擴(kuò)充不僅緩和了模型的過(guò)擬合現(xiàn)象,還有效地提升了3種重金屬的反演精度。

    表1 樣本擴(kuò)充前后的3種重金屬反演精度( 和 分別代表訓(xùn)練集和測(cè)試集的反演精度R2)Table 1 The estimation accuracies of three heavy metals before and after sample expansion( and represent the inversion accuracies R2 of the training and test sets,respectively)

    圖10 展示了樣本擴(kuò)充前后重金屬Ni 的反演結(jié)果散點(diǎn)圖,樣本中Ni含量大于2000 mg/kg的高含量樣本僅有兩個(gè),從圖10(a)(b)可以看出,兩個(gè)高含量樣本均被選入訓(xùn)練集,測(cè)試集中未包含高含量樣本,難以通過(guò)測(cè)試集精度評(píng)估該模型應(yīng)用于圖像后在高含量區(qū)域的反演能力。從圖10(c)(d)可以看出,樣本擴(kuò)充后訓(xùn)練集和測(cè)試集都包括了高含量的樣本,因此測(cè)試集精度可以代表在高含量區(qū)域的反演效果。在使用3.1 節(jié)所述的樣本擴(kuò)充條件下,使用KS算法劃分樣本集能夠確保高含量樣本存在于訓(xùn)練集和測(cè)試集,使本研究方法在反演含量偏分布的重金屬時(shí),能夠得到可靠性較高的反演模型。

    圖10 樣本擴(kuò)充前后重金屬Ni反演散點(diǎn)圖Fig.10 Scatter plot of Ni estimation results before and after sample expansion

    4.3 波段選擇算法的建模精度對(duì)比分析

    GA 在PLSR 建模中被認(rèn)為是一種有效的波段選擇算法(Leardi 和González,1998)。此外,相關(guān)系數(shù)法在土壤重金屬反演中被廣泛使用,因此本研究將GA+PLSR、CC+PLSR 兩種建模方法與CARS+PLSR 進(jìn)行對(duì)比分析。參考已有的GA+PLSR研究(Sun 等,2022),GA 的參數(shù)設(shè)置為染色體個(gè)數(shù)20,迭代次數(shù)150,代際間隙90%,基因變異概率10%。CC 的相關(guān)性篩選閾值通過(guò)在0—1 范圍內(nèi),以0.1 為間隔設(shè)置閾值進(jìn)行3 種重金屬的反演實(shí)驗(yàn),以CC+PLSR 的反演精度最優(yōu)為原則設(shè)置閾值的參數(shù)。

    通過(guò)分?jǐn)?shù)階微分對(duì)擴(kuò)充后的216條像元光譜進(jìn)行處理后,分別用CC、GA 和CARS 等3 種算法篩選波段,選出的3 組波段組合分別構(gòu)建PLSR 反演模型。

    表2 展示了3 種重金屬在使用不同波段選擇算法下的反演精度,表中的T(s)代表單次模型構(gòu)建所花費(fèi)的秒數(shù),分別代表訓(xùn)練集和測(cè)試集的反演精度R2。3種建模方法中,CC+PLSR 構(gòu)建的模型測(cè)試集精度最低,甚至在反演Pb 時(shí)R2<0.5,模型沒(méi)有估算能力,并且在3種重金屬反演中均出現(xiàn)過(guò)擬合現(xiàn)象,反演Pb 和Zn 時(shí)過(guò)擬合最為顯著,訓(xùn)練集精度遠(yuǎn)高于測(cè)試集精度。

    表2 不同波段選擇算法的重金屬反演精度Table 2 Heavy metal estimation accuracy of different band selection algorithms

    GA+PLSR 的模型精度較高,3 種重金屬反演結(jié)果僅略低于CARS+PLSR 模型,反演Zn時(shí)的測(cè)試集精度R2為0.8119,屬于良好模型,反演Pb 和Ni時(shí)測(cè)試集精度R2也均大于0.78,模型沒(méi)有顯示出過(guò)擬合現(xiàn)象。但由于GA 在篩選波段時(shí),需要多次迭代選擇最優(yōu)波段組合,每次迭代中所有波段均參與運(yùn)算,并且由于GA 算法中有種群機(jī)制的設(shè)置,在本研究的參數(shù)設(shè)置下,每次迭代需生成20種波段組合,意味著有20倍的總波段數(shù)參與運(yùn)算,導(dǎo)致整個(gè)模型的構(gòu)建時(shí)間較長(zhǎng),3種重金屬的GA+PLSR反演模型構(gòu)建時(shí)間都在300 s 以上,耗費(fèi)時(shí)間遠(yuǎn)遠(yuǎn)多于CARS和CC兩種篩選方法。

    CARS+PLSR 在3 種建模方法中獲得了最高的反演精度,Zn 和Ni 的反演精度R2都高于0.8,Pb的反演精度R2也達(dá)到0.7974,并且模型沒(méi)有呈現(xiàn)出明顯的過(guò)擬合現(xiàn)象,訓(xùn)練集和測(cè)試集都保持著較高的精度。此外,雖然CARS 和GA 算法都需要多次迭代尋找最優(yōu)波段組合,但由于CARS 內(nèi)的EDF 和ARS 算法在迭代初期就快速、強(qiáng)力的去除重要性低的波段,因此迭代中后期僅有少數(shù)波段參與運(yùn)算,所需的建模時(shí)間顯著縮短。因此CARS是3種方法中最優(yōu)的波段選擇算法,能夠更快篩選出最優(yōu)波段組合,使模型具有更好的反演能力。

    4.4 CARS的波段選擇結(jié)果分析

    圖11 展示了CARS 選出的最優(yōu)波段組合的直方圖,直方圖內(nèi)柱形圖越高代表該波長(zhǎng)范圍內(nèi)選擇的波段越多,3 種重金屬的CARS 篩選結(jié)果在不同波段區(qū)間選擇的波段數(shù)量不同,說(shuō)明各波長(zhǎng)范圍對(duì)不同重金屬反演的貢獻(xiàn)不同。土壤反射光譜分析表明土壤在VNIR-SWIR 區(qū)間的吸收特征主要由土壤有機(jī)質(zhì)、鐵氧化物等土壤光譜活性物質(zhì)引起(Kooistra 等,2003)。由于土壤光譜活性物質(zhì)在土壤中對(duì)重金屬的吸附具有主導(dǎo)作用,可據(jù)此推算土壤重金屬濃度,這是通過(guò)反射光譜間接反演重金屬的主要機(jī)理(Rathod 等,2013)。然而有機(jī)質(zhì)、鐵氧化物等土壤光譜活性物質(zhì)對(duì)于不同的重金屬吸附強(qiáng)度不同,因此反演不同重金屬時(shí)起主導(dǎo)作用的土壤組分不同(Covelo 等,2007),現(xiàn)有研究表明土壤鐵氧化物對(duì)重金屬Pb 的吸附作用是土壤反射光譜反演Pb 含量的主要機(jī)理,土壤有機(jī)質(zhì)對(duì)重金屬Ni、Zn 的吸附作用是土壤反射光譜反演Ni、Zn 的主要機(jī)理,鐵氧化物的吸收特征在500 nm(Wu 等,2007),而600—800 nm 附近的吸收峰被認(rèn)為是土壤有機(jī)質(zhì)的吸收特征(徐彬彬 等,1991)。已有研究證明500 nm 的鐵氧化物特征波段對(duì)反演重金屬Pb 是有效的(張霞 等,2022),Sun 等(Sun 和Zhang,2017;Sun 等,2018)證明了使用600—800 nm的有機(jī)質(zhì)特征譜段對(duì)反演Zn和Ni是有效的。根據(jù)圖11所示,在反演Pb時(shí),CARS篩選結(jié)果的直方圖中500 nm 范圍內(nèi)的柱形圖最高;在反演Zn 和Ni 時(shí),直方圖內(nèi)最高的柱形圖落在600—800 nm區(qū)間和1600 nm范圍內(nèi),Li等(2022)指出有機(jī)質(zhì)的官能團(tuán)對(duì)1600 nm 光譜影響顯著,證明1600 nm 同樣是對(duì)有機(jī)質(zhì)敏感的特征波段。綜上所述CARS的波段選擇結(jié)果與已有反演機(jī)理研究保持一致,證明CARS能夠篩選出對(duì)反演重金屬有益且合理的特征波段。

    圖11 3種重金屬反演模型的CARS優(yōu)選波段結(jié)果Fig.11 Results of CARS band selection for three heavy metal estimation models

    4.5 重金屬含量分布圖分析

    為了保證模型應(yīng)用到高光譜圖像上時(shí)反演結(jié)果的可靠性,通過(guò)光譜角匹配法(Ramirez-Lopez 等,2013)篩選出圖像上與土壤樣本像元光譜相似度高的區(qū)域,將3種重金屬的最優(yōu)模型分別應(yīng)用于該區(qū)域進(jìn)行重金屬含量反演,圖12(b)(c)(d)展示了依據(jù)此方法制作出的3 種重金屬含量分布圖。由于Ni 的含量區(qū)間跨度較大,連續(xù)的圖例難以體現(xiàn)出含量差距,圖12(d)的Ni含量分布圖中采用了分級(jí)圖例進(jìn)行展示,根據(jù)2018 年發(fā)布的《土壤環(huán)境質(zhì)量 建設(shè)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)》(GB 36600-2018)相關(guān)規(guī)定,礦區(qū)作為二類用地,Ni 含量的篩選值(可能存在風(fēng)險(xiǎn))為900 mg/kg,因此設(shè)置900—2500 mg/kg 的級(jí)別代表超過(guò)篩選值,用于表示對(duì)人體健康可能存在風(fēng)險(xiǎn)的區(qū)域,900 mg/kg 以下采用分位數(shù)法進(jìn)行分級(jí),確保每一級(jí)別內(nèi)所含的像元個(gè)數(shù)大致相等。

    圖12 重金屬含量分布圖Fig.12 Heavy metal content distribution map

    研究區(qū)屬于鎳銅礦山區(qū)域,可以看出Ni 在整個(gè)研究區(qū)范圍內(nèi)的含量偏高,但幾乎沒(méi)有超過(guò)篩選值的情況,大部分區(qū)域的Ni 含量明顯高于Pb 和Zn 的含量值,彭再華和蔣素芳(2018)研究也表明在哈密黃山南礦山中最主要的有價(jià)元素為銅和鎳,其他金屬元素含量較低。重金屬Pb 和Zn 具有相似的分布趨勢(shì),李玲等(2020)在新疆礦區(qū)的研究也指出Pb、Zn 含量表現(xiàn)出顯著相關(guān)性。3 種重金屬分布的共同點(diǎn)為:在北部以及南部區(qū)域的重金屬含量較高,在靠近居住區(qū)的東部區(qū)域含量較低。

    圖13展示了研究區(qū)的高程分布,可以看出研究區(qū)的地勢(shì)呈現(xiàn)出東北高,西南低。經(jīng)實(shí)地調(diào)研得知圖12(a)中的A 區(qū)域用于堆放采礦渣,圖12(a)中的B區(qū)域是采礦區(qū),采礦被認(rèn)為是重金屬污染的最重要因素之一,重金屬會(huì)伴隨尾礦渣,廢水等進(jìn)入土壤(王海洋 等,2022),該采礦區(qū)于2007年開(kāi)始勘探,已在此地進(jìn)行了多年的采礦活動(dòng),由于研究區(qū)西南地勢(shì)較低,重金屬在重力作用下不斷發(fā)生遷移,長(zhǎng)時(shí)間作用下導(dǎo)致研究區(qū)南部、西南部區(qū)域受到污染,因此3 種重金屬在南部的含量較高。

    圖13 研究區(qū)高程圖Fig.13 Elevation map of the study area

    研究區(qū)北部含量較高的區(qū)域?qū)?yīng)圖13 中高程出現(xiàn)顯著變化的部分,圖12(a)中可看出該區(qū)域圖像與周圍有明顯差異,根據(jù)實(shí)地調(diào)研得知該區(qū)域?yàn)樯襟w斜坡,坡上的土壤粒徑極小呈細(xì)砂狀,坡面平滑。在坡面上采集的4個(gè)土壤樣本重金屬含量均不低,其中3 個(gè)樣本的Pb、Zn 含量高于平均值,2 個(gè)樣本的Ni含量遠(yuǎn)高于平均值,且其中1 個(gè)達(dá)1220 mg/kg,表明該區(qū)域土壤的重金屬含量偏高,反演結(jié)果與采樣分析結(jié)果一致。

    5 結(jié)論

    目前針對(duì)高光譜圖像的重金屬反演研究較少,本文以3 種重金屬Pb、Zn、Ni 為例,開(kāi)展面向高光譜圖像的反演研究。由于土壤樣本點(diǎn)與圖像土壤像元數(shù)差距懸殊,通過(guò)提取采樣點(diǎn)鄰近像元的方法擴(kuò)充樣本數(shù)據(jù),同時(shí)增加樣本的光譜多樣性;采用FOD 增強(qiáng)對(duì)重金屬反演有益的光譜信息,提高反演精度;由于高光譜數(shù)據(jù)的波段眾多,使用CARS 算法選出特征波段以建立反演模型,從降維角度減弱小樣本的影響。以新疆哈密黃山南銅鎳礦區(qū)為研究區(qū),分析了該反演方法的有效性。

    研究結(jié)果表明,樣本擴(kuò)充有效地提升了3種重金屬的測(cè)試集反演精度,緩解了模型的過(guò)擬合問(wèn)題,擴(kuò)充前3種重金屬的訓(xùn)練集精度R2均在0.98以上,遠(yuǎn)高于測(cè)試集精度,模型明顯過(guò)擬合,而樣本擴(kuò)充后,3 種重金屬的測(cè)試集反演精度都得到了提升,同時(shí),訓(xùn)練集R2也在0.8 以上,與測(cè)試集精度相近,過(guò)擬合問(wèn)題得到顯著改善。

    3 種重金屬在最佳階數(shù)分?jǐn)?shù)階微分下的模型反演精度均高于使用整數(shù)階微分的反演精度,說(shuō)明分?jǐn)?shù)階微分能夠更加有效地突出對(duì)土壤重金屬反演有益的光譜信息。在最佳階數(shù)微分的基礎(chǔ)上,將CARS與常用的CC和GA兩種波段選擇算法進(jìn)行了對(duì)比,CARS+PLSR 模型的精度最高,并且建模耗時(shí)短,本研究認(rèn)為CARS 是3 種算法中最優(yōu)的波段選擇策略。

    對(duì)CARS的波段選擇結(jié)果進(jìn)行分析,反演Pb時(shí),選出的波段大多位于500 nm 范圍內(nèi),該區(qū)間是對(duì)Pb 敏感有效的鐵氧化物的特征波段,反演Zn、Ni時(shí),選出的波段大多位于600—800 nm 和1600 nm附近,該波長(zhǎng)范圍屬于對(duì)Zn、Ni 敏感的有機(jī)質(zhì)特征譜段,與現(xiàn)有反演機(jī)理研究相一致,證明CARS算法能夠有效地篩選出光譜中對(duì)反演重金屬有益的特征波段。

    本研究方法應(yīng)用于航空高光譜圖像上反演土壤重金屬Pb、Zn、Ni 的含量,反演精度較高,并且反演出的3種重金屬分布圖與實(shí)際相符,表明該方法有很好的魯棒性,具有反演多種土壤重金屬的能力。但該方法尚需要在其它礦區(qū)、農(nóng)作區(qū)驗(yàn)證其適用性,且土壤類型的影響也需要進(jìn)一步探討。

    猜你喜歡
    模型
    一半模型
    一種去中心化的域名服務(wù)本地化模型
    適用于BDS-3 PPP的隨機(jī)模型
    提煉模型 突破難點(diǎn)
    函數(shù)模型及應(yīng)用
    p150Glued在帕金森病模型中的表達(dá)及分布
    函數(shù)模型及應(yīng)用
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    3D打印中的模型分割與打包
    国产精品一区二区三区四区久久| 久久久久精品国产欧美久久久| 日日摸夜夜添夜夜添小说| 国产精品98久久久久久宅男小说| avwww免费| 欧美性猛交╳xxx乱大交人| 禁无遮挡网站| 观看免费一级毛片| 99视频精品全部免费 在线| 亚洲欧美日韩无卡精品| 少妇裸体淫交视频免费看高清| 亚洲欧美日韩高清专用| 日本在线视频免费播放| 日日摸夜夜添夜夜添小说| 国产精品一区www在线观看 | 欧美一区二区精品小视频在线| 亚洲国产精品合色在线| 美女高潮的动态| 久久精品综合一区二区三区| 黄色女人牲交| 一夜夜www| 久久久久性生活片| 中文字幕高清在线视频| 日日啪夜夜撸| 露出奶头的视频| 欧美色视频一区免费| 国产综合懂色| 波野结衣二区三区在线| 他把我摸到了高潮在线观看| 99riav亚洲国产免费| 国产精品久久电影中文字幕| av天堂中文字幕网| 男女做爰动态图高潮gif福利片| 91在线精品国自产拍蜜月| 麻豆av噜噜一区二区三区| 亚洲精华国产精华液的使用体验 | 99久久九九国产精品国产免费| 精品人妻视频免费看| 我要看日韩黄色一级片| 亚洲欧美激情综合另类| 国产伦一二天堂av在线观看| 97人妻精品一区二区三区麻豆| 听说在线观看完整版免费高清| 国产美女午夜福利| 久久久久久久久中文| 亚洲自拍偷在线| 91久久精品国产一区二区三区| 日韩 亚洲 欧美在线| 久久久久久久亚洲中文字幕| av在线蜜桃| 88av欧美| 亚洲成人久久爱视频| 欧美最黄视频在线播放免费| 亚洲精品亚洲一区二区| 久久久久久久午夜电影| 热99re8久久精品国产| 国产欧美日韩精品一区二区| 此物有八面人人有两片| 亚洲av一区综合| 免费看av在线观看网站| 黄色配什么色好看| 男女做爰动态图高潮gif福利片| 嫩草影院入口| 九色国产91popny在线| 亚洲人成网站在线播放欧美日韩| 亚洲国产欧洲综合997久久,| 欧美色视频一区免费| 亚洲三级黄色毛片| 国产精品久久久久久亚洲av鲁大| 18禁在线播放成人免费| 欧美另类亚洲清纯唯美| 九九在线视频观看精品| 国产成人福利小说| 熟女人妻精品中文字幕| 91狼人影院| 亚洲精品一区av在线观看| 嫩草影院精品99| aaaaa片日本免费| av在线亚洲专区| 成人三级黄色视频| 91久久精品国产一区二区成人| 亚洲av不卡在线观看| 99久久九九国产精品国产免费| 搡女人真爽免费视频火全软件 | 一本一本综合久久| 国产精品野战在线观看| 国产精品无大码| 最近中文字幕高清免费大全6 | 波多野结衣高清无吗| 国产精品自产拍在线观看55亚洲| 中出人妻视频一区二区| 99久久九九国产精品国产免费| 国产伦精品一区二区三区视频9| 97热精品久久久久久| 老女人水多毛片| 欧美三级亚洲精品| 亚洲精品国产成人久久av| 少妇人妻精品综合一区二区 | 久久精品国产亚洲网站| 特级一级黄色大片| 色综合亚洲欧美另类图片| a在线观看视频网站| 久久国产精品人妻蜜桃| 男女啪啪激烈高潮av片| 日本爱情动作片www.在线观看 | 免费无遮挡裸体视频| 国产蜜桃级精品一区二区三区| 国模一区二区三区四区视频| 可以在线观看毛片的网站| 欧美又色又爽又黄视频| 免费人成视频x8x8入口观看| 最近最新中文字幕大全电影3| 国产午夜福利久久久久久| 亚洲熟妇熟女久久| 美女高潮的动态| 日韩精品中文字幕看吧| 99热只有精品国产| 欧美xxxx性猛交bbbb| 国产亚洲精品av在线| 18禁黄网站禁片免费观看直播| 色播亚洲综合网| 日韩欧美国产在线观看| 一区二区三区高清视频在线| 我要搜黄色片| 精品久久久久久成人av| 不卡一级毛片| 在线观看美女被高潮喷水网站| 亚洲欧美精品综合久久99| 伦理电影大哥的女人| 高清日韩中文字幕在线| 精品人妻熟女av久视频| 中文资源天堂在线| www.色视频.com| 男人和女人高潮做爰伦理| 免费一级毛片在线播放高清视频| 国产真实伦视频高清在线观看 | 欧美一区二区精品小视频在线| 乱人视频在线观看| 久久精品国产亚洲av天美| 日本黄色视频三级网站网址| 淫秽高清视频在线观看| 久久精品国产亚洲av天美| 日本黄大片高清| 亚洲无线在线观看| 色尼玛亚洲综合影院| 亚洲一区二区三区色噜噜| 91久久精品国产一区二区三区| 国产单亲对白刺激| 国产极品精品免费视频能看的| 长腿黑丝高跟| 免费看日本二区| 国产精品亚洲一级av第二区| 午夜福利欧美成人| 黄色视频,在线免费观看| 天美传媒精品一区二区| 男插女下体视频免费在线播放| 一卡2卡三卡四卡精品乱码亚洲| 看黄色毛片网站| 亚洲三级黄色毛片| 人妻制服诱惑在线中文字幕| 99国产极品粉嫩在线观看| 在线免费观看的www视频| 97超视频在线观看视频| 国产美女午夜福利| 免费无遮挡裸体视频| 精品午夜福利在线看| 美女 人体艺术 gogo| 欧美激情久久久久久爽电影| 欧美精品国产亚洲| 性欧美人与动物交配| 国产高清不卡午夜福利| 露出奶头的视频| 亚洲av电影不卡..在线观看| 亚洲成人久久性| 男女之事视频高清在线观看| 亚洲在线观看片| 日本五十路高清| 免费搜索国产男女视频| 国内精品久久久久精免费| 欧美中文日本在线观看视频| 中文亚洲av片在线观看爽| 免费不卡的大黄色大毛片视频在线观看 | 天堂√8在线中文| 桃色一区二区三区在线观看| av.在线天堂| 男女边吃奶边做爰视频| 精品99又大又爽又粗少妇毛片 | 干丝袜人妻中文字幕| 麻豆成人av在线观看| 亚洲av日韩精品久久久久久密| 亚洲欧美精品综合久久99| 国产午夜福利久久久久久| 最近视频中文字幕2019在线8| 欧美国产日韩亚洲一区| 村上凉子中文字幕在线| 国产欧美日韩一区二区精品| 欧美+日韩+精品| 他把我摸到了高潮在线观看| 国产视频内射| 国产午夜福利久久久久久| 国内毛片毛片毛片毛片毛片| 99精品在免费线老司机午夜| 亚洲精品成人久久久久久| 久久九九热精品免费| 欧美日本视频| 亚州av有码| 97超级碰碰碰精品色视频在线观看| 国产黄片美女视频| 欧美性猛交╳xxx乱大交人| 午夜久久久久精精品| 97碰自拍视频| 国内精品宾馆在线| 黄色配什么色好看| 精品久久久久久久久av| 如何舔出高潮| 日韩人妻高清精品专区| 在线观看午夜福利视频| 亚洲欧美日韩东京热| 亚洲国产欧洲综合997久久,| 亚洲第一电影网av| 精品人妻视频免费看| 午夜精品在线福利| 婷婷亚洲欧美| 超碰av人人做人人爽久久| 国产aⅴ精品一区二区三区波| av天堂在线播放| 久久午夜福利片| 小蜜桃在线观看免费完整版高清| 亚洲av美国av| 九九爱精品视频在线观看| 久久久久久九九精品二区国产| 欧美zozozo另类| 国产精品电影一区二区三区| 草草在线视频免费看| 中文字幕精品亚洲无线码一区| 国产精品98久久久久久宅男小说| 嫩草影院新地址| 亚洲三级黄色毛片| 久久中文看片网| 干丝袜人妻中文字幕| 国产免费一级a男人的天堂| 一区二区三区四区激情视频 | 真实男女啪啪啪动态图| 五月伊人婷婷丁香| aaaaa片日本免费| 成年女人看的毛片在线观看| 欧美成人a在线观看| 91在线观看av| 色哟哟·www| 国产精品永久免费网站| 欧美日韩精品成人综合77777| 国产成人aa在线观看| 欧美最新免费一区二区三区| 校园人妻丝袜中文字幕| 色综合婷婷激情| 香蕉av资源在线| 最新在线观看一区二区三区| 嫩草影视91久久| 最后的刺客免费高清国语| 一个人看视频在线观看www免费| 床上黄色一级片| 国产精品永久免费网站| 午夜福利欧美成人| 女人被狂操c到高潮| 老熟妇仑乱视频hdxx| 欧美成人a在线观看| 国产午夜精品论理片| 久久精品影院6| 国产精品久久电影中文字幕| 一级a爱片免费观看的视频| 人妻少妇偷人精品九色| 麻豆久久精品国产亚洲av| 久久久久久久精品吃奶| 搞女人的毛片| 18禁在线播放成人免费| 99精品在免费线老司机午夜| 日本精品一区二区三区蜜桃| 黄色丝袜av网址大全| 成人av一区二区三区在线看| 日日撸夜夜添| 国产极品精品免费视频能看的| 国产真实伦视频高清在线观看 | 久99久视频精品免费| 99精品在免费线老司机午夜| 日韩一本色道免费dvd| 他把我摸到了高潮在线观看| 久久久久久久久久久丰满 | 亚洲中文日韩欧美视频| 亚洲精品亚洲一区二区| 噜噜噜噜噜久久久久久91| 日韩 亚洲 欧美在线| 啦啦啦韩国在线观看视频| 身体一侧抽搐| 深夜精品福利| 麻豆成人午夜福利视频| 欧美另类亚洲清纯唯美| 天堂av国产一区二区熟女人妻| 欧美+亚洲+日韩+国产| 男人舔奶头视频| 五月玫瑰六月丁香| 午夜福利在线观看吧| 乱码一卡2卡4卡精品| 日本-黄色视频高清免费观看| 国产高清视频在线观看网站| 国产精品伦人一区二区| 国产乱人伦免费视频| 久久精品影院6| 免费看美女性在线毛片视频| 哪里可以看免费的av片| 天天躁日日操中文字幕| 成人鲁丝片一二三区免费| 黄色配什么色好看| 色综合亚洲欧美另类图片| 欧美中文日本在线观看视频| 国产高清视频在线播放一区| 午夜福利在线观看吧| 又爽又黄a免费视频| 亚洲精品亚洲一区二区| or卡值多少钱| 麻豆精品久久久久久蜜桃| 亚洲精华国产精华液的使用体验 | 伦理电影大哥的女人| 他把我摸到了高潮在线观看| 国产精华一区二区三区| 久久久久久久久中文| 无遮挡黄片免费观看| 男女视频在线观看网站免费| 天堂影院成人在线观看| 日韩中文字幕欧美一区二区| a级毛片免费高清观看在线播放| 夜夜夜夜夜久久久久| 欧美不卡视频在线免费观看| 亚洲精品色激情综合| 亚洲四区av| 中国美白少妇内射xxxbb| 亚洲精品色激情综合| 色尼玛亚洲综合影院| 男女下面进入的视频免费午夜| 男人狂女人下面高潮的视频| av国产免费在线观看| 干丝袜人妻中文字幕| 婷婷精品国产亚洲av在线| 18禁黄网站禁片午夜丰满| 精品一区二区三区人妻视频| 久久久色成人| 亚洲第一区二区三区不卡| 性色avwww在线观看| 黄色日韩在线| 午夜福利在线观看吧| 国语自产精品视频在线第100页| а√天堂www在线а√下载| 亚洲成av人片在线播放无| 色综合站精品国产| 日韩欧美在线乱码| 免费大片18禁| 搞女人的毛片| 99国产精品一区二区蜜桃av| 久久久久国产精品人妻aⅴ院| 日本a在线网址| 熟妇人妻久久中文字幕3abv| 久久6这里有精品| 两性午夜刺激爽爽歪歪视频在线观看| 观看美女的网站| 尾随美女入室| а√天堂www在线а√下载| 欧美绝顶高潮抽搐喷水| 岛国在线免费视频观看| 色综合色国产| 国模一区二区三区四区视频| 欧美三级亚洲精品| 精品一区二区三区视频在线| 91午夜精品亚洲一区二区三区 | 97人妻精品一区二区三区麻豆| 欧美3d第一页| 黄片wwwwww| 色av中文字幕| 久久精品国产亚洲av香蕉五月| 亚洲av电影不卡..在线观看| www.色视频.com| 级片在线观看| 天堂动漫精品| 麻豆一二三区av精品| 久久精品综合一区二区三区| 亚洲国产日韩欧美精品在线观看| 亚洲性夜色夜夜综合| 中国美白少妇内射xxxbb| 美女xxoo啪啪120秒动态图| 欧美日韩国产亚洲二区| 亚洲欧美日韩卡通动漫| 一个人看的www免费观看视频| 色播亚洲综合网| 波野结衣二区三区在线| 日本一二三区视频观看| av在线蜜桃| 久99久视频精品免费| 又粗又爽又猛毛片免费看| 亚洲人与动物交配视频| 亚洲美女黄片视频| 亚洲av电影不卡..在线观看| 日日撸夜夜添| 99精品在免费线老司机午夜| 乱人视频在线观看| 欧美激情久久久久久爽电影| 午夜老司机福利剧场| 亚洲中文字幕日韩| 别揉我奶头~嗯~啊~动态视频| 国产乱人视频| 免费看av在线观看网站| 久久精品夜夜夜夜夜久久蜜豆| 内地一区二区视频在线| 高清日韩中文字幕在线| 成人一区二区视频在线观看| 欧美性猛交╳xxx乱大交人| 国产精品99久久久久久久久| 欧美在线一区亚洲| 日本黄色视频三级网站网址| 精品人妻熟女av久视频| 亚洲成人精品中文字幕电影| 亚洲人成网站在线播| 久久这里只有精品中国| 亚洲av成人精品一区久久| 九色成人免费人妻av| 久久久久久久久中文| 国内毛片毛片毛片毛片毛片| 看十八女毛片水多多多| 亚洲四区av| 亚洲专区中文字幕在线| 又爽又黄无遮挡网站| 日韩精品有码人妻一区| 国产又黄又爽又无遮挡在线| 亚洲国产欧美人成| 国产成人一区二区在线| 亚洲五月天丁香| 久久香蕉精品热| 麻豆国产av国片精品| 日韩精品青青久久久久久| 女人被狂操c到高潮| 久久精品久久久久久噜噜老黄 | 亚洲精品粉嫩美女一区| 亚洲,欧美,日韩| 成人午夜高清在线视频| 在线观看免费视频日本深夜| 欧美色视频一区免费| 波多野结衣高清无吗| 免费看日本二区| 久久婷婷人人爽人人干人人爱| 男女那种视频在线观看| 欧美性猛交黑人性爽| 国产免费男女视频| 日韩欧美精品免费久久| 中文字幕av成人在线电影| 又爽又黄无遮挡网站| 男女做爰动态图高潮gif福利片| 国产精品一区二区免费欧美| a级毛片免费高清观看在线播放| 欧美激情国产日韩精品一区| 久久亚洲真实| 永久网站在线| 两人在一起打扑克的视频| 婷婷色综合大香蕉| 国产综合懂色| 露出奶头的视频| 黄色配什么色好看| 舔av片在线| 国内精品久久久久精免费| 国产精品爽爽va在线观看网站| 男女那种视频在线观看| 色吧在线观看| 久久久久久久亚洲中文字幕| 中文字幕熟女人妻在线| 久久精品国产亚洲av涩爱 | 黄片wwwwww| 国产大屁股一区二区在线视频| 国产精品一区二区三区四区久久| 久久久久免费精品人妻一区二区| 校园人妻丝袜中文字幕| 亚洲成人免费电影在线观看| 亚洲人与动物交配视频| 少妇丰满av| 色尼玛亚洲综合影院| 搡女人真爽免费视频火全软件 | 色吧在线观看| 午夜福利成人在线免费观看| 大型黄色视频在线免费观看| 最近中文字幕高清免费大全6 | 我的老师免费观看完整版| 欧美在线一区亚洲| 少妇熟女aⅴ在线视频| 深爱激情五月婷婷| 亚洲精华国产精华液的使用体验 | 日韩精品中文字幕看吧| 偷拍熟女少妇极品色| 亚洲中文日韩欧美视频| 欧美性猛交黑人性爽| 欧美日韩国产亚洲二区| 成人午夜高清在线视频| 日本免费a在线| 成人欧美大片| 免费高清视频大片| 国产大屁股一区二区在线视频| 热99在线观看视频| 麻豆av噜噜一区二区三区| 久久久午夜欧美精品| 国产美女午夜福利| 日韩人妻高清精品专区| 午夜精品久久久久久毛片777| 亚洲第一区二区三区不卡| 韩国av在线不卡| 日韩精品有码人妻一区| 成人欧美大片| 欧美在线一区亚洲| 日韩欧美 国产精品| 成人鲁丝片一二三区免费| 国产私拍福利视频在线观看| 在线播放无遮挡| 中国美女看黄片| 12—13女人毛片做爰片一| 成年免费大片在线观看| 99热网站在线观看| 在线播放无遮挡| 午夜福利在线观看吧| 又粗又爽又猛毛片免费看| 国产精品久久电影中文字幕| 亚洲国产欧美人成| 网址你懂的国产日韩在线| 99在线视频只有这里精品首页| 日本色播在线视频| 亚洲av成人精品一区久久| 日韩欧美免费精品| 日本-黄色视频高清免费观看| 啦啦啦啦在线视频资源| 深夜精品福利| 人妻丰满熟妇av一区二区三区| 亚洲内射少妇av| 中文字幕免费在线视频6| 黄色丝袜av网址大全| 欧美+亚洲+日韩+国产| 51国产日韩欧美| 成人特级黄色片久久久久久久| 乱人视频在线观看| 赤兔流量卡办理| 国产成年人精品一区二区| 精品久久久久久成人av| 亚洲欧美日韩无卡精品| 亚洲精品456在线播放app | 69人妻影院| 免费av观看视频| 国产精品三级大全| 国产av麻豆久久久久久久| 男女边吃奶边做爰视频| 国产色爽女视频免费观看| 免费不卡的大黄色大毛片视频在线观看 | 久久这里只有精品中国| 国产久久久一区二区三区| 久久久成人免费电影| 内地一区二区视频在线| 99久久精品热视频| 欧美中文日本在线观看视频| 国产 一区精品| 成熟少妇高潮喷水视频| 日韩人妻高清精品专区| 成年免费大片在线观看| 亚洲自拍偷在线| 啦啦啦啦在线视频资源| 一a级毛片在线观看| a级一级毛片免费在线观看| 中亚洲国语对白在线视频| 午夜免费成人在线视频| 国产在线精品亚洲第一网站| 91在线精品国自产拍蜜月| 久久久久久九九精品二区国产| 亚洲国产精品久久男人天堂| 我要看日韩黄色一级片| 亚洲第一区二区三区不卡| 亚洲欧美精品综合久久99| 最新中文字幕久久久久| 美女被艹到高潮喷水动态| av视频在线观看入口| 日韩一本色道免费dvd| av黄色大香蕉| 国产精品女同一区二区软件 | 女人被狂操c到高潮| АⅤ资源中文在线天堂| 婷婷精品国产亚洲av| 97超级碰碰碰精品色视频在线观看| 亚洲av免费在线观看| 两性午夜刺激爽爽歪歪视频在线观看| 99精品久久久久人妻精品| 国产主播在线观看一区二区| 午夜免费男女啪啪视频观看 | 淫妇啪啪啪对白视频| 久久久久久久精品吃奶| 国产精品久久久久久亚洲av鲁大| 国产色爽女视频免费观看| 日日干狠狠操夜夜爽| 亚洲精品亚洲一区二区| 真实男女啪啪啪动态图| 久久6这里有精品| 国产精品久久久久久亚洲av鲁大| 国产高清三级在线| 69av精品久久久久久| 男女啪啪激烈高潮av片| 九色国产91popny在线| 成人午夜高清在线视频| 久久欧美精品欧美久久欧美| 又黄又爽又刺激的免费视频.| 久久久久久久精品吃奶| 午夜免费成人在线视频| 国内精品久久久久精免费| 亚洲精华国产精华精| 无遮挡黄片免费观看| 免费一级毛片在线播放高清视频| 18禁裸乳无遮挡免费网站照片| 亚洲成人中文字幕在线播放| 久久久久久久久中文|