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

基于三維數(shù)字化的小麥植株表型參數(shù)提取方法

2022-04-29 00:00:00鄭晨曦溫維亮盧憲菊郭新宇趙春江
智慧農(nóng)業(yè)(中英文) 2022年2期

摘要:針對小麥植株分蘗多、器官間交叉遮擋嚴(yán)重,難以用圖像或點云準(zhǔn)確提取植株和器官表型的問題,本研究提出了基于三維數(shù)字化的小麥植株表型參數(shù)提取方法。首先提出了小麥植株各器官數(shù)字化表達(dá)方法,制定了適用于小麥全生育期的三維數(shù)字化數(shù)據(jù)獲取規(guī)范,并依據(jù)該規(guī)范進(jìn)行數(shù)據(jù)獲取。根據(jù)三維數(shù)字化數(shù)據(jù)的空間位置語義信息和表型參數(shù)的定義,提出了小麥植株表型參數(shù)計算方法,實現(xiàn)了小麥植株和器官長度、粗度和角度等3類共11個常規(guī)可測表型參數(shù)的計算。進(jìn)一步提出了定量描述小麥株型和葉形的表型指標(biāo)。其中,植株圍度通過基于最小二乘法擬合三維離散坐標(biāo)計算,用于定量化描述小麥植株松散/緊湊程度;小麥葉片卷曲和扭曲程度為定量化葉形的指標(biāo),根據(jù)葉面向量方向變化計算得到。利用豐抗13號、西農(nóng)979號和濟(jì)麥44號三個品種小麥起身期、拔節(jié)期、抽穗期三個時期的人工測量值和提取值進(jìn)行驗證。結(jié)果表明,在保持植株原始三維形態(tài)結(jié)構(gòu)的前提下,提取的莖長、葉長、莖粗、莖葉夾角與實測數(shù)據(jù)精度相對較高,R2 分別為0.93、0.98、0.93、0.85;葉寬和葉傾角與實測數(shù)據(jù)的 R2分別為0.75、0.73。本方法能便捷、精確地提取小麥植株和器官形態(tài)結(jié)構(gòu)表型參數(shù),為小麥表型相關(guān)研究提供了有效技術(shù)支撐。

關(guān)鍵詞:小麥;三維數(shù)字化;可視化;表型參數(shù)提取

中圖分類號:S512.1;TP391.9"""" 文獻(xiàn)標(biāo)志碼:A"""""""" 文章編號:SA202203009

引用格式:鄭晨曦, 溫維亮, 盧憲菊, 郭新宇, 趙春江.基于三維數(shù)字化的小麥植株表型參數(shù)提取方法[J].智慧農(nóng)業(yè)(中英文), 2022, 4(2):150-162.

ZHENG Chenxi, WEN Weiliang, LU Xianju, GUO Xinyu, ZHAO Chunjiang. Phenotypic traits extraction of wheat plants using 3D digitization[J]. Smart Agriculture, 2022, 4(2):150-162.(in Chinese with English abstract)

1引言

植物表型是能夠反映植物細(xì)胞、組織、器官、植株和群體結(jié)構(gòu)及功能特性的物理、生理和生化性狀,其本質(zhì)是植物基因圖譜的時序三維表達(dá)及其地域分異特征和代際演進(jìn)規(guī)律[1]。植物表型研究是遺傳育種[2]和精準(zhǔn)栽培[3] 的關(guān)鍵,也是植物功能-結(jié)構(gòu)模型[4] 以及數(shù)字植物技術(shù)[5]等研究的核心內(nèi)容。隨著各學(xué)科的發(fā)展和數(shù)據(jù)處理能力的提升,植物表型研究已邁入了組學(xué)時代[6,7],傳統(tǒng)測量方式已無法滿足植物表型高通量獲取需求[8],迫切需要針對不同植物開展精準(zhǔn)、高效表型高通量獲取研究。小麥?zhǔn)鞘澜缛蠹Z食作物之一,開展小麥植株形態(tài)結(jié)構(gòu)表型參數(shù)提取方法研究對于小麥株型育種[9]具有重要意義。

傳統(tǒng)的植物形態(tài)結(jié)構(gòu)表型獲取以人工測量為主,所獲取表型參數(shù)精度受環(huán)境和人為因素影響較大。近年來,基于圖像的植物表型提取方法獲得廣泛應(yīng)用[10-13],但利用單幅圖像難以解決植物器官復(fù)雜的遮擋問題,部分形態(tài)結(jié)構(gòu)參數(shù)的提取精度受圖像獲取角度影響較大。相比之下,基于三維數(shù)據(jù)提取植物形態(tài)結(jié)構(gòu)表型可以得到較高的精度[14, 15]。在小麥三維表型提取方面,付晶波等[16]利用多視角圖像采集系統(tǒng)獲取小麥幼苗圖片并提取了小麥的株高和葉長,但僅針對小麥分蘗前的幼苗期進(jìn)行了試驗;翟苗苗[17] 利用 Li‐ DAR 獲取小麥植株三維數(shù)據(jù)并解析出高度和體積性狀,用于構(gòu)建小麥生物量模型。可以發(fā)現(xiàn),對于小麥這種分蘗多、葉片等器官遮擋嚴(yán)重的作物,利用常規(guī)的二維圖像和三維點云解析方法都難以提取到高精度的植株表型參數(shù),且無法利用一組數(shù)據(jù)提取到所需的多個表型參數(shù),尤其是小麥生育后期植株內(nèi)部的形態(tài)結(jié)構(gòu)信息。三維數(shù)字化儀是一種可以測量物體特征點精確空間坐標(biāo)位置的電磁設(shè)備,特別適用于分枝結(jié)構(gòu)復(fù)雜的植物[18, 19]。為解決器官間存在大量交叉遮擋植物形態(tài)結(jié)構(gòu)表型高精度提取問題,研究者利用三維數(shù)字化儀(3D Digitizer)直接還原植物植株在三維空間中形態(tài)結(jié)構(gòu),并在此基礎(chǔ)上提取植株表型參數(shù)。例如,溫維亮等[20,21]針對玉米形態(tài)結(jié)構(gòu)特征,制定了玉米三維數(shù)字化數(shù)據(jù)獲取標(biāo)準(zhǔn),并在此基礎(chǔ)上實現(xiàn)了玉米株型參數(shù)的準(zhǔn)確解析; Fournier等[22]利用三維數(shù)字化儀對小麥葉片出生方向、葉脈、葉形進(jìn)行準(zhǔn)確測量并構(gòu)建了數(shù)據(jù)集,并用于構(gòu)建小麥植株生長發(fā)育模型。然而,目前未見利用三維數(shù)字化儀對小麥植株進(jìn)行精確的形態(tài)結(jié)構(gòu)表型解析的報道。

本研究以此為出發(fā)點,針對不同生育期小麥植株形態(tài)結(jié)構(gòu)特征,利用三維數(shù)字化技術(shù),制定小麥植株三維數(shù)字化數(shù)據(jù)獲取規(guī)范,從“單株-單莖-器官”不同尺度一次性快速提取所需的多個小麥植株表型參數(shù)并進(jìn)行可視化。為更加深入分析小麥三維形態(tài)結(jié)構(gòu),本研究在玉米株型參數(shù)提取方法[21] 的基礎(chǔ)上,針對小麥的復(fù)雜分蘗結(jié)構(gòu)進(jìn)一步詳細(xì)了數(shù)據(jù)獲取規(guī)范。利用該規(guī)范獲取的小麥三維數(shù)字化數(shù)據(jù)相較于骨架信息,包含了葉片精細(xì)形態(tài)信息,能提取更多三維表型參數(shù)。同時,還對植株松散程度和葉片卷曲、扭曲程度進(jìn)行了定量化計算。在此基礎(chǔ)上,進(jìn)行小麥植株的下一步三維重建和表型解析能夠解決器官間遮擋嚴(yán)重的問題,從而獲得更準(zhǔn)確、完整的信息。

2小麥植株三維數(shù)字化數(shù)據(jù)獲取與表型參數(shù)提取方法

在不破壞小麥原始三維形態(tài)結(jié)構(gòu)的前提下,對大田種植的小麥植株進(jìn)行取樣,并依據(jù)制定好的數(shù)據(jù)獲取規(guī)范進(jìn)行小麥形態(tài)數(shù)據(jù)獲取。獲取完畢后,在 Python3.8環(huán)境下利用Open3D 0.13.0對不同生育時期的小麥植株進(jìn)行可視化。同時進(jìn)行相關(guān)表型參數(shù)的自動化提取,包括高度和長度參數(shù)、粗度參數(shù)、角度參數(shù),以及葉片卷曲和扭曲程度四類。圖1給出了小麥表型參數(shù)提取方法的整個流程。

2.1小麥三維數(shù)字化數(shù)據(jù)獲取

2.1.1數(shù)據(jù)獲取

冬小麥試驗地點為北京市農(nóng)林科學(xué)院試驗田( N39°56′,E116°16′)。選取了豐抗13號( FK13)、西農(nóng)979號( XN979)和濟(jì)麥44號( JM44)三個株型、葉形差異較大的品種。其中,豐抗13號植株較高,株型較松散,葉片具有細(xì)長的形態(tài)特征;濟(jì)麥44號植株較矮,株型較緊湊,葉片具有寬和短的特點;西農(nóng)979號株型和葉形特征均居中。播種時間為2020年10月4日,每個品種一個小區(qū),各小區(qū)大小設(shè)定為2.25 m×1.5 m,小區(qū)植株行距和株距分別為0.2和0.05 m 。田間播種前撒施復(fù)合肥( N-P2O5-K2O :12-18-15),施肥量為75 g/m2,并于拔節(jié)期隨水追施30 g/m2尿素。越冬期澆凍水一次,拔節(jié)期、灌漿期、成熟期分別澆水一次,灌水量為45 m3/m2。

采用機(jī)械臂式數(shù)字化儀 Microscribe" i ( Revware ,USA ,圖2)進(jìn)行數(shù)據(jù)獲取,其有效獲取范圍達(dá)1.27 m ,精度為0.178 mm ,滿足小麥植株三維數(shù)字化數(shù)據(jù)獲取要求。

于小麥生長的三個關(guān)鍵生育期(起身期、拔節(jié)期和抽穗期)對各品種植株進(jìn)行三維數(shù)字化數(shù)據(jù)獲取,每個品種重復(fù)3次。通過整株移植至花盆、搬至室內(nèi)的方式取樣,盡量保證所取樣品為葉片表面完整的植株,在無風(fēng)的環(huán)境下根據(jù)制定的數(shù)據(jù)獲取規(guī)范對小麥地上部進(jìn)行三維數(shù)字化數(shù)據(jù)采集。數(shù)據(jù)獲取前,基于電磁定位原理使用定標(biāo)塊進(jìn)行位置標(biāo)定,植株擺放在機(jī)械臂能觸達(dá)的范圍內(nèi)。為保持同一植株各器官在同一坐標(biāo)系下,數(shù)據(jù)獲取過程中植株不能移動。同時,對各器官形態(tài)結(jié)構(gòu)參數(shù)進(jìn)行人工測量并記錄。在數(shù)據(jù)采集過程中,要求無風(fēng)且儀器不能觸碰小麥植株,以避免小麥植株葉片在數(shù)據(jù)獲取過程中發(fā)生位置變化形成誤差。每株小麥植株數(shù)據(jù)獲取完成后,需利用自主開發(fā)的程序?qū)?shù)據(jù)以可視化的方式進(jìn)行檢查,若發(fā)現(xiàn)數(shù)據(jù)缺失或有誤,及時重新獲取。

2.1.2 小麥植株數(shù)字化表達(dá)及數(shù)據(jù)獲取規(guī)范

基于小麥生長發(fā)育規(guī)律制定適用于整個生育期的數(shù)據(jù)獲取規(guī)范,以確保整個植株三維數(shù)字化數(shù)據(jù)獲取的完整性和各結(jié)構(gòu)單元之間連接的正確性。曹衛(wèi)星和李存東[23]利用英文名稱首字母對小麥各器官命名,系統(tǒng)地提出了地上部器官序列化命名方案,且通過數(shù)字標(biāo)記表明植株上各器官發(fā)育的時空序位。本研究結(jié)合冬小麥的結(jié)構(gòu)單元劃分,在此基礎(chǔ)上以器官為單位,制定了小麥植株各器官數(shù)字化表達(dá)及數(shù)據(jù)獲取規(guī)范。

(1) 莖。由 S ( Stem )表示,主莖( MainStem )為 SM ;分蘗( Tiller )為 STi,j , i 為著生在主莖上一級分蘗序號,j為著生在一級分蘗上的二級分蘗序號。由于一般三級分蘗存活率不高,本研究暫不考慮。首先在主莖出生處環(huán)繞式獲取3個坐標(biāo)點用來記錄莖的粗度信息,然后從莖的出生處開始由下至上獲取若干點(圖3(a)),莖上最后一片葉與莖的連接處作為結(jié)束點。小麥拔節(jié)后,要求每個點與小麥的節(jié)點位置重合。抽穗后,需以穗的出生處作為結(jié)束點。主莖數(shù)據(jù)點集由{Sk(M)}表示,N SM 為所獲取主莖三維數(shù)字化點的數(shù)量,因此0lt; k ≤ N SM,k ∈ Z ,Z 為整數(shù);將主莖出生點 S4(M)視為整個小麥植株的出生點。分蘗數(shù)據(jù)點集由{Sk(T)i,j }表示,Ni,j(ST)為對應(yīng)分蘗包含數(shù)據(jù)點的數(shù)量。

(2) 葉。由 L ( Leaf)表示,主莖葉片為LM ;分蘗葉片為 LTi,j 。按照從葉片出生處沿葉尖方向,逐行獲取葉片左邊緣、葉脈、右邊緣三個點(圖3(b)) 數(shù)據(jù),以葉尖處一個點作為結(jié)束。要求間隔距離不宜過長,需包含能明顯表示葉片形態(tài)變化處以及葉片最高點信息,葉片發(fā)生扭曲時仍須以葉片平展開之后的初始橫向方向進(jìn)行數(shù)據(jù)獲取。主莖葉片數(shù)據(jù)點集由{Lk(M),n }表示,n 為葉片生長發(fā)育順序, Q0為主莖包含葉片數(shù)量,因此0lt; n ≤ Q0,n ∈ Z ;Nn(L)M 為所獲取主莖第 n 片葉三維數(shù)字化點的數(shù)量,因此0lt; k ≤ Nn(L)M,k ∈ Z ;分蘗葉片數(shù)據(jù)點集由{L k(T)i,j,n }表示,Qi,j 為對應(yīng)分蘗包含的葉片數(shù)量;Ni,j,(LT)n 為對應(yīng)分蘗第 n片葉包含數(shù)據(jù)點的數(shù)量。

(3) 穗。由 E ( Ear )表示,主莖的穗為 EM;分蘗的穗為 ETi,j 。數(shù)據(jù)獲取從穗的出生處開始由下至上獲取若干點,以穗尖處作為結(jié)束點。主莖的穗數(shù)據(jù)點集由{Ek(M)}表示,NEM 為所獲取主莖穗的三維數(shù)字化點的數(shù)量,因此0lt; k ≤ NEM,k ∈ Z ;分蘗穗的數(shù)據(jù)點集由{Ek(T)i,j }表示,N 為對應(yīng)分蘗的穗包含數(shù)據(jù)點的數(shù)量。

2.2小麥表型參數(shù)提取方法

2.2.1 數(shù)據(jù)標(biāo)準(zhǔn)化

數(shù)據(jù)可視化和提取小麥植株表型參數(shù)前,需對三維數(shù)字化數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理,從而使小麥植株均位于原點,且主莖生長方向豎直向上重合于 Z軸,便于植株的可視化顯示以及高度參數(shù)的提取。

由于主莖數(shù)據(jù)點集前3個數(shù)據(jù)點為粗度信息,第4個點 S4(M)為主莖出生點。首先將 S4(M)平移至坐標(biāo)(0,0,0),并對所有剩余數(shù)據(jù)點做同等平移操作。然后根據(jù)主莖生長方向?qū)χ仓赀M(jìn)行旋轉(zhuǎn),生長方向由主莖上第一個和最后一個記錄長度信息的坐標(biāo)點形成的向量表示," = SNS(M)M .z -Si(M) .x、Si(M) .y、Si(M) .z 分別表示三維空間點 x、y、z 坐標(biāo)值。記θ為與今v-z 的夾角,以向量今axis = ×今v-z作為旋轉(zhuǎn)軸,所有數(shù)據(jù)點繞今axis 旋轉(zhuǎn)角度θ完成數(shù)據(jù)標(biāo)準(zhǔn)化。完成標(biāo)準(zhǔn)化之后,針對不同數(shù)據(jù)點集使用參數(shù)提取方法進(jìn)行表型參數(shù)的計算與提取。

2.2.2 高度、長度參數(shù)提取

(1) 株高、莖長、節(jié)間長的計算。根據(jù)數(shù)據(jù)獲取規(guī)則,拔節(jié)后莖的數(shù)據(jù)點的獲取與節(jié)點一致。d ( pi,pj )表示三維空間中兩點的歐式距離,以主莖為例,則從主莖出生處開始,第 i個節(jié)間長的計算方法如公式(1) 所示。node為節(jié)間長度,cm。

莖長為各節(jié)間長的總和,計算方法如公式(2)所示。s為莖長,cm。

數(shù)據(jù)標(biāo)準(zhǔn)化后,三維空間內(nèi)的 Z坐標(biāo)就表示高度信息。小麥植株株高為主莖和分蘗的高度最大值,抽穗前,高度最大值由最上部葉葉尖點的Z坐標(biāo)決定。NQ(L)0(M)為主莖最上部葉數(shù)據(jù)點的數(shù)量,令 n0 = N" 則 Ln0(M),Q0為主莖最上部葉的最后一個坐標(biāo)點;Ni,j,(LT)Qi,j 為分蘗 Ti,j 的最上部葉的最后一個坐標(biāo)點,令 ni,j = Ni,j,(LT)Qi,j ,則 L" 為分蘗最上部葉的最后一個坐標(biāo)點,株高的計算方法如公式 (3) 所示。H為株高,cm。

抽穗后,高度最大值根據(jù)旗葉姿態(tài)決定,選擇旗葉和穗最后一個點的 Z 坐標(biāo)的較大值。ENE(M)M 為主莖上穗的最后一個坐標(biāo)點,EN(T) 為分蘗 Ti,j 上穗的最后一個坐標(biāo)點,則株高的計算方法如公式 (4) 所示。

(2) 葉長、葉寬、葉片出生位置的計算。數(shù)據(jù)獲取時,葉片每行數(shù)據(jù)點的間隔距離保證了曲率沒有較大變化,因此葉長采用計算葉脈上相鄰數(shù)據(jù)點間距離和的方式。但折線段的累積和小于葉長實際值,根據(jù)不同品種實際值和提取值確定經(jīng)驗比例系數(shù)ξ調(diào)整計算公式。由于集合內(nèi)葉脈數(shù)據(jù)點的間隔為3,最后葉尖處與前一點的間隔為2,以主莖第一片葉為例,數(shù)據(jù)點集合大小為 N" 則葉長的計算方法如公式(5) 所示。l為葉長,cm。

葉寬的計算為左、右邊緣點到葉脈點的距離之和,以和葉長同樣方式確定經(jīng)驗比例系數(shù)τ。以主莖第一片葉為例,從葉片出生處開始,第 i 排位置處葉寬的計算方法如公式(6) 所示。w 為葉寬,cm。

其最大葉寬 w為 max (wi )。

葉片的出生位置即該葉片與莖的連接點的位置,具體為第1行數(shù)據(jù)點的中間葉脈點( L x,L y,L z ),且出生高度為 L z。

(3) 穗長、穗出生位置的計算。穗長也采用穗上相鄰數(shù)據(jù)點間距離和的方式計算,以主莖為例,穗長的計算方法如公式(7) 所示, e 為穗長,cm。

穗的出生位置為穗的數(shù)據(jù)點中的第一個坐標(biāo)點,具體為(E 1(M) .x,E 1(M) .y,E 1(M) .z ),且出生高度為E1(M) .z。

2.2.3 粗度參數(shù)提取

(1) 莖粗的計算。通過計算莖數(shù)據(jù)點集前3個點的外接圓半徑確定莖粗。根據(jù)三點共面、三點到圓心距離相等兩個約束條件求得3個平面方程如(8) 所示。

其中,A 為表示方程系數(shù)的三階方陣; B 為常數(shù)項;P為外心坐標(biāo)矩陣。以主莖為例,通過S 1(M)、S2(M)、S3(M)三點可以求得外心po" (xo ,yo ,zo ),則莖粗 o 為 S 1(M)、S2(M)、S3(M)三點到圓心距離的兩倍加權(quán)平均值,cm ,計算方法如公式(9) 所示。

(2) 圍度的計算。小麥植株的緊湊和松散程度與分蘗產(chǎn)生后莖稈下半部分生長情況有著直接聯(lián)系,將小麥下部1/3處各莖組成的多邊形外接圓的直徑定義為小麥的圍度。株高為H的小麥植株通過擬合 H/3高度處各莖共 m 個坐標(biāo)點pi ,由于坐標(biāo)點不在確定的圓周上,所以利用拉格朗日乘數(shù)法求解圓心的最優(yōu)解[24]來計算小麥植株的圍度 u (圖4)。圍度大小可直接說明小麥植株的松散程度,定量描述小麥的株型。

通過構(gòu)建出拉格朗日函數(shù) F (Pu ,λ),分別對 Pu ,λ求一階偏導(dǎo)數(shù),并令 F’(Pu ,λ)等于零,解出 Pu = [xu yu" zu ] T 即為最后圓心坐標(biāo)的最優(yōu)解,可構(gòu)建方程如公式(10)所示。

其中,D 和 G 分別為表示線性方程組f (xu,yu ,zu )方程系數(shù)的×3階矩陣和常數(shù)項;f (xu,yu ,zu )為在假設(shè)所有坐標(biāo)點都在圓周上此前提下,利用其中任意兩點的連線垂直于連線中點與圓心連線此性質(zhì)所構(gòu)建的×3維線性方程組; U = (a, b, c ) T,其中 a , b , c 為所有離散點(包括圓心)所在的平面 ax + by + cz -1=0的方程系數(shù)。

根據(jù)各莖坐標(biāo)點pi 到圓心pu" (xu,yu ,zu )的兩倍加權(quán)平均距離確定圍度 u ,計算方法如公式(11)所示。

2.2.4角度參數(shù)提取

(1) 莖葉夾角的計算。根據(jù)莖葉夾角定義并參考人工測量方法,可通過葉片出生處生長方向和莖局部生長方向這兩個空間向量的夾角計算莖葉夾角。葉片出生處生長方向可由前2個葉脈數(shù)據(jù)點來表示,莖局部生長方向通過葉片出生高度定義此葉片在第幾節(jié)間,計算節(jié)點間向量。則莖葉夾角的計算如公式(12)所示,α為莖葉夾角,(°)。

(2) 葉傾角的計算。一般情況下,通過直接連接葉耳與葉尖,測量此向量與 XOY 平面的夾角確定葉傾角。但小麥葉片在不同時期存在較大的差異性,彎曲程度很大,所以需要分段計算。根據(jù)數(shù)據(jù)獲取規(guī)范,葉片最高點一定在數(shù)據(jù)點集合中。計算相鄰葉脈點連成的向量和 XOY 平面之間的夾角取平均值,為最后的葉傾角結(jié)果(圖5)。葉片后期受重力原因或衰老出現(xiàn)下垂情況,不納入葉傾角計算范圍內(nèi)。以主莖第一片葉為例,葉片數(shù)據(jù)點集合大小為 N 所以(N 1)/3段葉脈向量中只取 Z為正值的向量。如圖5所示,葉片最高點之前的有效分段向量(顏色標(biāo)記段)納入計算范圍。假設(shè)前 m段葉脈向量符合要求,取第一段葉脈向量在 XOY 平面上的投影向量 e-z,則葉傾角的計算如公式(13)所示,β為葉傾角,(°)。

2.2.5葉片卷曲、扭曲程度參數(shù)提取

在小麥生長過程中,部分葉片會發(fā)生扭曲、卷曲等形變。通過參數(shù)對葉片扭曲、卷曲程度進(jìn)行定量化,有利于分析不同品種間葉片的形態(tài)差異。

(1) 葉片卷曲程度的計算。葉面卷曲程度刻畫以葉脈為中心,左右兩側(cè)葉面形成的卷曲。通過提取葉面有關(guān)向量,定義小麥扭曲和卷曲程度參數(shù)。通過右手定則計算葉面數(shù)字化數(shù)據(jù)中葉脈左右面元法向 u left和 u right的夾角γ(圖6),從而確定此葉長處的卷曲程度。以主莖第一片葉為例,葉片數(shù)據(jù)點集合大小為N 則能確定(N 1)/3個葉脈左右面元夾角,卷曲程度c的計算如公式(14)所示。

c取值范圍在[0, 1],值越大證明葉片卷曲越明顯。

(2) 葉片卷曲程度的計算。葉面扭曲程度刻畫沿葉脈從基部到葉尖方向葉面方向的整體變化趨勢。未發(fā)生扭曲的小麥葉片,對應(yīng)的左右邊緣點確定的向量在三維空間中應(yīng)該具有平行的性質(zhì)。通過計算相鄰向量的夾角累積和(圖7),可以確定葉片的扭曲程度。以主莖第一片葉為例,葉片數(shù)據(jù)點集合大小為N 則能確定(N 1)/3

個向量,向量計算如公式(15)所示。

扭曲程度 q的計算方法如公式(16)所示。

q取值范圍在[0, 1],數(shù)值越大表明葉片扭曲越明顯。圖8給出了對小麥葉片數(shù)據(jù)點集進(jìn)行可視化的結(jié)果,其扭曲度分別為0.12、0.56和0.89。

2.3數(shù)據(jù)可視化

用 Open3D 庫對小麥植株數(shù)據(jù)進(jìn)行處理并可視化,數(shù)據(jù)處理流程可以分為以下3個步驟。

(1) 三角網(wǎng)格化。利用數(shù)字化數(shù)據(jù)包含點的坐標(biāo)和順序信息,生成三角網(wǎng)格。

(2) 網(wǎng)格細(xì)分。對三角網(wǎng)格化后的葉片,采用 Loop網(wǎng)格細(xì)分[25]方法對各葉片網(wǎng)格加細(xì),針對小麥葉片迭代細(xì)分2次后可達(dá)到最佳效果。

(3) 網(wǎng)格著色。根據(jù)小麥器官顏色使用paint_uniform_color函數(shù)對網(wǎng)格進(jìn)行上色。

圖9分別給出了小麥單個葉片和整個植株原始數(shù)字化數(shù)據(jù)、生成三角網(wǎng)格、進(jìn)行網(wǎng)格細(xì)分,和網(wǎng)格著色的效果圖。

2.4數(shù)據(jù)分析方法

采用以下幾個標(biāo)準(zhǔn)評價參數(shù)精度。

(1) 均方根誤差(Root Mean Square Error,RMSE )是實測值與預(yù)測值之間誤差平均值的平方根,計算方法如公式(17)所示。

(2) 決定系數(shù)( R-square ,R2)是實測值與預(yù)測值之間的模型擬合度,取值范圍在[0, 1],計算方法如公式(18)所示。

其中,Mean 為實測值均值,計算方法如公式(19)所示。

3結(jié)果與分析

3.1可視化結(jié)果

利用 Open3D 對小麥三維數(shù)字化數(shù)據(jù)進(jìn)行可視化,以最直觀的方式體現(xiàn)小麥植株品種間以及不同生育期的差異,在此基礎(chǔ)上進(jìn)行小麥模型構(gòu)建,可解決小麥器官遮擋嚴(yán)重的問題。

圖10給出了在起身期、拔節(jié)期、抽穗期獲取的3個品種數(shù)字化數(shù)據(jù)的可視化結(jié)果。由結(jié)果可知, JM44分蘗最多, XN979次之, FK13最少。在株型方面,XN979在拔節(jié)前植株相對更高;但拔節(jié)后 FK13節(jié)間迅速伸長,直至形態(tài)不再發(fā)生變化時,F(xiàn)K13植株均為最高,XN979相對較矮;同時 FK13圍度值較大,更加松散; XN979和 JM44圍度值較小,更加緊湊。在葉片形態(tài)方面,F(xiàn)K13葉片具有細(xì)長的特點,葉長最大,葉寬最小,卷曲程度較大;XN979和 JM44葉寬較大,卷曲程度較小;但由于品種特性, JM44葉片擁有更大的葉片扭曲程度,XN979葉片最為平展。

3.2參數(shù)提取結(jié)果精度分析

為評價所提出的小麥植株表型參數(shù)提取方法的精度,選取了起身期、拔節(jié)期、抽穗期對應(yīng)數(shù)據(jù)進(jìn)行驗證。其他方法無法精準(zhǔn)提取到各時期所需表型參數(shù),因此利用實測數(shù)據(jù)進(jìn)行對比分析。表1給出了抽穗期3個品種小麥植株長度、粗度和角度參數(shù)的提取值和實測值均值。由結(jié)果可知,提取值與實測值相差不大,且能體現(xiàn)出品種差異。其中, FK13葉長值最大,葉寬值最小;JM44葉片擁有相對較小的葉長值,葉寬值最大;XN979居中。在莖長、莖粗方面,F(xiàn)K13莖長且細(xì),XN979莖偏短且粗,JM44居中。FK13葉片的莖葉夾角較大,但葉傾角最小;JM44莖葉夾角最小,但葉傾角最大;XN979居中。圖 11進(jìn)一步給出了3個品種小麥植株在3個生育期提取的葉長、葉寬、莖長、莖粗、莖葉夾角、葉傾角值與實測值對比結(jié)果,以及驗證結(jié)果均方根誤差RMSE和決定系數(shù) R2。圖中粉色、黑色和藍(lán)色分別表示起身期、拔節(jié)期和抽穗期三個生育期的擬合效果。

在小麥表型參數(shù)提取方面,其他方法都無法在保持植株原有三維形態(tài)結(jié)構(gòu)的前提下,進(jìn)行小麥以上表型參數(shù)的完整、準(zhǔn)確獲取[26];特別在三葉一心后的主要關(guān)鍵生育時期,小麥形態(tài)結(jié)構(gòu)復(fù)雜,器官遮擋嚴(yán)重。由對比結(jié)果可知,葉長、莖長、莖粗、莖葉夾角的計算誤差較小,三個品種的平均 R2分別為0.93、0.98、0.93、0.85,葉寬和葉傾角與實測值的相關(guān)系數(shù) R2分別為0.75、0.73,但在不進(jìn)行植株拆除測量的情況下,誤差不大。

3.3參數(shù)提取時間效率分析

為評價本方法的重建效率,計算了單個小麥植株數(shù)據(jù)處理和參數(shù)提取的 CPU 所用時間。以抽穗期小麥為例,其平均分蘗數(shù)和平均葉片數(shù)分別為7和27,在此情況下數(shù)據(jù)讀取和處理所耗時間約為121 ms,表型參數(shù)提取所耗時間約為12 ms,總時間為133 ms 。以上方法是在 windows10操作系統(tǒng)下 VS2019中實現(xiàn),使用的計算機(jī)配置為: Intel" ( R )" Core" ( TM )" i7-10700 CPU @2.90GHz ,16G RAM。

4討論

本研究針對小麥分蘗多、器官遮擋嚴(yán)重等特征,基于三維數(shù)字化技術(shù)構(gòu)建了一套快速、精確、自動提取小麥表型參數(shù)的方法,具有一致性強(qiáng)、精度高、操作便捷的特點。此方法在保持植株原始三維形態(tài)結(jié)構(gòu)的前提下,對小麥三維數(shù)字化數(shù)據(jù)進(jìn)行可視化,并一次性計算提取小麥包含長度、粗度和角度的11個常規(guī)可測表型參數(shù)。同時,定量分析小麥植株的松散、緊湊程度和葉片卷曲、扭曲程度,分別從植株和器官等不同尺度更好地描述了小麥的形態(tài)特征。此方法有助于分析小麥品種的形態(tài)結(jié)構(gòu)差異,對小麥株型理想株型培育具有一定意義。

利用三維數(shù)字化儀獲取單株小麥數(shù)據(jù),獲取工作量與對所有參數(shù)進(jìn)行人工測量相近。但利用三維數(shù)字化數(shù)據(jù)可對如植株圍度、葉片卷曲和扭曲程度等參數(shù)定量化,人工則難以定量測量。同時,人工測量主觀性強(qiáng),不同操作、個人讀數(shù)習(xí)慣等都容易影響測量結(jié)果,而三維數(shù)字化數(shù)據(jù)更為客觀,不因上述因素產(chǎn)生數(shù)據(jù)誤差,具有更好的一致性。

基于圖像和基于三維點云的植物株型參數(shù)提取是目前應(yīng)用較多的兩類方法,可實現(xiàn)玉米[12]、棉花[13]、甜菜[14]等植物表型參數(shù)的提取。但由于小麥屬于多分蘗作物,具有分蘗和葉片多、葉片細(xì)軟的特點,同時葉片相互遮擋嚴(yán)重,基于圖像和點云的方法都無法適用于小麥整個生育期,且提取的表型參數(shù)有限。目前未見其他方法在不破壞植株原始三維形態(tài)結(jié)構(gòu)的情況下,自動批量提取本系統(tǒng)中包含的所有小麥表型參數(shù)。在已有的小麥表型參數(shù)提取研究中,利用多視角圖像系統(tǒng)采集的小麥幼苗圖片,結(jié)合 Mask R-CNN 方法[16]提取了小麥葉長和株高,相關(guān)系數(shù) R2分別為0.87和0.98。本方法提取的小麥葉長和莖長 R2 分別為0.93和0.98,誤差更小。基于 LiDAR的方法[17]對 2016—2018年連續(xù)兩年小麥的分蘗數(shù)進(jìn)行了估算,相關(guān)系數(shù) R2分別為0.61和0.56,同時利用高度指標(biāo)和體積指標(biāo)構(gòu)建小麥生物量模型,但缺乏長度和角度等更精準(zhǔn)的表型參數(shù)的獲取。本方法可以直接根據(jù)數(shù)據(jù)集數(shù)量得知分蘗數(shù),并提取長度、粗度、角度相關(guān)的表型參數(shù)。提取的葉長、莖長、莖粗、莖葉夾角的 R2分別為0.93、0.98、0.93、0.85,葉寬和葉傾角的 R2分別為0.75、0.73。其中,葉寬和葉傾角的誤差相較其他表型參數(shù)略大。主要原因是小麥葉片較窄,同時葉片會出現(xiàn)卷曲和扭曲現(xiàn)象,所獲取的葉邊緣點可能不在葉脈的垂線上,人工測量時也不易確定最大葉寬處。葉傾角誤差較大是由于人工測量對測量位置選取不一致導(dǎo)致測量精度不高;而利用三維數(shù)字化數(shù)據(jù)計算的葉傾角是多處提取結(jié)果的均值,各種彎曲程度葉片都能實現(xiàn)精度較高的計算結(jié)果。

本方法制定了包含植株語義信息的數(shù)據(jù)獲取規(guī)范,利用三維數(shù)字化儀進(jìn)行數(shù)據(jù)獲取精確定位了各個莖、葉片、穗等器官的空間位置,解決了小麥器官間遮擋問題,可以實現(xiàn)小麥表型參數(shù)的精準(zhǔn)提取。該數(shù)據(jù)獲取規(guī)范可直接應(yīng)用到水稻等形態(tài)結(jié)構(gòu)相似的禾本科分蘗作物上,其他作物也可以在此基礎(chǔ)上根據(jù)自身形態(tài)結(jié)構(gòu)特征進(jìn)行修改和應(yīng)用。

同時,本研究實現(xiàn)了小麥植株的可視化仿真,能完整還原小麥三維形態(tài)結(jié)構(gòu)并進(jìn)行表型參數(shù)的精確提取。雖然已有研究實現(xiàn)了小麥三維模型的快速重建,但由于利用三維掃描或多視角三維重建得到的植株點云存在如葉緣和植株內(nèi)部信息缺失問題[27,28],不適用于復(fù)雜植株的三維重建,且難以高精度地提取器官尺度表型參數(shù)[29,30]。本方法獲取數(shù)據(jù)精確,但依賴于儀器設(shè)備,數(shù)據(jù)獲取的時間成本也相對較高,對于實現(xiàn)高通量的小麥表型參數(shù)獲取還需進(jìn)一步研究。

參考文獻(xiàn):

[1] 趙春江.植物表型組學(xué)大數(shù)據(jù)及其研究進(jìn)展[J].農(nóng)業(yè)大數(shù)據(jù)學(xué)報, 2019, 1(2):5-18.

ZHAO C. Big data of plant phenomics and its researchprogress[J]. Journal of Agricultural Big Data, 2019, 1(2):5-18.

[2] VOS J, EVERS J B, BUCK-SORLIN G H, et al. Func‐tional-structural plant modelling: A new versatile toolin" crop" science[J]. Journal" of" Experimental" Botany,2010, 61(8):2101-2115.

[3] ROBERT T, FURBANK, TESTERMARK. Phenomics-technologies to relieve the phenotyping bottleneck[J].Trends in Plant Science, 2011, 16(12):635-644.

[4] COBB" J" N," DECLERCK" G," GREEBERG" A," et "al.Next-generation phenotyping: Requirements and strate‐ gies for enhancing our understanding of genotype-phe‐ notype relationships and its relevance to crop improve‐ ment[J]. Theoretical and Applied Genetics, 2013, 126:867-887.

[5] 趙春江, 陸聲鏈, 郭新宇, 等.數(shù)字植物及其技術(shù)體系探討[J].中國農(nóng)業(yè)科學(xué), 2010, 43(10):2023-2030.

ZHAO C, LU S, GUO X, et al. Exploration of digital plant and its technology system[J]. Scientia Agricultura Sinica, 2010, 43(10):2023-2030.

[6] FIORANI F," SCHURR U. Future" scenarios" for plantphenotyping[J]. Annual" Review" of" Plant" Biology, 2013, 64(1):267-291.

[7] HOULE D, GOVINDARAJU D R, OMHOLT S. Phe‐nomics: The next challenge[J]. Nature Reviews Genet‐ ics, 2010, 11(12):855-866.

[8] KUMAR" J," PRATAP A," KUMAR" S. Phenomics" incrop plants: Trends, options and limitations[M]. Berlin: Springer India, 2015.

[9] GRO?KINSKY DOMINIK K, JESPER S, SVEND C,et al. Plant phenomics and the need for physiological phenotyping" across" scales to narrow the genotype-to- phenotype knowledge gap[J]. Journal of Experimental Botany, 2015, 66(18):5429-5440.

[10]朱榮勝, 李帥, 孫永哲, 等.作物三維重構(gòu)技術(shù)研究現(xiàn)狀及前景展望[J].智慧農(nóng)業(yè)(中英文), 2021, 3(3):94-115.

ZHU R, LI" S," SUN Y," et" al. Research" advances" and prospect" of" crop 3D" reconstruction" technology[J]. Smart Agriculture, 2021, 3(3):94-115.

[11] GUO W, RAGE U K, NINOMLYA S. Illumination in‐variant" segmentation" of" vegetation" for" time" series wheat images based on decision tree model[J]. Com ‐ puter and Electronics in Agriculture, 2013, 96:58-66.

[12]宗澤, 張雪, 郭彩玲, 等.基于骨架提取算法的作物表型參數(shù)提取方法[J].農(nóng)業(yè)工程學(xué)報 , 2015, 31(S2):180-185.

ZONG Z, ZHANG X, GUO C, et al. Crop phenotypic parameters extraction method based on skeleton extrac‐ tion algorithm[J]. Transactions of the CSAE, 2015, 31(S2):180-185.

[13] PAPROKI A, SIRAULT X, BERRY S, et al. A novelmesh" processing" based" technique" for" 3D" plant analysis[J]. BMC Plant Biology, 2012, 12(1):63-63.

[14]柴宏紅, 邵科, 于超, 等.基于三維點云的甜菜根表型參數(shù)提取與根型判別[J].農(nóng)業(yè)工程學(xué)報 , 2020, 36(10):181-188.

CHAI H, SHAO K, YU C, et al. Extraction of pheno‐ typic parameters and discrimination of beet root types based on 3D point cloud[J]. Transactions of the CSAE, 2020, 36(10):181-188.

[15] PAULUS S, SCHUMANN H, KUHLMANN H, et al.High-precision laser scanning system for capturing 3Dplant" architecture" and" analysing" growth" of" cerealplants[J]. Biosystems Engineering, 2014, 121:1-11.

[16]付晶波, 施家偉, 張俊, 等.基于 Mask R-CNN的盆栽小麥單片葉長和株高提取研究[J].吉林農(nóng)業(yè)大學(xué)學(xué)報, 2021, 43(2):163-170.

FU J, SHI J, ZHANG J, et al. Extraction of leaf lengthand plant height from potted wheat based on Mask R-CNN [J]. Journal of Jilin Agricultural University, 2021,43(2):163-170.

[17]翟苗苗.基于地面激光雷達(dá)的小麥生物量和分蘗數(shù)估算研究[D].南京:南京農(nóng)業(yè)大學(xué), 2018.

ZHAI M. Estimating wheat biomass and tiller numberbased" on terrestrial" laser" scanning[D]. Nanjing: Nan‐jing Agricultural University, 2018.

[18] SINOQUET" H," "THANISAWANYANGKURA"" S,MABROUK H, et al. Characterization of the light envi‐ronment in canopies using 3D digitising and image pro‐cessing[J]. Annals of Botany, 1998, 82(2):203-212.

[19]郭焱 , 李保國.虛擬植物的研究進(jìn)展[J].科學(xué)通報 ,2001(4):273-280.

GUO Y," LI" B. Research" progress" of virtual" plant[J].Chinese Science Bulletin, 2001(4):273-280.

[20]溫維亮, 郭新宇, 盧憲菊, 等.玉米器官三維模板資源庫構(gòu)建[J].農(nóng)業(yè)機(jī)械學(xué)報, 2016, 47(8):266-272.

WEN W, GUO X, LU X, et al. Three-dimensional tem ‐plate resource library construction of maize organs[J].Transactions of the CASM, 2016, 47(8):266-272.

[21]溫維亮, 郭新宇, 趙春江, 等.基于三維數(shù)字化的玉米株型參數(shù)提取方法研究[J].中國農(nóng)業(yè)科學(xué) , 2018, 51(6):1034-1044.

WEN W, GUO X, ZHAO C, et al. Research on maizeplant type parameter extraction by using three dimen‐sional" digitizing" data[J]. Scientia" Agricultura" Sinica,2018, 51(6):1034-1044.

[22] FOUTNIER" C," ANDRIEU" B," LJUTOVAC" S," et" al.ADEL-Wheat: A 3D Architectural Model of wheat de‐velopment[C]// International" Symposium" on" PlantGrowth Modeling. Beijing, China: Tsinghua UniversityPress-Springer Verlag, 2003.

[23]曹衛(wèi)星 , 李存東.小麥器官發(fā)育序列化命名方案[J].中國農(nóng)業(yè)科學(xué), 1997, 30(5):67-71.

CAO W, LI C. A Sequential naming scheme for devel‐opmental organs in wheat[J]. Scientia Agricultura Sini‐ca, 1997, 30(5):67-71.

[24]李英碩, 楊帆, 袁兆奎.空間圓形擬合檢測新方法[J].測繪科學(xué), 2013, 38(6):147-148.

LI Y, YANG F, YUAN Z. A new method of spatial cir‐cle fitting detection[J]. Science of Surveying and Map‐ping, 2013, 38(6):147-148.

[25]欒婉娜 , 劉成明.基于逆 Loop 細(xì)分的半正則網(wǎng)格簡化算法[J].圖學(xué)學(xué)報, 2020, 41(6):980-986.

LUAN W, LIU C. A semi-regular mesh simplificationalgorithm based on inverse Loop subdivision[J]. Jour‐ nal of Graphics, 2020, 41(6):980-986.

[26]諸葉平, 李世娟, 李書欽.作物生長過程模擬模型與形態(tài)三維可視化關(guān)鍵技術(shù)研究[J].智慧農(nóng)業(yè), 2019, 1(1):53-66.

ZHU Y, LI S, LI S. Research on key technologies of crop growth process simulation model and morphologi‐ cal 3D visualization[J]. Smart Agriculture, 2019, 1(1):53-66.

[27]吳倩, 孫颯爽, 趙哲民, 等.基于3DSOM 的植株三維重建方法研究[J].農(nóng)機(jī)化研究, 2017, 39(9):148-153.

WU Q, SUN Y, ZHAO Z, et al. Study on the 3D recon‐ struction method of plants based on 3DSOM[J]. Jour‐ nal of Agricultural Mechanization Research, 2017, 39(9):148-153.

[28]史維杰, 張吳平, 郝雅潔, 等.基于視覺三維重建的作物表型分析[J].湖北農(nóng)業(yè)科學(xué), 2019, 58(16):125-128.

SHI W, ZHANG W, HAO Y," et" al. Crop phenotypicanalysis" based" on" visual 3D" reconstruction[J]. HubeiAgricultural Sciences, 2019, 58(16):125-128.

[29]方偉, 馮慧, 楊萬能, 等.表型檢測中用于小麥株型研究的快速三維重建方法[J].中國農(nóng)業(yè)科技導(dǎo)報 ,2016, 18(2):95-101.

FANG W, FENG H, YANG W, et al. A fast 3D recon‐struction for wheat plant architecture studies in pheno‐typing[J]. Journal of Agricultural Science and Technol‐ogy, 2016, 18(2):95-101.

[30] YANG Y, ZHANG J, WU K, et al.3D point cloud onsemantic information for wheat reconstruction[J]. Agri‐culture, 2021, 11(5): ID 450.

Phenotypic Traits Extraction of Wheat Plants Using3D Digitization

ZHENG Chenxi1,2,3 , WEN Weiliang1,2* , LU Xianju1,2 , GUO Xinyu1,2 , ZHAO Chunjiang1,2,3

(1. Information Technology Research Center, Beijing Academy of Agriculture and Forestry Sciences, Beijing, 100097, China;2. Beijing Key Lab of Digital Plant, National Engineering Research Centerfor Information Technol‐ogy in Agriculture, Beijing, 100097, China;3. College of Information Engineering, Northwest Aamp;F University, Yan‐gling, 712100, China )

Abstract: Aiming at the difficulty of accurately extract the phenotypic traits of plants and organs from images or point clouds caused by the multiple tillers and serious cross-occlusion among organs of wheat plants, to meet the needs of accurate phenotyp‐ ic analysis of wheat plants, three-dimensional (3D) digitization was used to extract phenotypic parameters of wheat plants. First‐ ly, digital representation method of wheat organs was given and a 3D digital data acquisition standard suitable for the whole growth period of wheat was formulated. According to this standard, data acquisition was carried out using a 3D digitizer. Based on the definition of phenotypic parameters and semantic coordinates information contained in the 3D digitizing data, eleven con‐ ventional measurable phenotypic parameters in three categories were quantitative extracted, including lengths, thicknesses, and angles of wheat plants and organs. Furthermore, two types of new parameters for shoot architecture and 3D leaf shape were de‐ fined. Plant girth was defined to quantitatively describe the looseness or compactness by fitting 3D discrete coordinates based on the least square method. For leaf shape, wheat leaf curling and twisting were defined and quantified according to the direc‐ tion change of leaf surface normal vector. Three wheat cultivars including FK13, XN979, and JM44 at three stages (rising stage, jointing stage, and heading stage) were used for method validation. The Open3D library was used to process and visualize wheat plant data. Visualization results showed that the acquired 3D digitization data of maize plants were realistic, and the data acquisition approach was capable to present morphological differences among different cultivars and growth stages. Validation results showed that the errors of stem length, leaf length, stem thickness, stem and leaf angle were relatively small. The R2 were 0.93, 0.98, 0.93, and 0.85, respectively. The error of the leaf width and leaf inclination angle were also satisfactory, the R2 were 0.75 and 0.73. Because wheat leaves are narrow and easy to curl, and some of the leaves have a large degree of bending, the er‐ ror of leaf width and leaf angle were relatively larger than other parameters. The data acquisition procedure was rather time-con‐ suming, while the data processing was quite efficient. It took around 133 ms to extract all mentioned parameters for a wheat plant containing 7 tillers and total 27 leaves. The proposed method could achieve convenient and accurate extraction of wheat phenotypes at individual plant and organ levels, and provide technical support for wheat shoot architecture related research.

Key words: wheat; three-dimensional digitization; visualization; phenotypic traits extraction

主站蜘蛛池模板: 亚洲国产精品成人久久综合影院 | 老司机午夜精品网站在线观看| 欧美日韩中文国产va另类| 女人18一级毛片免费观看 | 欧美精品H在线播放| 在线精品欧美日韩| AV不卡国产在线观看| 久久综合伊人77777| 色吊丝av中文字幕| 韩日无码在线不卡| 中文国产成人精品久久| 伊人91在线| 精品久久久久久久久久久| 国产va欧美va在线观看| 这里只有精品免费视频| 国产成人综合亚洲欧美在| 国产在线无码av完整版在线观看| 人人澡人人爽欧美一区| 亚洲第一精品福利| 国产91无码福利在线| 亚洲中文字幕手机在线第一页| 一区二区三区四区精品视频| 亚洲AⅤ综合在线欧美一区| 2019年国产精品自拍不卡| 五月婷婷综合色| 国产无码高清视频不卡| www.狠狠| 日韩欧美中文| 成年网址网站在线观看| 国产迷奸在线看| 欧美日韩精品在线播放| 夜夜高潮夜夜爽国产伦精品| 成人中文在线| 黄色国产在线| 成人国产免费| 精品一区二区三区自慰喷水| 精品精品国产高清A毛片| 日本久久久久久免费网络| 精品国产自在在线在线观看| 国产精品综合色区在线观看| 国产女人在线视频| 久久青草精品一区二区三区| 九九热在线视频| www.91在线播放| 强奷白丝美女在线观看| 国产福利一区二区在线观看| 伊人色在线视频| 精品亚洲麻豆1区2区3区| 国产18在线播放| 中国国语毛片免费观看视频| 国产真实乱子伦精品视手机观看| 一级片免费网站| 青青青国产免费线在| 日韩中文字幕免费在线观看| 中国精品自拍| 精品无码一区二区在线观看| 亚洲系列无码专区偷窥无码| 亚洲an第二区国产精品| 伊人久久福利中文字幕| 午夜精品国产自在| www亚洲天堂| 丁香亚洲综合五月天婷婷| 91精品国产综合久久香蕉922| 又黄又湿又爽的视频| 99成人在线观看| 国产欧美精品一区aⅴ影院| AV色爱天堂网| AV无码无在线观看免费| 91精品综合| yy6080理论大片一级久久| 亚洲人成色在线观看| 91娇喘视频| 白浆视频在线观看| 国产99精品久久| 成人小视频在线观看免费| 毛片免费视频| 欧美成一级| 中国一级特黄视频| 乱色熟女综合一区二区| 性色生活片在线观看| 成人精品免费视频| 日本一区二区三区精品国产|