国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

基于RF-VR的紫丁香葉片葉綠素含量高光譜反演

2021-11-27 13:29肖志云王伊凝
關(guān)鍵詞:紫丁香波段葉綠素

肖志云, 王伊凝

(內(nèi)蒙古工業(yè)大學(xué) 電力學(xué)院,內(nèi)蒙古機(jī)電控制重點(diǎn)實(shí)驗(yàn)室,內(nèi)蒙古 呼和浩特 010051)

綠色植物的生長過程離不開光合作用,葉片葉綠素含量及其動態(tài)變化與光合作用能力密不可分,檢測葉片葉綠素含量,對植物的長勢監(jiān)測和精準(zhǔn)農(nóng)業(yè)的實(shí)施具有重要意義[1-3]。目前葉綠素含量測定方法主要分兩種:化學(xué)測定法[4]需要破壞植物樣本,耗時(shí)且費(fèi)力;SPAD(soil plant analysis development)葉綠素測定儀通過測量葉片對兩個(gè)波段(紅波段和近紅外波段)的吸收率,計(jì)算當(dāng)前葉片中葉綠素的相對含量,但其只可實(shí)時(shí)測出局部點(diǎn)光譜對應(yīng)的SPAD值,無法反映整張葉片各像素點(diǎn)的葉綠素分布差異,對儀器精度有很高的依賴性。葉綠素含量的變化會引起植物反射光譜特征的變化。高光譜技術(shù)[5]具有圖譜信息合一的優(yōu)勢,既可以利用多波段光譜對葉片葉綠素含量進(jìn)行定量反演,又可以利用圖像像素點(diǎn)分布進(jìn)行葉綠素分布可視化研究,這就為利用高光譜技術(shù)獲取植物生化參數(shù)提供了理論基礎(chǔ),但由于高光譜所含波段數(shù)量大,波段間相關(guān)性強(qiáng)導(dǎo)致數(shù)據(jù)中冗余信息增多,當(dāng)下關(guān)鍵問題是,對高維的高光譜數(shù)據(jù)降維,達(dá)到簡化模型的目的,同時(shí)保持甚至提高模型的預(yù)測能力。

早期學(xué)者們利用相關(guān)性分析法[6-8](correlation analysis,CA)探究植物生理參數(shù)與其光譜反射率(或經(jīng)不同數(shù)學(xué)變換后的光譜反射率)的關(guān)系,選取相關(guān)系數(shù)高的波段作為敏感波長。但其只考慮了單波段與植物生理參數(shù)間的相關(guān)性,未考慮各波段間的共線性,難以解決光譜數(shù)據(jù)的冗余問題,而且所選波段較集中,只考慮某一段波長范圍的重要性,未考慮到其他波段,造成光譜數(shù)據(jù)的浪費(fèi)。而后學(xué)者們嘗試采用敏感變量優(yōu)選方法[9]從全波段內(nèi)剔除無關(guān)變量,優(yōu)選出敏感變量,減少數(shù)據(jù)量從而簡化模型。常用的變量優(yōu)選方法包括競爭性自適應(yīng)重加權(quán)算法(competitive adaptive reweighted sampling,CARS)、無信息變量消除算法(uninformative variable elimination, UVE)、移動窗口偏最小二乘法(moving window partial least square ,MWPLS)等。Li等[10]通過小波變換結(jié)合UVE技術(shù)簡化模型并提高了偏最小二乘回歸(PLSR) 模型預(yù)測的穩(wěn)定性。趙艷茹等[11]、邵園園等[12]利用CARS方法篩選敏感波段,簡化模型后得到比全波段還要高的PLSR預(yù)測精度。結(jié)果表明,對原始光譜進(jìn)行敏感波段優(yōu)選既可以降低模型復(fù)雜度,又能很好地提高模型的精度和穩(wěn)定性。隨機(jī)蛙跳算法(random frog,RF)通過在特征空間模擬一條平穩(wěn)分布的馬爾科夫鏈來計(jì)算每個(gè)變量被選擇的概率,從而進(jìn)行重要變量的篩選,被證明是一種較優(yōu)的變量優(yōu)選算法。如龍燕等[13]利用連續(xù)投影法(SPA)結(jié)合RF優(yōu)選出最佳波段,用于建立偏最小二乘回歸模型(PLSR)預(yù)測番茄的硬度;孫紅等[14]用CA和RF算法篩選到的敏感波段建立PLSR模型,結(jié)果表明,相比于CA法,RF算法篩選的敏感波段分布范圍更廣且對馬鈴薯葉片含水率預(yù)測性能更優(yōu);孫紅等[15]基于馬鈴薯葉片成像高光譜數(shù)據(jù),利用RF-PLSR模型反演出不同位置葉片逐像素點(diǎn)的葉綠素值。研究中大部分使用線性回歸(LR)、神經(jīng)網(wǎng)絡(luò)(NN)、偏最小二乘回歸(PLSR)等算法來建立回歸模型,各模型均有其特點(diǎn)和優(yōu)勢,因?yàn)镻LSR可有效解決高光譜數(shù)據(jù)波段間的共線性和信息冗余問題[16-18],應(yīng)用最為普遍。隨著機(jī)器學(xué)習(xí)方法[19]日漸成熟,被廣泛應(yīng)用于高光譜反演中。金秀等[20]優(yōu)選并組合了4個(gè)單模型,對集成算法進(jìn)行了優(yōu)化,結(jié)果顯示梯度提升樹算法(GBT)對土壤磷含量的預(yù)測精度最高。研究表明,采用機(jī)器學(xué)習(xí)算法可有效提高植物生理參數(shù)反演精度,明顯優(yōu)于傳統(tǒng)方法,但將機(jī)器學(xué)習(xí)方法融合的建模方法并不多見,與單個(gè)機(jī)器學(xué)習(xí)建模算法相比,融合建模方法對異常值和噪聲的敏感度更低,預(yù)測穩(wěn)定性能更優(yōu)。

本研究利用高光譜成像技術(shù)獲取紫丁香葉片光譜信息,針對葉片葉綠素含量基于隨機(jī)蛙跳(random frog,RF)方法篩選敏感波段,建立具有低復(fù)雜度和高穩(wěn)定性的投票回歸器(vote regressor,VR) 模型,并與全波段以及其他經(jīng)典變量提取方法篩選出的敏感波段建立的偏最小二乘回歸(partial least squares regression,PLSR)和投票回歸(VR)模型的預(yù)測結(jié)果進(jìn)行比較,同時(shí)結(jié)合偽彩色技術(shù)繪制紫丁香葉片葉綠素含量可視化分布圖,探索RF算法結(jié)合VR模型快速估測葉片葉綠素含量的可行性,以期為大面積監(jiān)測紫丁香冠層葉片養(yǎng)分分布和生長狀況提供技術(shù)支持。

1 材料與方法

1.1 試驗(yàn)樣本

實(shí)驗(yàn)對象為紫丁香,研究區(qū)選定在內(nèi)蒙古工業(yè)大學(xué)校園內(nèi)(呼和浩特市),研究對象為2020年5月采集于開花期的紫丁香葉片。根據(jù)校園紫丁香樹的分布情況,同時(shí)保證實(shí)驗(yàn)結(jié)果具代表性,在每棵紫丁香樹的東、西、南、北四個(gè)方位,隨機(jī)采集100個(gè)完整無損葉片樣本入袋密封并編號,后帶回實(shí)驗(yàn)室低溫冷藏的同時(shí)進(jìn)行實(shí)驗(yàn)測定。

1.2 丁香葉片高光譜圖像信息獲取

本文采用芬蘭Specim IQ高光譜成像系統(tǒng),一款帶有集成操作系統(tǒng)和控制裝置的手持式掃帚系統(tǒng),采集丁香葉片光譜成像數(shù)據(jù),結(jié)構(gòu)如圖1所示,該系統(tǒng)主要由高光譜相機(jī)、可控載物臺、植物葉片樣品、2個(gè)鹵素?zé)綦娫?、?jì)算機(jī)及相應(yīng)配套控制軟件組成。高光譜相機(jī)攝像頭分辨率為512×512像素,光譜范圍為400~1 000 nm,光譜分辨率為7 nm,設(shè)置載物臺和鏡頭之間的距離為20 cm,系統(tǒng)曝光時(shí)間為15 ms,為消除基線漂移須測量前預(yù)熱20 min,然后將已編號丁香葉片放于載物臺正對相機(jī),最終得到一個(gè)同時(shí)包含圖譜信息的三維數(shù)據(jù)塊。

1,樣品;2、3,鹵素電源;4,可控載物臺;5,高光譜相機(jī);6,數(shù)據(jù)傳輸線;7,計(jì)算機(jī);8,三腳架。1, Leaf sample; 2, 3, Light source; 4, Storage platform; 5, The Specim IQ hyperspectral camera; 6, Transmission data line; 7, Computer; 8, Tripod.圖1 高光譜成像系統(tǒng)Fig.1 Hyperspectral imaging monitoring system

在數(shù)據(jù)采集實(shí)驗(yàn)過程中,光照強(qiáng)度不均勻或暗電流等因素都會對實(shí)驗(yàn)結(jié)果產(chǎn)生影響,故需要對采集好的高光譜圖像進(jìn)行黑白板校正,最終得到校正后的原始光譜數(shù)據(jù)Rraw,校正公式如下:

(1)

式中:Rraw為黑白板校正后圖像數(shù)據(jù);W為白板數(shù)據(jù);B為黑板數(shù)據(jù);I為原始圖像數(shù)據(jù)。

1.3 葉綠素含量測定

本文采用手持式植物參數(shù)檢測儀對劃分區(qū)域進(jìn)行無損檢測,以SPAD值作為葉綠素含量參考指標(biāo)[21]。測量時(shí)避開葉脈和不平整區(qū)域,每片葉片主葉脈左右各選取3個(gè)感興趣區(qū)域(ROI),并對3個(gè)感興趣區(qū)域求平均,每片葉子可得2個(gè)SPAD值,最終通過對100個(gè)丁香葉片樣本的測量,獲得200個(gè)SPAD值。在測定SPAD值過程中對測量區(qū)域用馬克筆標(biāo)記測量范圍并編號,以便獲取相應(yīng)位置光譜。

1.4 光譜預(yù)處理

卷積平滑(savitzky golay,SG)濾波算法可以減少噪聲干擾,使光譜曲線更加平滑。光譜微分技術(shù)(spectral differentiation technology)通過計(jì)算光譜的n(n取1,2,3,…)階微分值來確定光譜曲線的極值點(diǎn)選取光譜響應(yīng)波段。應(yīng)用光譜微分技術(shù)能夠消除大氣效應(yīng)和植物背景的影響,將光譜曲線間的微小差異放大,可以更明顯地反映出不同葉綠素含量的植物的光譜響應(yīng)差異。故在波段篩選和建模前,選用卷積平滑(SG)和二階微分處理(second derivative,SD)對原始光譜數(shù)據(jù)進(jìn)行預(yù)處理[22-23],獲得SG-SD預(yù)處理后的葉片光譜反射率RSG-SD。

1.5 基于RF的特征波長提取方法

作為一種高效降維方法,隨機(jī)蛙跳[24-25](random frog,RF)算法在特征空間建立一條具有平穩(wěn)分布特性的馬爾科夫鏈,計(jì)算得到一個(gè)一維概率矩陣,每個(gè)概率值代表每個(gè)波段被選擇的概率大小。相比較于經(jīng)典變量優(yōu)選算法,該算法具有隨機(jī)搜索的特性,能夠利用較少的變量迭代建模。RF算法主要的運(yùn)算步驟包括以下4步:

(1)輸入一個(gè)初始波段子集F0,初始化時(shí)包含K個(gè)隨機(jī)波段,設(shè)定迭代次數(shù)N;

(2)在原始波段子集F0基礎(chǔ)上選出一個(gè)候選波段子集F*,包含K*個(gè)波段;對初始波段子集F0建立PLS模型,計(jì)算并降序排列各波段的絕對回歸系數(shù):若K*=K,則F*=F0;若K*K,前Q個(gè)波段構(gòu)成候選子集F*;

(3)選擇F*替代原波段子集F0,迭代N次后完成計(jì)算;

(4)計(jì)算N次迭代后每個(gè)波段被選擇的概率值,此概率值大小被作為變量是否被選取的標(biāo)準(zhǔn),概率值越大說明此波段越被優(yōu)先篩選。

1.6 基于VR的葉綠素含量預(yù)測方法

本研究提出的投票回歸器(vote regressor,VR)[26]是一種分步非參數(shù)方法,融合了線性模型和非線性模型兩大類建模方法,與傳統(tǒng)回歸模型和單個(gè)機(jī)器學(xué)習(xí)算法相比,能更好地處理偏離點(diǎn)和噪聲,平衡它們各自的弱點(diǎn)。建模流程如圖2,具體運(yùn)算步驟如下:

圖2 投票回歸器算法流程圖Fig.2 Voting regression algorithm flow chart

(1)多元線性回歸(multivariate linear regression,MLR)通過最小化誤差平方尋找最佳擬合函數(shù)預(yù)測葉綠素含量fMLR。

(2)利用梯度提升回歸(gradient boosting regression,GBR),通過串行地生成多個(gè)弱學(xué)習(xí)器,來擬合各分類器先前累加模型的損失函數(shù)的負(fù)梯度,使加上該弱學(xué)習(xí)器后的累積模型損失往負(fù)梯度的方向減少來預(yù)測葉綠素含量fGBT。

(3)利用隨機(jī)森林回歸[27](random forest regressor,RFR)模型,以決策樹為基學(xué)習(xí)器構(gòu)建Bagging,在生成樹的時(shí)候, 每個(gè)樹的每個(gè)節(jié)點(diǎn)都是隨機(jī)生成的,每個(gè)節(jié)點(diǎn)的拆分變量由少量隨機(jī)選擇的變量生成,形成多個(gè)決策樹,對結(jié)果取平均可得葉綠素含量的預(yù)測值fRFR。

(4)計(jì)算并返回三者的平均預(yù)測值,最終模型對葉綠素含量預(yù)測結(jié)果如式(2):

(2)

建立定量分析模型后,需選擇有效評估模型預(yù)測能力的指標(biāo),本研究選用決定系數(shù)R2進(jìn)行模型的評估,R2值在0和1之間,值越接近1,表明預(yù)測精度越高。本文數(shù)據(jù)處理與建模在Matlab2016和Python3.7軟件中完成。

2 結(jié)果與分析

2.1 丁香葉片光譜曲線分析

本文利用ENVI5.3提取丁香葉片編號區(qū)域的平均光譜值。Rraw曲線走勢總體符合綠色植物葉片的光譜響應(yīng)特性(圖3-A),主要特征包括在可見光的綠波段對葉綠素的強(qiáng)反射現(xiàn)象導(dǎo)致出現(xiàn)“波峰”現(xiàn)象,由于可見光紅波段對葉綠素的強(qiáng)吸收現(xiàn)象導(dǎo)致出現(xiàn)“波谷”現(xiàn)象,且在近紅外區(qū)(750~800 nm)光譜反射率急劇上升后不再變化,形成近紅外強(qiáng)反射平臺。從圖3-B中可以看出,Rraw經(jīng)SG-SD預(yù)處理后,由葉綠素含量差異導(dǎo)致的光譜曲線的等級差異得到有效消除,光譜曲線的微小細(xì)節(jié)特征被放大了。

圖3 紫丁香葉片樣本原始光譜曲線Rraw(A)和預(yù)處理后光譜曲線RSG-SD(B)Fig.3 original spectral curve Rraw (A) and pretreated spectral curve RSG-SD (B) of syringa oblata leaves

2.2 丁香葉片SPAD值統(tǒng)計(jì)

200條樣本SPAD值的分布范圍在18.3~44.3,平均值為32.4,其中SPAD值主要集中在22.9~41.3,采用SPXY(sample set partitioning based on joint X-Y distance)算法將200個(gè)樣本中160個(gè)樣本劃分為建模集,其余40個(gè)劃分為驗(yàn)證集,它是基于統(tǒng)計(jì)學(xué)角度的一種樣本集劃分方法,綜合考慮光譜和化學(xué)性質(zhì)的差異來選擇建模集,劃分結(jié)果如表1所示。

表1 樣本SPAD值統(tǒng)計(jì)與劃分

2.3 特征波長篩選

2.3.1 CA、CARS、MA-UVE、MWPLS算法特征變量篩選

全波段光譜數(shù)據(jù)量大且波譜間信息重疊現(xiàn)象嚴(yán)重,需要進(jìn)行敏感波段優(yōu)選以簡化模型提升模型效率。本研究利用隨機(jī)蛙跳算法(RF)和相關(guān)系數(shù)法(CA)、自適應(yīng)重加權(quán)算法(CARS)、無信息變量消除算法(UVE)、移動窗口偏最小二乘法(MWPLS)進(jìn)行敏感波段的篩選,而后對其篩選結(jié)果進(jìn)行建模,并對比和分析其建模精度。圖3為CA、CARS、UVE、MWPLS算法從SG-SD光譜數(shù)據(jù)中提取的特征波長在一條原始光譜曲線上的分布。

基于CA算法從RSG-SD中提取出相關(guān)系數(shù)絕對值大于閾值0.8的31個(gè)波段(圖4-A),CA算法的優(yōu)點(diǎn)是計(jì)算簡單,計(jì)算公式直觀且比較容易理解,但其選擇的特征波長較集中,所選敏感波段集中在485~710 nm,如果只考慮此段波長范圍的重要性,難以解決光譜數(shù)據(jù)的冗余問題,未考慮到其他波段,造成光譜數(shù)據(jù)的浪費(fèi)。

基于CARS算法從RSG-SD中共選擇29個(gè)波段(圖4-B),所選敏感波段集中在420~690 nm、800~920 nm,CARS算法選擇的敏感波段不穩(wěn)定,而且存在大量冗余信息,這些無用信息會影響重要變量的優(yōu)選。

圖4 CA(A)、CARS(B)、UVE(C)和MWPLS(D)算法篩選變量分布圖Fig.4 Distribution map of variables selected by CA(A), CARS(B), UVE(C) and MWPLS(D) method

基于UVE算法從RSG-SD中篩選敏感波段,通過交叉驗(yàn)證,RSG-SD的變量個(gè)數(shù)為40個(gè)(圖4-C)時(shí),模型精度達(dá)到最高,所選敏感波段集中在420~620 nm、700~750 nm、800~1000 nm。UVE可以盡可能消除無用波段,但其并未確定對紫丁香葉綠素敏感性強(qiáng)的波段,計(jì)算量大,耗時(shí)長,且篩選出的變量建模精度相對較低。

基于MWPLS 算法篩選最優(yōu)波段區(qū)間,但本研究中其篩選出波段數(shù)量為190個(gè)(圖4-D),所選區(qū)間為420~980 nm,只剔除掉14個(gè)波段,未解決光譜數(shù)據(jù)冗余問題。

2.3.2 RF算法特征變量篩選

隨機(jī)蛙跳(random frog,RF)算法與PLSR方法相結(jié)合,計(jì)算PLSR模型中所有變量的回歸系數(shù),將各變量絕對值大小作為迭代過程中每次該變量是否被選擇或者提出的依據(jù),后基于不同的波長點(diǎn)具有不同的概率值進(jìn)行敏感波段的選擇,運(yùn)行結(jié)果如圖5-A,橫軸代表波段,縱軸代表某波段被選擇的概率值,概率值越大說明該波段越重要。設(shè)定0.3作為葉片葉綠素含量對應(yīng)的篩選敏感波段的閾值,最終基于Rraw和RSG-SD分別篩選出49個(gè)和35個(gè)敏感波段,圖5-B為RF算法對RSG-SD篩選出的35個(gè)敏感波段在一條原始光譜曲線上的分布狀況,分布在420~450 nm、500~590 nm、620~650 nm、700~800 nm、850~900 nm、950~1 000 nm。RF算法既降低了波段間的多重共線性,又能更全面地提取與葉綠素含量相關(guān)的敏感波段,所選出的波段范圍更分散,跨度更廣。

圖5 RF運(yùn)行結(jié)果每個(gè)波段被選擇概率(A)和波長篩選結(jié)果(B)Fig.5 Results of of RF, Probability of each wavelength selected (A) and wavelength selected results (B)

2.4 PLSR建模

表2 不同變量篩選方法PLSR建模精度

2.5 VR建模

基于RF、CA、CARS、MWPLS和UVE算法選擇出的敏感波段和全波段(FULL),對Rraw和RSG-SD建立VR模型。表3中可看出,經(jīng)SG-SD處理后的光譜數(shù)據(jù)的建模精度較原始光譜數(shù)據(jù)均有不同程度提高,但對RSG-SD建立的CA-VR、RF-VR、CARS-VR、UVE-VR、MWPLS-VR模型的精度相比于FULL-VR模型提高不多,各精度值相差微小,說明采用VR建模對于葉綠素含量預(yù)測精度的提高效果不大,但建模過程中輸入的波段數(shù)卻大大減少,表明VR模型可以更好地解決變量間復(fù)雜的非線性關(guān)系,VR模型對異常值和光譜噪聲的敏感度更低,使模型預(yù)測穩(wěn)定性能更優(yōu)。在建模前對原始光譜數(shù)據(jù)進(jìn)行預(yù)處理和敏感變量篩選,在保證模型預(yù)測度的同時(shí)大大降低了模型的復(fù)雜度。圖7為對RSG-SD建立RF-VR模型后建模集和驗(yàn)證集樣本葉綠素含量實(shí)測值和預(yù)測值的散點(diǎn)圖。

圖6 RF-PLSR模型預(yù)測值和實(shí)測值散點(diǎn)圖Fig.6 Scatter diagram of predicted and measured values for the RF-PLSR model

表3 不同變量篩選方法VR建模精度

圖7 RF-VR模型預(yù)測值和實(shí)測值散點(diǎn)圖Fig.7 Scatter plot of predicted and measured values for the RF-VR model

2.6 紫丁香葉片葉綠素含量分布反演圖

由上得, RF波段優(yōu)選算法結(jié)合VR模型可有效預(yù)測紫丁香葉片上各個(gè)像素點(diǎn)的葉綠素含量。具體步驟如下:

(1)獲取敏感波段下的純紫丁香葉片高光譜圖像;

(2)提取圖像每個(gè)像素點(diǎn)的反射率;

(3)將(2)中結(jié)果代入RF-VR模型中求出各像素點(diǎn)SPAD值,得到灰度圖像;

(4)利用偽彩圖技術(shù)將灰度圖轉(zhuǎn)化為彩色圖,得到紫丁香葉片葉綠素分布圖(圖8)。

圖8可以直觀地看出紫丁香葉片上葉綠素的分布情況,偽彩色圖中顏色的差異和相同顏色深淺程度差異代表了紫丁香葉片中葉綠素濃度的差異。圖中,葉脈兩側(cè)葉綠素分布均勻,葉脈部分主要顯示為藍(lán)色(SPAD值為10~20),葉肉部分主要顯示為綠色(SPAD值為20~40),葉肉中的葉綠素含量較葉脈整體偏高。葉片首端顏色主要顯示紅色(SPAD值為40~50),葉片末端顯示黃色(SPAD值為30~40),首端葉綠素含量高于末端。圖中葉片外黑色區(qū)域?yàn)閿?shù)據(jù)采集背景,葉片內(nèi)黑色區(qū)域是由于葉片邊緣光照強(qiáng)度不均勻或陰影導(dǎo)致,故不可代表此區(qū)域SPAD值??傊?,最終根據(jù)RF-VR模型比較準(zhǔn)確地預(yù)測出了葉片葉綠素含量分布。

3 結(jié)論

綜上,利用高光譜成像技術(shù)結(jié)合光譜預(yù)處理技術(shù)和RF-VR模型能夠較好地對紫丁香葉片葉綠素含量反演和葉綠素分布可視化表達(dá)。但仍存在不足,本實(shí)驗(yàn)數(shù)據(jù)僅采集于紫丁香開花期,未考慮本模型是否同樣適用于紫丁香其他時(shí)期,采樣范圍也局限于校園內(nèi)。因此,今后應(yīng)擴(kuò)大采樣范圍和采樣時(shí)段以期提升模型精度和普適性。

猜你喜歡
紫丁香波段葉綠素
Ku波段高隔離度雙極化微帶陣列天線的設(shè)計(jì)
最佳波段組合的典型地物信息提取
新型X波段多功能EPR譜儀的設(shè)計(jì)與性能
紫丁香
最佳波段選擇的遷西縣土地利用信息提取研究
提取葉綠素
人教版小學(xué)語文二年級上冊第5課《一株紫丁香》教學(xué)設(shè)計(jì)
鋅對白菜和香蔥葉綠素含量的影響研究
綠茶保存要精心
紫丁香
东阿县| 曲阳县| 西充县| 兴安县| 邹城市| 齐河县| 手游| 招远市| 蒲城县| 固镇县| 客服| 赤峰市| 卢氏县| 图木舒克市| 美姑县| 高台县| 全南县| 政和县| 贵溪市| 泗水县| 浙江省| 喀喇沁旗| 台州市| 陈巴尔虎旗| 东丽区| 正安县| 安阳市| 九台市| 安塞县| 宁国市| 都兰县| 鄄城县| 宁强县| 山阴县| 东兴市| 久治县| 望奎县| 临朐县| 绍兴县| 马鞍山市| 夏津县|