基于無(wú)監(jiān)督學(xué)習(xí)的三維肺部CT圖像配準(zhǔn)方法研究

姜?杉，張紅運(yùn)，楊志永，張國(guó)彬

基于無(wú)監(jiān)督學(xué)習(xí)的三維肺部CT圖像配準(zhǔn)方法研究

姜?杉，張紅運(yùn)，楊志永，張國(guó)彬

(天津大學(xué)機(jī)械工程學(xué)院，天津 300350)

三維肺部電子計(jì)算機(jī)斷層掃描(computed tomography，CT)圖像非剛性配準(zhǔn)是醫(yī)學(xué)圖像配準(zhǔn)領(lǐng)域中最重要的任務(wù)之一．但是，肺部組織受呼吸運(yùn)動(dòng)影響而產(chǎn)生的非線性形變與大尺度位移給三維肺部CT圖像的非剛性配準(zhǔn)帶來(lái)巨大挑戰(zhàn)．針對(duì)這一難題，設(shè)計(jì)開(kāi)發(fā)了一種基于無(wú)監(jiān)督學(xué)習(xí)端到端的配準(zhǔn)方法．通過(guò)改進(jìn)現(xiàn)有U-Net神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，在跳接之間引入Inception模塊，充分融合多尺度深層特征生成高精度的稠密位移向量場(chǎng)．為保證位移向量場(chǎng)光滑，在損失函數(shù)中加入雅可比正則化項(xiàng)，以達(dá)到訓(xùn)練中顯式懲罰位移向量場(chǎng)中奇點(diǎn)的目的．另外，為緩解現(xiàn)有公開(kāi)數(shù)據(jù)資源有限導(dǎo)致的過(guò)擬合問(wèn)題，提出了一種基于三維薄板樣條(3D-thin plate spline，3D-TPS)變換的數(shù)據(jù)增強(qiáng)方法實(shí)現(xiàn)對(duì)訓(xùn)練數(shù)據(jù)的擴(kuò)充，將具有60套三維肺部CT圖像的訓(xùn)練數(shù)據(jù)集EMPIRE10擴(kuò)充為6060套以滿足卷積神經(jīng)網(wǎng)絡(luò)訓(xùn)練的需要．設(shè)計(jì)驗(yàn)證實(shí)驗(yàn)，通過(guò)與基于學(xué)習(xí)的Voxelmorph方法和兩個(gè)包含傳統(tǒng)方法配準(zhǔn)工具包ANTs和Elastix進(jìn)行比較．實(shí)驗(yàn)結(jié)果表明：在公開(kāi)可用的DIR-Lab 4DCT數(shù)據(jù)集上，所提出的方法在目標(biāo)配準(zhǔn)誤差(target registration error，TRE)上達(dá)到次優(yōu)的2.09mm，平均Dice得分達(dá)到最優(yōu)的0.987，同時(shí)所生成的扭曲圖像中幾乎不存在折疊體素．

非剛性配準(zhǔn)；卷積神經(jīng)網(wǎng)絡(luò)；無(wú)監(jiān)督學(xué)習(xí)；數(shù)據(jù)增強(qiáng)；雅可比正則化

三維肺部CT圖像已廣泛應(yīng)用于圖像引導(dǎo)下的放射治療[1]．肺作為典型的運(yùn)動(dòng)器官，在面向肺癌的放射治療之前，為估計(jì)所需的精確劑量必須追蹤肺呼吸運(yùn)動(dòng)[2]．術(shù)中需要將不同呼吸狀態(tài)的三維肺部CT圖像與參考狀態(tài)的三維肺部CT圖像進(jìn)行對(duì)齊，跟蹤每個(gè)體素的劑量．三維肺部CT圖像非剛性配準(zhǔn)是準(zhǔn)確跟蹤肺部運(yùn)動(dòng)的一種可行方法．

非剛性配準(zhǔn)是指在一對(duì)待配準(zhǔn)圖像間建立非線性稠密映射關(guān)系．現(xiàn)有的傳統(tǒng)算法大多通過(guò)幾何方法解決每個(gè)圖像對(duì)的非剛性配準(zhǔn)優(yōu)化問(wèn)題，如SyN (standard symmetric normalization)[3]和Diffeomorphic Demons[4]．傳統(tǒng)方法計(jì)算開(kāi)銷很大，每當(dāng)配準(zhǔn)一對(duì)新圖像時(shí)，需重新迭代運(yùn)算．近年來(lái)，越來(lái)越多的學(xué)者傾向于利用基于學(xué)習(xí)方法完成配準(zhǔn)任務(wù)．訓(xùn)練后的卷積神經(jīng)網(wǎng)絡(luò)模型可以在幾秒鐘內(nèi)完成一對(duì)三維醫(yī)學(xué)圖像的配準(zhǔn)，且其精度足以媲美較成熟的傳統(tǒng)方法．深度學(xué)習(xí)下的圖像配準(zhǔn)大致可分為兩類：基于監(jiān)督學(xué)習(xí)的配準(zhǔn)方法和基于無(wú)監(jiān)督學(xué)習(xí)的配準(zhǔn)方法．

基于監(jiān)督學(xué)習(xí)的配準(zhǔn)方法在訓(xùn)練網(wǎng)絡(luò)時(shí)，需要提供與訓(xùn)練樣本相對(duì)應(yīng)的標(biāo)簽(ground truth，GT)．Fan等[5]使用傳統(tǒng)方法對(duì)圖像進(jìn)行配準(zhǔn)，將獲取到位移向量場(chǎng)作為GT，然而此方法獲得的GT與真實(shí)GT有所偏差．Hu等[6]提出用分割的解剖結(jié)構(gòu)掩膜作為GT來(lái)訓(xùn)練卷積神經(jīng)網(wǎng)絡(luò)．在該方法中，卷積神經(jīng)網(wǎng)絡(luò)以固定和運(yùn)動(dòng)圖像對(duì)作為輸入，通過(guò)對(duì)解剖結(jié)構(gòu)掩膜進(jìn)行對(duì)齊達(dá)到配準(zhǔn)圖像的目的．Miao等[7]將原始圖像作為浮動(dòng)圖像，將被模擬位移向量場(chǎng)扭曲的原始圖像作為固定圖像，將模擬位移向量場(chǎng)作為GT．盡管監(jiān)督學(xué)習(xí)在圖像配準(zhǔn)領(lǐng)域有很大的潛力，但是通過(guò)傳統(tǒng)的配準(zhǔn)工具獲取GT過(guò)程較為繁瑣．

基于無(wú)監(jiān)督學(xué)習(xí)的配準(zhǔn)無(wú)需額外手工生成的GT，在訓(xùn)練過(guò)程中利用固定圖像和扭曲圖像之前的差異來(lái)指導(dǎo)網(wǎng)絡(luò)參數(shù)的優(yōu)化．VoxelMorph為一種典型的基于無(wú)監(jiān)督學(xué)習(xí)的三維醫(yī)學(xué)圖像端到端配準(zhǔn)方法，將配準(zhǔn)中發(fā)生的扭曲變換定義為參數(shù)函數(shù)，并對(duì)其參數(shù)進(jìn)行了優(yōu)化[8]．Zhao等[9]采用級(jí)聯(lián)配準(zhǔn)子網(wǎng)絡(luò)預(yù)測(cè)位移向量場(chǎng)，取得了較好的效果．上述方法在腦、肝圖像上取得了良好的效果，而神經(jīng)網(wǎng)絡(luò)在肺部CT圖像配準(zhǔn)中應(yīng)用較少．相比較于腦、肝圖像，不同時(shí)刻采集到的肺部圖像往往會(huì)由于人體的呼吸運(yùn)動(dòng)產(chǎn)生大范圍復(fù)雜非線性形變和位移，這將導(dǎo)致配準(zhǔn)肺部圖像較為困難．

本文提出一種基于無(wú)監(jiān)督學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)方法用于三維肺部CT圖像非剛性配準(zhǔn)．在現(xiàn)有U-Net神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)上進(jìn)行改進(jìn)，將Inception模塊引入至跳接之間以提升網(wǎng)絡(luò)的特征檢測(cè)能力，充分融合多尺度深層特征生成高精度的稠密位移向量場(chǎng)．為了抑制不可逆變形，對(duì)稠密位移向量場(chǎng)中的具有負(fù)雅可比行列式值的變換進(jìn)行懲罰．實(shí)驗(yàn)表明，在變形后的圖像中幾乎沒(méi)有折疊體素．訓(xùn)練后的卷積神經(jīng)網(wǎng)絡(luò)模型，可實(shí)現(xiàn)端到端的三維肺部CT圖像非剛性配準(zhǔn)，同時(shí)具有較強(qiáng)的泛化能力．

1?方?法

圖1 基于無(wú)監(jiān)督學(xué)習(xí)的三維肺部CT圖像配準(zhǔn)流程

Fig.1 Flow chart of 3D lung CT image registration based on unsupervised learning

1.1?卷積網(wǎng)絡(luò)結(jié)構(gòu)

如圖2所示，使用的卷積網(wǎng)絡(luò)是基于U-Net的改進(jìn)．通常情況下，為生成高精度的稠密位移向量場(chǎng)，可在U-Net上堆疊卷積層加深網(wǎng)絡(luò)來(lái)實(shí)現(xiàn)．但是隨著網(wǎng)絡(luò)深度的增加，可能會(huì)出現(xiàn)梯度消失現(xiàn)象，進(jìn)而導(dǎo)致模型過(guò)擬合．另外，不同個(gè)體的肺部CT圖像中信息位置存在較大差異，因此選取合適尺寸的卷積核變得十分困難．針對(duì)上述問(wèn)題，本文不再通過(guò)盲目加深U-Net來(lái)優(yōu)化網(wǎng)絡(luò)性能，而是通過(guò)引入Inception模塊加寬網(wǎng)絡(luò)，在保證U-Net網(wǎng)絡(luò)一定深度的基礎(chǔ)上，通過(guò)Inception模塊所特有的多尺度卷積操作，促進(jìn)網(wǎng)絡(luò)捕捉并學(xué)習(xí)更加豐富的多尺度特征．配合U-Net中跳躍連接所融合的上、下文信息，以達(dá)到生成高精度稠密位移向量場(chǎng)的目的．

圖2?3D卷積網(wǎng)絡(luò)模型結(jié)構(gòu)

1.2?損失函數(shù)

損失函數(shù)整體可以表示為

2?實(shí)?驗(yàn)

2.1?數(shù)據(jù)集

本文使用3個(gè)包含同一病人不同時(shí)刻的三維肺部圖像數(shù)據(jù)集：EMPIRE10[12]、DIR-Lab 4DCT[13]和POPI[14]．DIR-Lab 4DCT數(shù)據(jù)集包含10個(gè)四維肺部CT圖像，每個(gè)四維肺部CT圖像中包含一個(gè)完整呼吸周期中的10個(gè)呼吸相位．選擇具有最大相對(duì)變形的兩個(gè)階段數(shù)據(jù)：呼氣末時(shí)刻和吸氣末時(shí)刻，兩階段分別標(biāo)注有300個(gè)地標(biāo)點(diǎn)．EMPIRE10數(shù)據(jù)集由30對(duì)三維肺部CT組成．POPI數(shù)據(jù)集由6個(gè)四維肺部CT圖像組成，使用最大相對(duì)變形的兩個(gè)階段：呼氣末時(shí)刻和吸氣末時(shí)刻，兩階段分別標(biāo)注有100個(gè)地標(biāo)點(diǎn).

2.2?數(shù)據(jù)增強(qiáng)

3D-TPS插值函數(shù)的參數(shù)能夠通過(guò)求解下面線性方程組得出．

(9)

圖3?3D-TPS生成圖像示例

2.3?評(píng)價(jià)指標(biāo)

由于CT圖像中肺實(shí)質(zhì)區(qū)域往往會(huì)占據(jù)大量的體素，因此使用地標(biāo)點(diǎn)間的TRE來(lái)評(píng)估肺內(nèi)部的配準(zhǔn)質(zhì)量．

Dice雖然不能很好地評(píng)價(jià)肺實(shí)質(zhì)區(qū)域的配準(zhǔn)質(zhì)量，但可用于評(píng)估肺邊界的對(duì)齊程度．Dice分?jǐn)?shù)表達(dá)式為

此外，圖像折疊在解剖學(xué)上是不合理的．第1.2節(jié)中提到的折疊體素的數(shù)量也作為算法的評(píng)價(jià)指標(biāo)之一．

2.4?對(duì)比實(shí)驗(yàn)

首先，將所提出的算法與ANTs[15]配準(zhǔn)工具包中提供的SyN算法進(jìn)行了比較．實(shí)驗(yàn)表明，使用SyN默認(rèn)參數(shù)對(duì)肺部CT圖像進(jìn)行配準(zhǔn)，達(dá)不到最佳配準(zhǔn)性能．通過(guò)在多個(gè)數(shù)據(jù)上開(kāi)展大量實(shí)驗(yàn)，可得到更優(yōu)參數(shù)設(shè)置：SyN步長(zhǎng)為0.25，高斯參數(shù)為(3，0)，在4種尺度下，每一種最多進(jìn)行219次迭代．同時(shí)，將本文算法與基于Elastix[16]工具包中提供的B樣條(B-Spline)配準(zhǔn)方法進(jìn)行了比較．使用歸一化交叉相關(guān)作為相似性度量，進(jìn)行了5個(gè)分辨率的1000次迭代．另外，還與基于學(xué)習(xí)的VoxelMorph算法進(jìn)行了比較．VoxelMorph提出了兩種以U-Net為網(wǎng)絡(luò)框架的Vm1和Vm2網(wǎng)絡(luò)結(jié)構(gòu)，本文分別在Vm1和Vm2上進(jìn)行了對(duì)比實(shí)驗(yàn)．最初的VoxelMorph算法被應(yīng)用于腦組織配準(zhǔn)，為保證公平對(duì)比，在EMPIRE10數(shù)據(jù)集上重新對(duì)VoxelMorph進(jìn)行了訓(xùn)練．最后，為驗(yàn)證Inception模塊的有效性，筆者將本文方法與未添加Inception模塊的U-Net進(jìn)行了對(duì)比實(shí)驗(yàn)．

2.5?實(shí)現(xiàn)細(xì)節(jié)

3?實(shí)驗(yàn)結(jié)果

通過(guò)使用測(cè)試數(shù)據(jù)集中包含的3000對(duì)地標(biāo)點(diǎn)計(jì)算TRE．TRE結(jié)果如表1所示，括號(hào)中為標(biāo)準(zhǔn)偏差.本文提出的方法取得的平均配準(zhǔn)誤差為2.09mm，標(biāo)準(zhǔn)偏差為1.55mm，取得了次優(yōu)結(jié)果．

表1?不同算法在DIR-Lab4DCT數(shù)據(jù)集上的TRE

Tab.1?TRE of different algorithms on the DIR-Lab 4DCT dataset

注：括號(hào)內(nèi)為標(biāo)準(zhǔn)偏差，無(wú)量綱．

如圖5所示，通過(guò)計(jì)算測(cè)試數(shù)據(jù)集上固定圖像掩膜與扭曲圖像掩膜間的平均Dice得分，相比于VoxelMorph、SyN和B-Spline，本文所提出的方法可實(shí)現(xiàn)最優(yōu)的配準(zhǔn)性能(平均Dice＝0.987)．

能夠進(jìn)行主觀性評(píng)價(jià)的實(shí)驗(yàn)結(jié)果如圖6所示.稠密位移向量場(chǎng)以RGB圖像表示，3個(gè)通道對(duì)應(yīng)位移場(chǎng)的3個(gè)維度，越亮位置表示該位置發(fā)生的變形越大．圖中也給出被稠密位移向量場(chǎng)扭曲的形變網(wǎng)格圖像．值得注意的是，本文提出的方法與VoxelMorph使用的采樣方法相同，與ANTs和Elastix的采樣方式不同．

在稠密位移向量場(chǎng)的雅可比行列式圖像中，比0大的位置映射為藍(lán)色，而小于等于0的位置映射為紅色，即紅色標(biāo)記了出現(xiàn)奇點(diǎn)的位置．在雅可比圖像行列式圖像中，本文方法并未顯現(xiàn)出折疊體素．在RGB圖像中和網(wǎng)格圖像中，能夠看出本文方法獲得了更加豐富的變形．

表2總結(jié)了所有的配準(zhǔn)實(shí)驗(yàn)結(jié)果，并給出了所有方法在GPU和CPU上的運(yùn)行時(shí)間對(duì)比．到目前為止，還沒(méi)有針對(duì)ANTs和Elastix的GPU實(shí)現(xiàn)．結(jié)果表明，本文所提方法的配準(zhǔn)速度相比傳統(tǒng)方法(SyN和BSpline)有著指數(shù)級(jí)提升，與VoxelMorph相近．值得注意的是，本文方法在獲得高配準(zhǔn)精度的同時(shí)幾乎不產(chǎn)生折疊體素．另外，表2中的差異和是絕對(duì)差異圖像中所有體素值的總和，配準(zhǔn)效果越好則差異和越?。畬?shí)驗(yàn)結(jié)果表明本文方法在差異圖像的定量評(píng)價(jià)上取得了最優(yōu)．

表2?各算法在DIR-Lab4DCT數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果

Tab.2?Experimental results of different algorithms on the DIR-Lab 4DCT dataset

注：括號(hào)內(nèi)為各變量的標(biāo)準(zhǔn)偏差，無(wú)量綱．

4?討?論

本文提出了一種基于無(wú)監(jiān)督的三維肺部CT圖像配準(zhǔn)方法．利用浮動(dòng)圖像和固定圖像對(duì)之間的NCC來(lái)優(yōu)化卷積神經(jīng)網(wǎng)絡(luò)參數(shù)，而不需要額外GT.在損失函數(shù)中引入雅可比正則項(xiàng)，有效地減少了扭曲圖像的折疊體素?cái)?shù)量．采用基于3D-TPS變換的數(shù)據(jù)增強(qiáng)方法人工生成大規(guī)模訓(xùn)練數(shù)據(jù)，滿足了卷積神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練需求．實(shí)驗(yàn)結(jié)果表明，該方法取得TRE值為2.09mm，平均Dice得分為0.987，且扭曲圖像中僅包含極少量的折疊體素．測(cè)試集上的實(shí)驗(yàn)結(jié)果也表明該方法在三維肺部CT圖像配準(zhǔn)中具有良好的魯棒性．

5?結(jié)?語(yǔ)

本文提出了一種基于無(wú)監(jiān)督學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)方法實(shí)現(xiàn)三維肺部CT圖像的非剛性配準(zhǔn)．基于3D-TPS變換的數(shù)據(jù)增強(qiáng)能夠解決數(shù)據(jù)資源有限導(dǎo)致的過(guò)擬合問(wèn)題．在跳接之間引入Inception模塊，充分融合多尺度深層特征生成高精度的位移向量場(chǎng)，以獲得更佳配準(zhǔn)效果．實(shí)驗(yàn)結(jié)果表明，相較于傳統(tǒng)方法，本文算法可實(shí)現(xiàn)可靠的配準(zhǔn)性能，且配準(zhǔn)效率得到了指數(shù)級(jí)提升．

［1］Kaus M R，Brock K K，Pekar V，et al. Assessment of a model-based deformable image registration approach for radiation therapy planning[J]. International Journal of Radiation Oncology Biology Physics，2007，68(2)：572-580.

［2］Yu Z H，Lin S H，Balter P，et al. A comparison of tumor motion characteristics between early stage and locally advanced stage lung cancers[J]. Radiotherapy & Oncology，2012，104(1)：33-38.

［3］Avants B B，Epstein C L，Grossman M，et al. Symmetric diffeomorphic image registration with cross-correlation：Evaluating automated labeling of elderly and neurodegenerative brain[J]. Medical Image Analysis，2008，12(1)：26-41.

［4］Thirion J P. Image matching as diffusion process：An analogy with Maxwell’s demons[J]. Medical Image Analysis，1998，2(3)：243-260.

［5］Fan J，Cao X，Yap P T，et al. BIRNet：Brain image registration using dual-supervised fully convolutional networks[J]. Medical Image Analysis，2018，54(1)：193-206.

［6］Hu Yipeng ，Marc M，Eli G，et al. Label-driven weakly-supervised learning for multimodal deformable image registration[EB/OL]. https://arxiv.org/abs/1711. 01666v2，2017-11-05.

［7］Miao S，Wang Z J，Liao R. A CNN regression approach for real-time 2D/3D registration[J]. IEEE Transactions on Medical Imaging，2016，35(5)：1352-1363.

［8］Balakrishnan G，Zhao A，Sabuncu M R，et al. VoxelMorph：A learning framework for deformable medical image registration[J]. IEEE Transactions on Medical Imaging，2019，38(8)：1788-1800.

［9］Zhao S，Dong Y，Chang E，et al. Recursive cascaded networks for unsupervised medical image registration[C]// IEEE/CVF International Conference on Computer Vision. Seoul，Korea，2019：10599-10609.

［10］ Max J，Karen S，Andrew Z，et al. Spatial transformer networks[EB/OL]. https：//arxiv.org/abs/1506.02025，2016-02-04.

［11］ Balakrishnan G，Zhao A，Sabuncu M R，et al. An unsupervised learning model for deformable medical image registration[C]// IEEE/CVF Conference on Computer Vision and Pattern Recognition. Salt Lake City，USA，2018：9252-9260.

［12］ Murphy K，van Ginneken B，Reinhardt J M，et al. Evaluation of registration methods on thoracic CT：The EMPIRE10 challenge[J]. IEEE Transactions on Medical Imaging，2011，30(11)：1901-1920.

［13］ Castillo E，Castillo R，Martinez J，et al. Four-dimensional deformable image registration using trajectory modeling[J]. Physics in Medicine and Biology，2010，55(1)：305-327.

［14］ Vandemeulebroucke J，Rit S，Kybic J，et al. Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs：Spatiotemporal motion estimation for 4D CT[J]. Medical Physics，2010，38(1)：166-178.

［15］ Avants B B，Tustison N J，Song G，et al. A reproducible evaluation of ANTs similarity metric performance in brain image registration[J]. NeuroImage，2011，54(3)：2033-2044.

［16］ Klein S，Staring M，Murphy K，et al. Elastix：A toolbox for intensity-based medical image registration[J]. IEEE Transactions on Medical Imaging，2010，29(1)：196-205.

Research on a 3D Lung Computed Tomography Image Registration Method Based on Unsupervised Learning

Jiang Shan，Zhang Hongyun，Yang Zhiyong，Zhang Guobin

(School of Mechanical Engineering，Tianjin University，Tianjin 300350，China)

Deformable registration of 3D lung CT images is crucial in medical image registration. However，nonlinear deformation and large-scale displacement of lung tissues caused by respiratory motion pose great challenges in the deformable registration of 3D lung CT images. Thus，we present a fast end-to-end registration method based on unsupervised learning. We optimized the classic U-Net model and added Inception modules between skip connections. The Inception module aims to capture and merge information at different spatial scales for generating a high-precision dense displacement vector field. To ensure a smooth displacement vector field，we introduced the Jacobian regularization term into the loss function to directly penalize the singularity of the displacement field during training. The existing publicly available datasets cannot implement model training. To address over-fitting caused by limited data resources and to expand the training data，we proposed a data augmentation method based on a 3D thin plate spline transform. Moreover，6060 CT scans will be generated based on the EMPIRE10 dataset，which contains 60 original CT scans to meet the requirement of convolution neural network training. Regarding the DIR-Lab 4DCT dataset，we achieved a target registration error of 2.09mm，an optimal Dice score of 0.987，and almost no folding voxels in comparison with the experimental results obtained using the deep learning method Voxelmorph and registration packages，such as advanced normalization tools (ANTs) and Elastix.

deformable registration；convolution neural network；unsupervised learning；data augmentation；Jacobian regularization

TP242

A

0493-2137(2022)03-0247-08

10.11784/tdxbz202010040

2020-10-21；

2020-12-10.

姜?杉（1973—??），女，博士，教授.

姜?杉，shanjmri@tju.edu.cn.

國(guó)家自然科學(xué)基金資助項(xiàng)目(51775368，81871457，51811530310)；天津市科技資助項(xiàng)目(18YFZCSY01300)；天津市津南區(qū)科技計(jì)劃資助項(xiàng)目(20200110)．

Supported by the National Natural Science Foundation of China(No.51775368，No.81871457，No.51811530310)，Tianjin Science and Tech-nology Project(No.18YFZCSY01300)，Tianjin Jinnan District Science and Technology Planning Project(No. 20200110).

(責(zé)任編輯：王曉燕)