基于肌肉的人體運(yùn)動仿真研究

馬文超，侯 進(jìn)，曹寧哲

(西南交通大學(xué) 信息科學(xué)與技術(shù)學(xué)院，四川 成都 611756 )

基于肌肉的人體運(yùn)動仿真研究

馬文超，侯 進(jìn)，曹寧哲

(西南交通大學(xué) 信息科學(xué)與技術(shù)學(xué)院，四川 成都 611756 )

針對當(dāng)前國內(nèi)缺少滿足牛頓運(yùn)動定律的人體運(yùn)動仿真問題，提出了基于肌肉的人體運(yùn)動仿真方法。該方法首先建立了基于解剖學(xué)的人體骨骼肌肉模型，然后進(jìn)行運(yùn)動捕獲數(shù)據(jù)與人體模型的匹配，其次構(gòu)建了主要關(guān)節(jié)點(diǎn)的局部坐標(biāo)系，接著運(yùn)用多剛體運(yùn)動學(xué)與動力學(xué)進(jìn)行運(yùn)動學(xué)及動力學(xué)參數(shù)的計算，之后施加空間力學(xué)平衡約束得到了有效的關(guān)節(jié)運(yùn)動約束機(jī)制。實(shí)驗(yàn)結(jié)果表明，利用該算法實(shí)現(xiàn)的人體運(yùn)動仿真真實(shí)、自然、逼真，證明了算法的有效性。關(guān)鍵詞: 肌肉；骨骼肌肉模型；運(yùn)動捕獲；多剛體運(yùn)動學(xué)

人體運(yùn)動仿真是生物力學(xué)、機(jī)器人學(xué)和計算機(jī)科學(xué)交叉產(chǎn)生的研究領(lǐng)域，是一種利用計算機(jī)模擬真實(shí)人體運(yùn)動過程的技術(shù)，具體包括建立數(shù)學(xué)模型，計算給定虛擬人約束條件下自然真實(shí)的物理運(yùn)動過程，并在利用計算機(jī)生成的虛擬環(huán)境中以三維圖像方式逼真呈現(xiàn)該運(yùn)動過程[1]。與人體骨骼動畫相比，該技術(shù)實(shí)現(xiàn)的人體運(yùn)動符合牛頓運(yùn)動定律，既可增強(qiáng)虛擬環(huán)境的真實(shí)性與沉浸感，還可進(jìn)行現(xiàn)實(shí)任務(wù)的完成和操作培訓(xùn)等。

計算機(jī)技術(shù)與虛擬人應(yīng)用的發(fā)展，人們對符合人體生理特性以及運(yùn)動特征的虛擬人模型需求越來越高。Chadwick等人[2]重點(diǎn)對人體手臂和肩部的運(yùn)動特性進(jìn)行了深入分析，為了減少算法運(yùn)算時間，達(dá)到模型對控制參數(shù)改變進(jìn)行實(shí)時反應(yīng)的效果，采用多項(xiàng)式來表示肌肉的行為狀態(tài)，取得了很好的效果。Ayusawa等人[3]利用非線性規(guī)劃，通過遞推牛頓-歐拉方程進(jìn)行梯度計算，解決了傳統(tǒng)方法需要在給定時間序列情況下進(jìn)行大規(guī)模的逆運(yùn)動學(xué)計算的問題，從而可以快速從運(yùn)動捕獲數(shù)據(jù)中獲得人體運(yùn)動學(xué)參數(shù)。Jung等人[4]分別在僅提供關(guān)節(jié)運(yùn)動學(xué)參數(shù)以及同時提供腳部壓力的情況下，對反作用力模型進(jìn)行了深入探討，并在ANYBODY中進(jìn)行了驗(yàn)證。對在某些不適合放置測力墊的情景下進(jìn)行地面作用力的分析，提供了極大的便利。Xu等人[5]對肌肉系統(tǒng)進(jìn)行了研究，從肌肉的生理形狀出發(fā)，提出虛擬肌肉控制模型，將研究成果應(yīng)用到仿生領(lǐng)域，對機(jī)器人手臂的控制取得了良好效果。Dustin等人[6]針對目前控制算法缺乏魯棒性等缺點(diǎn)，建立了新的人體手部的肌肉骨骼模型，通過肌電信號來預(yù)測手和腕關(guān)節(jié)的運(yùn)動。戚福州[7]利用國際人體測量學(xué)數(shù)據(jù)及層級建模的方法，構(gòu)造了具有骨骼層與皮膚層的雙層虛擬人模型，實(shí)現(xiàn)了虛擬人角色建模。黃永坤[8]針對虛擬人運(yùn)動過程中皮膚出現(xiàn)塌陷失真的問題，運(yùn)用構(gòu)造二次貝塞爾曲線的方法，使該問題得到明顯改善。Assassi[9]利用人體解剖學(xué)數(shù)據(jù)結(jié)合人體運(yùn)動學(xué)及動力學(xué)原理，對髖關(guān)節(jié)軟骨運(yùn)動進(jìn)行了深入分析。Stanev[10]使用kinect捕獲人體的關(guān)節(jié)運(yùn)動軌跡，其次運(yùn)用逆向運(yùn)動學(xué)及動力學(xué)對人體各關(guān)節(jié)運(yùn)動參數(shù)進(jìn)行了求解，然后采用靜態(tài)優(yōu)化算法對肌肉力進(jìn)行了估算。

目前國內(nèi)外眾多組織和研究者都把工作重點(diǎn)放到各自應(yīng)用領(lǐng)域的虛擬人研究，只對虛擬人特定部位(上肢、頸部等)進(jìn)行仿真評價，缺乏對虛擬人在滿足牛頓運(yùn)動定律下的整體運(yùn)動仿真的研究。本文首先通過分析真實(shí)人體結(jié)構(gòu)特征，建立了基于人體解剖學(xué)和人體測量學(xué)的參數(shù)化全身人體骨肌模型。然后對運(yùn)動捕獲數(shù)據(jù)進(jìn)行整理，得到實(shí)驗(yàn)對象的運(yùn)動學(xué)模型，并將這些運(yùn)動學(xué)參數(shù)與所建立人體模型進(jìn)行匹配。其次建立主要關(guān)節(jié)點(diǎn)的局部坐標(biāo)系。接著將各體段視為剛體，具有質(zhì)量、質(zhì)心和轉(zhuǎn)動慣量，進(jìn)而運(yùn)用多剛體運(yùn)動學(xué)與動力學(xué)進(jìn)行質(zhì)心加速度、角加速度、力與力矩的計算，之后利用空間力學(xué)平衡約束，得到了有效的關(guān)節(jié)運(yùn)動約束機(jī)制，實(shí)現(xiàn)了人體運(yùn)動仿真。系統(tǒng)流程如圖1所示。

圖1 人體運(yùn)動仿真系統(tǒng)整體流程圖

1 人體幾何模型的建立

1.1 人體骨骼系統(tǒng)構(gòu)建

人體骨骼系統(tǒng)十分復(fù)雜，包含206塊骨骼，約占人體體重的20%，主要分為軀干骨、顱骨、上肢骨以及下肢骨。因此，在保證人體結(jié)構(gòu)真實(shí)性的情況下，為了降低仿真的復(fù)雜度，需要對人體結(jié)構(gòu)進(jìn)行必要的簡化。目前，H-Anim是常用的表示虛擬人的國際標(biāo)準(zhǔn)。通過對人體關(guān)節(jié)的分析可得，各個關(guān)節(jié)符合層級結(jié)構(gòu)特征。各個節(jié)點(diǎn)采用繼承機(jī)構(gòu)，形成兄弟節(jié)點(diǎn)和兒子節(jié)點(diǎn)。一般根據(jù)實(shí)際需要，只保留能體現(xiàn)人體運(yùn)動特征的關(guān)節(jié)點(diǎn)。如圖2所示。本文所使用的人體骨骼模型數(shù)據(jù)均由三維圖形軟件獲取，并對骨骼模型數(shù)據(jù)進(jìn)行處理，獲得了簡化的人體骨骼模型，按人體骨性標(biāo)志將人體分割成若干段，每段稱為體段。以骨性標(biāo)志為分界點(diǎn)對體段進(jìn)行劃分，得到一個15個體段，29個自由度(Degree of Freedom， DOF)的人體解剖骨骼模型。如圖3所示。

圖2 人體關(guān)節(jié)樹形結(jié)構(gòu)圖

圖3 人體解剖模型

1.2 人體關(guān)節(jié)模型

各骨之間借助結(jié)締組織相互連接，形成骨連結(jié)。根據(jù)連結(jié)方式，分為直接連結(jié)和間接連結(jié)。間接連結(jié)即關(guān)節(jié)。人體的關(guān)節(jié)是骨骼的骨塊(如股骨和髖骨)之間的連接點(diǎn)，也是伸展力、扭轉(zhuǎn)力的傳輸以及各類運(yùn)動所必須的結(jié)構(gòu)。筆者按照關(guān)節(jié)實(shí)際運(yùn)動中的自由度對人體的關(guān)節(jié)進(jìn)行分類同時確定活動范圍：不能旋轉(zhuǎn)，只能上下彎曲的1個DOF的關(guān)節(jié)，如膝關(guān)節(jié)。既能旋轉(zhuǎn)又能擺動的3個DOF的關(guān)節(jié)，如腰椎關(guān)節(jié)。在旋轉(zhuǎn)和平移上均有3個DOF的關(guān)節(jié)，如根關(guān)節(jié)。各關(guān)節(jié)自由度及活動范圍分配情況如表1所示。

表1 各關(guān)節(jié)自由度及活動范圍

關(guān)節(jié)名稱自由度XYZmaxminmaxminmaxmin關(guān)節(jié)名稱自由度XYZmaxminmaxminmaxmin左肩關(guān)節(jié)3180-18090-45135-45右肩關(guān)節(jié)3180-18045-90135-45左髖關(guān)節(jié)325-12560-3030-60右髖關(guān)節(jié)3125-2530-6030-60左肘關(guān)節(jié)15-145180-1355-145右肘關(guān)節(jié)1145-5135-1805-145左腕關(guān)節(jié)125-50110-9090-80右腕關(guān)節(jié)150-2590-11090-80左膝關(guān)節(jié)10040-301200右膝關(guān)節(jié)10030-401200腰椎關(guān)節(jié)640-4070-70100-30根關(guān)節(jié)645-4550-50160-10

1.3 肌肉幾何建模的實(shí)現(xiàn)

目前肌肉的幾何建模主要是用肌肉力作用線來表示。第1種為直線段(straight-line)法。利用肌肉在骨骼上的附著點(diǎn)之間的連線進(jìn)行表示，此方法適合肌肉附著點(diǎn)單一的肌肉。第2種為質(zhì)心線段(centroid-line)法。在充分考慮肌纖維方向的基礎(chǔ)上，假定肌肉橫截面上的質(zhì)心構(gòu)成了作用在肌肉上的力的曲線[11]。本文根據(jù)所需肌肉的不同種類，分別采用這兩種肌肉的幾何建模方法，共建立92條肌肉。如圖4所示。

圖4 人體骨肌模型

2 運(yùn)動數(shù)據(jù)的處理

2.1 運(yùn)動數(shù)據(jù)的獲取

近年來在人體運(yùn)動仿真領(lǐng)域中，常采用運(yùn)動捕獲設(shè)備來獲取運(yùn)動數(shù)據(jù)。運(yùn)動捕獲設(shè)備通過使用位置傳感器測量到物體的運(yùn)動軌跡，并將運(yùn)動數(shù)據(jù)保存下來。然而，目前具有較高精度的運(yùn)動捕獲設(shè)備均為國外生產(chǎn)，價格昂貴?？突仿〈髮W(xué)提供了免費(fèi)的運(yùn)動捕獲數(shù)據(jù)，本文所使用的運(yùn)動數(shù)據(jù)就來自此數(shù)據(jù)庫。

2.2 骨肌系統(tǒng)的匹配

為了將運(yùn)動數(shù)據(jù)作用到本系統(tǒng)所建立的人體骨肌模型上，需要完成被測對象骨肌系統(tǒng)與本系統(tǒng)人體模型的匹配。匹配因子的計算方法如下：定義兩個相鄰的標(biāo)記點(diǎn)(experimental marker)為一組，它們的距離為e1，骨肌系統(tǒng)對應(yīng)的一組虛擬標(biāo)記點(diǎn)(virtual marker)的距離為m1，則定義此組比例為

(1)

其中：e1由運(yùn)動捕獲數(shù)據(jù)獲得，m1由人體骨肌模型獲得。取比例因子為所有相鄰標(biāo)記點(diǎn)比例的均值，即

(2)

其中：n為所有兩個相鄰標(biāo)記點(diǎn)構(gòu)成的組數(shù)。

2.3 關(guān)節(jié)坐標(biāo)系的構(gòu)建

本文按照人體解剖結(jié)構(gòu)，分別在主要關(guān)節(jié)點(diǎn)建立了局部坐標(biāo)系。其中定義：同時與地面垂直與身體長軸平行的為垂直軸。與垂直軸垂直同時與地面平行，方向?yàn)樯眢w前后方向的軸為矢狀軸。與身體長軸垂直同時與地面平行，身體左右方向的軸為冠狀軸。矢狀面是沿矢狀軸方向?qū)⑷梭w分為左、右兩部分的面。冠狀面是沿冠狀軸方向?qū)⑷梭w分為前、后兩部分的面。水平面是同時垂直于矢狀面與冠狀面的面，將人體分為上、下兩部分。

(3)

則完成從X,Y,Z到Xi,Yi,Zi的坐標(biāo)變換的矩陣為T·R，坐標(biāo)變換公式為

(xi,yi,zi,1)=(xi,yi,zi,1)·Ti·Ri

(4)

經(jīng)推導(dǎo)化簡得從Xi+1Yi+1Zi+1到下一時刻坐標(biāo)Xi+1Yi+1Zi+1的平移矩陣為(Ti+2)-1·Ti，旋轉(zhuǎn)矩陣為Ri·(Ri+1)-1。

2.4 多剛體運(yùn)動學(xué)參數(shù)計算

運(yùn)動學(xué)參數(shù)的計算包括瞬時速度、加速度、瞬時角速度和角加速度。若要求解剛體的瞬時角速度以及角加速度，需要將旋轉(zhuǎn)矩陣用剛體旋轉(zhuǎn)角α,β,γ表示，需完成旋轉(zhuǎn)變換操作，方法如下：

(5)

整體旋轉(zhuǎn)矩陣R為3個坐標(biāo)軸的基本旋轉(zhuǎn)矩陣的乘積，即

(6)

則取任意坐標(biāo)(xi,yi,zi)及其對應(yīng)的坐標(biāo)(x,y,z)有以下轉(zhuǎn)換關(guān)系

(xi,yi,zi,1)=(x,y,z,1)·R

(7)

將式(7)與式(6)聯(lián)立求解可得在任意坐標(biāo)下剛體的旋轉(zhuǎn)角α,β,γ的值。

該剛體質(zhì)心的瞬時位移可由式(3)獲得，即

(8)

進(jìn)而可得剛體的瞬時速度vi,vi+1及加速度ai為

(9)

剛體角加速度可根據(jù)所求得旋轉(zhuǎn)角α,β,γ的值進(jìn)行求解，即

(10)

式中：i表示運(yùn)動數(shù)據(jù)的幀數(shù)；ω為質(zhì)心瞬時角速度。

2.5 多剛體動力學(xué)參數(shù)計算

則根據(jù)牛頓第二定律及牛頓-歐拉平衡方程可得矢量關(guān)系如下

(11)

(12)

假設(shè)力的作用點(diǎn)Q在世界坐標(biāo)系下的坐標(biāo)值為X,Y,Z，則對于任一力對軸之矩的解析式為

Mx(F)=YFz-ZFy,My(F)=ZFx-XFz,

Mz(F)=XFy-YFx

(13)

式中：Fx,Fy,Fz分別為力F在X,Y,Z軸上的分量；Mx(F),My(F),Mz(F)分別為力矩M(F)在X,Y,Z軸上的分量。

通過聯(lián)立求解代入所有參量后的式(11)～(13)，可得

(14)

(15)

(16)

(17)

2.6 多剛體系統(tǒng)平衡約束

則根據(jù)牛頓第二定律及牛頓-歐拉平衡方程可得矢量關(guān)系

(18)

(19)

體段1的空間力系平衡方程可由式(18)～(19)可獲得，即

(20)

(21)

(22)

(23)

(24)

(25)

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

本次實(shí)驗(yàn)數(shù)據(jù)的采集對象為一名健康男性，年齡38歲，身高1.83 m，體重65.9 kg，速度為1.96 m/s。采樣頻率為60 Hz，共采集500幀數(shù)據(jù)。由于運(yùn)動捕獲數(shù)據(jù)存在噪聲干擾，因此在實(shí)驗(yàn)仿真前對數(shù)據(jù)進(jìn)行了低通濾波，截止頻率為50 Hz。依據(jù)濾波后的運(yùn)動數(shù)據(jù)，本文實(shí)現(xiàn)了運(yùn)動仿真，并選取一個步態(tài)周期進(jìn)行展示。如圖5和圖6所示，該動作始于左腳與地面接觸，終于右腳與地面接觸，動作完成耗時0.912 s。圖5分別為第0幀，第5幀和第10幀。圖6分別為第15幀、第20幀和第25幀。

圖5 第0、5、10幀運(yùn)動圖像

圖6 第15、20、25幀運(yùn)動圖像

從圖中可看出，第0幀時，人體運(yùn)動狀態(tài)發(fā)生較大改變使加速度增大，從而使得右小腿骨及膝關(guān)節(jié)承受的力和扭矩比較大。第25幀時，在人體運(yùn)動狀態(tài)改變結(jié)束的情況下，此時右腿的加速度減小，則右小腿骨及膝關(guān)節(jié)承受的力和扭矩也相應(yīng)的較小。圖7分別為第2幀、第12幀和第22幀。圖8分別為第28幀、第38幀和第48幀，分別展示右腿股內(nèi)側(cè)肌以及左腿股內(nèi)側(cè)肌在運(yùn)動過程中的變化情況。由圖可見，肌肉收縮正常，肌肉作用線隨著膝關(guān)節(jié)角度的增大而加長且經(jīng)圖7和圖8對比可得，左右肌肉在相同關(guān)節(jié)角度下肌肉作業(yè)線長短相同。

圖7 右側(cè)大腿股內(nèi)側(cè)肌變化情況

圖8 左側(cè)大腿股內(nèi)側(cè)肌變化情況

結(jié)果與文獻(xiàn)[6]和[9]相比，本文實(shí)現(xiàn)了人體整體運(yùn)動的仿真，而文獻(xiàn)[6]和[9]只是對人體特定部位進(jìn)行了局部的仿真分析。文獻(xiàn)[7]則只是完成了具有骨骼層與皮膚層雙層結(jié)構(gòu)的人體建模，并實(shí)現(xiàn)人體的運(yùn)動仿真。與文獻(xiàn)[8]的仿真結(jié)果對比可得，本文實(shí)現(xiàn)的虛擬人的仿真更加真實(shí)，真實(shí)地再現(xiàn)了人體結(jié)構(gòu)，虛擬人動作更加逼真，并且增加了肌肉層面的仿真與分析，運(yùn)動更加符合人類特性。而且各個關(guān)節(jié)約束更加合理，動作更加逼真。文獻(xiàn)[8]則僅僅是完成了作者的指定動作，不能實(shí)現(xiàn)運(yùn)動捕獲數(shù)據(jù)的驅(qū)動，而且動作并不符合牛頓運(yùn)動定律。由仿真結(jié)果可觀察到，虛擬標(biāo)記點(diǎn)與實(shí)驗(yàn)數(shù)據(jù)采集點(diǎn)在空間上的位置基本一致，很好地還原了實(shí)驗(yàn)對象的真實(shí)運(yùn)動。各體段約束協(xié)調(diào)，關(guān)節(jié)角度合理，未出現(xiàn)失真，整個運(yùn)動過程真實(shí)、自然、逼真，符合人體生理特性及運(yùn)動規(guī)律。

4 小結(jié)

本文采用上面提到的人體骨骼和肌肉模型的建模方法以及運(yùn)動捕獲數(shù)據(jù)的處理方法，實(shí)現(xiàn)了基于肌肉的人體運(yùn)動仿真，同時滿足牛頓運(yùn)動定律且符合人體運(yùn)動特征。仿真效果真實(shí)、自然逼真。實(shí)驗(yàn)結(jié)果表明，該方法可靠、有效。但是，由于數(shù)據(jù)的計算量較大，對計算機(jī)硬件有較高要求。下一步筆者的工作重點(diǎn)是進(jìn)一步優(yōu)化數(shù)據(jù)處理算法以及對如何將人體模型運(yùn)用于特定場景下進(jìn)行各關(guān)節(jié)和肌肉力的分析展開研究。

[1] 夏時洪，魏毅，王兆其.人體運(yùn)動仿真綜述[J].計算機(jī)研究與發(fā)展，2010，47(4)：1354-1361.

[2] CHADWICK E K， BLANA D， KIRSCH R F. Real-time simulation of three-dimenssional shoulder girdle and arm dynamics[J]. IEEE transactions on biomedical engineering， 2014， 61(7)：1947-1956.

[3] AYUSAWA K， IKEGAMI Y， NAKAMURA Y. Simultaneous global inverse kinematics and geometric parameter identification of human skeletal model from motion capture data[J]. Mechanism and machine theory， 2014， 74(4)：274-284.

[4] JUNG Y， JUNG M， LEE K， et al. Ground reaction force estimation using an insole-type pressure mat and joint kinematics during walking[J]. Journal of biomechanis， 2014， 47(11)：2693-2699.

[5] XU D， ZHANG S G， WEI H X. Virtual musculoskeletal control model with a spindle-like fuzzy algorithm for robotic compliance[J].Applied mathematical modelling， 2015， 39(12)：3265-3279.

[6] CROUCH D L， HUANG H. Musculoskeletal model predicts multi-joint wrist and hand movement from limited EMG control signals[C]//Proc. 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. Milan：[s.n.]，2015：1132-1135.

[7] 戚福州，侯進(jìn)，黃永坤.一種基于國標(biāo)人體測量學(xué)數(shù)據(jù)的虛擬人建模方法[J].計算機(jī)工程與科學(xué)， 2015， 37(4)：783-789.

[8] 黃永坤，侯進(jìn)，戚福州.三維虛擬人皮膚變形與塌陷部位的調(diào)整[J].計算機(jī)應(yīng)用與軟件，2015，32(9)：184-187.

[9] ASSASSI L，THALMANN N M. Biomechanical modeling of the human articulation for osteoarthritis analysis[J]. Computer methodsin biomechanics and biomedical engineering， 2012，28(1)： 573-578.

[10] STAVEN D，MOUSTAKAS K. Virtual human behavioural profile extraction using Kinect based motion tracking[C] //International Conference on Cyberworlds. Santander：[s.n.]，2014：411-414.

[11] TANG G， WANG C T. A muscle-path-plane method for representing muscle contraction during joint movement[J]. Computer methods in biomechanics and biomedical engineering，2010， 14(1)： 59-69.

馬文超(1989— )，碩士生，主研虛擬現(xiàn)實(shí)；

侯 進(jìn)(1969— )，副教授，主要研究方向?yàn)槿藱C(jī)交互、虛擬現(xiàn)實(shí)、人工智能和自然語言處理；

曹寧哲(1990— )，碩士生，主研虛擬現(xiàn)實(shí)。

責(zé)任編輯：閆雯雯

Study of simulations of virtual human movement based on muscle

MA Wenchao，HOU Jin，CAO Ningzhe

(SchoolofInformationScienceandTechnology，SouthwestJiaotongUniversity，Chengdu611756，China)

Aiming at the questions that the current domestic is lacking in the human motion simulations which obey Newton’s laws of motion， a method of simulations of human movement based on muscle is proposed. The method firstly establishes a musculoskeletal model based on anatomy； then， processes the motion capture data to match the musculoskeletal model； next， calculates the kinematic and kinetic parameters by multi-body kinematics and dynamics and applies the mechanical balance of space constraints in order to obtain the effective joint motion constraint mechanism. Experimental results show that the simulation of human motion is real， natural and realistic by using the proposed algorithm， and demonstrates the effectiveness of the algorithm.

muscle； musculoskeletal model； motion capture； multi-body kinematic

馬文超，侯進(jìn)，曹寧哲. 基于肌肉的人體運(yùn)動仿真研究[J].電視技術(shù)，2017，41(2)：98-103. MA W C， HOU J， CAO N Z. Study of simulations of virtual human movement based on muscle[J]. Video engineering，2017，41(2)：98-103.

TN919.3

A

10.16280/j.videoe.2017.02.020

國家自然科學(xué)基金項(xiàng)目(61371165)；四川省動漫研究中心2015年度科研項(xiàng)目(DM201504)；成都市科技項(xiàng)目(科技惠民技術(shù)研發(fā)項(xiàng)目)(2015-HM01-00050-SF)

2016-05-24