莫俊超，吳孝槐，舒耀皋

（1. 上?；ぱ芯吭河邢薰?檢測(cè)中心，上海 200062；2. 上?；瘜W(xué)品公共安全工程技術(shù)研究中心，上海 200062；3. 中國(guó)合格評(píng)定國(guó)家認(rèn)可委員會(huì)，北京 100062）

當(dāng)前，我國(guó)環(huán)境管理政策逐漸從末端治理轉(zhuǎn)向風(fēng)險(xiǎn)管理［1-2］。多介質(zhì)模型是評(píng)估化學(xué)物質(zhì)環(huán)境風(fēng)險(xiǎn)的常用方法，按照模型的結(jié)構(gòu)可分為一級(jí)、二級(jí)、三級(jí)和四級(jí)模型［3-4］。在實(shí)際環(huán)境系統(tǒng)中，化學(xué)物質(zhì)的含量受人類活動(dòng)排放、降解和遷移等過(guò)程的影響，常隨時(shí)間發(fā)生改變，因此，能描述化學(xué)物質(zhì)含量變化的多介質(zhì)動(dòng)力學(xué)模型（即四級(jí)模型）正得到越來(lái)越多的研究和應(yīng)用。

多介質(zhì)動(dòng)力學(xué)模型的求解方法可分為數(shù)值解法［5］和解析解法。當(dāng)多介質(zhì)動(dòng)力學(xué)模型中微分方程數(shù)量超過(guò)3個(gè)，或者模型中化學(xué)物質(zhì)的排放速率、相體積、相流速和溫度不恒定時(shí)，模型的解析解很難求出［6］，這時(shí)求解模型的數(shù)值解便成為了唯一的求解方法。

求解多介質(zhì)動(dòng)力學(xué)模型的數(shù)值解需要使用科學(xué)計(jì)算軟件，目前使用較多的軟件包括Matlab［7-10］、Mathematica［11］、C++［12］、Anylogic［13］、VB［14-15］、Python［16］和R［17］等。其中，Matlab、Mathematica、Anylogic和VB，以及部分C++編譯器屬于商業(yè)軟件，價(jià)格昂貴，更新較慢，且對(duì)版權(quán)要求較多，限制了其應(yīng)用范圍。近年來(lái)，免費(fèi)、開(kāi)源且更新迅速的科學(xué)計(jì)算軟件發(fā)展很快，其豐富的包和庫(kù)等功能讓用戶可使用前人已有的成果，不必從零開(kāi)始編寫計(jì)算程序，從而大大減少了編程的工作量。該類軟件以Python和R為代表，但目前未見(jiàn)系統(tǒng)討論其在多介質(zhì)動(dòng)力學(xué)模型建模方面應(yīng)用的研究報(bào)道。

Julia是近10年內(nèi)開(kāi)發(fā)的一種主要面向科學(xué)計(jì)算的軟件，首個(gè)版本發(fā)布于2013年。同Python和R一樣，Julia也是一款自由軟件，并且結(jié)合了已有科學(xué)計(jì)算軟件的優(yōu)點(diǎn)。目前未見(jiàn)應(yīng)用Julia進(jìn)行多介質(zhì)動(dòng)力學(xué)模型建模的報(bào)道。

考慮到Python中算法的劣勢(shì)［18］，本研究選擇Julia和R，通過(guò)構(gòu)建3個(gè)不同復(fù)雜程度的多介質(zhì)動(dòng)力學(xué)模型，對(duì)Julia和R在程序編寫、運(yùn)行時(shí)間、算法選擇、敏感性分析和便捷性等方面進(jìn)行比較，得出針對(duì)不同的多介質(zhì)動(dòng)力學(xué)模型適用的軟件，以期拓展Julia和R在多介質(zhì)動(dòng)力學(xué)模型建模中的應(yīng)用。

1 材料與方法

1.1 軟件和包

使用R 4.0.2和Julia 1.5.2進(jìn)行模型構(gòu)建。R的下載地址為https：//cran.r-project.org/，使用的包有deSolve（版本1.28）、ODEsensitivity（版本1.1.2）和diffeqr（版本1.0.0）；Julia的下載地址為https：//julialang.org/，使用的包有DifferentialEquations.jl［19］（版本6.15.0）、LSODA.jl（版本0.6.2）、Sundials.jl（版本4.3.0）和DiffEqSensitivity.jl（版本6.31.6）。

1.2 EQC模型簡(jiǎn)介

選擇一個(gè)通用模型EQC（equilibrium criterion，平衡判據(jù)）和EQC模型中的標(biāo)準(zhǔn)環(huán)境進(jìn)行討論。EQC模型于1996年被首次提出［20］，現(xiàn)已廣泛應(yīng)用于化學(xué)物質(zhì)的風(fēng)險(xiǎn)篩查和評(píng)估中［21-25］。EQC模型中包括水、大氣、土壤和沉積物4個(gè)環(huán)境介質(zhì)，EQC動(dòng)力學(xué)模型即對(duì)每個(gè)環(huán)境介質(zhì)中化學(xué)物質(zhì)的物質(zhì)的量編寫微分平衡方程：dM/dt=v輸入-v輸出，式中M為化學(xué)物質(zhì)的物質(zhì)的量，t為時(shí)間，v輸入和v輸出分別為輸入速率和輸出速率。其構(gòu)建過(guò)程詳見(jiàn)參考文獻(xiàn)［26］。

1.3 不同復(fù)雜程度動(dòng)力學(xué)模型的構(gòu)建

為了構(gòu)建不同復(fù)雜程度的動(dòng)力學(xué)模型，選擇不同函數(shù)作為EQC動(dòng)力學(xué)模型中的排放速率函數(shù)，同時(shí)將模型的結(jié)束時(shí)間延長(zhǎng)，具體參數(shù)見(jiàn)表1。一般來(lái)說(shuō)，化學(xué)物質(zhì)不會(huì)直接排放進(jìn)入沉積物中，因此將沉積物中的排放速率設(shè)為0。所有模型中步長(zhǎng)均設(shè)為0.1 h。

表1 不同復(fù)雜程度的模型參數(shù)

1.4 物質(zhì)選擇

為了對(duì)比不同性質(zhì)的物質(zhì)在模型中的運(yùn)行情況，選擇4種性質(zhì)相差較大的物質(zhì)進(jìn)行模擬，其參數(shù)詳見(jiàn)表2［26-27］，表中PCB-209為十氯聯(lián)苯，KOW為正辛醇-水分配系數(shù)。

1.5 函數(shù)、算法和計(jì)算機(jī)配置

R中，使用deSolve包中的ode函數(shù)進(jìn)行模型計(jì)算，計(jì)算方法選擇ode函數(shù)中的lsoda（即ode函數(shù)的默認(rèn)算法）、bdf和radau，使用system.time函數(shù)獲取程序的運(yùn)行時(shí)間。Julia中，使用DifferentialEquations.jl包中的solve函數(shù)進(jìn)行模型計(jì)算，計(jì)算方法選擇lsoda（）（即對(duì)應(yīng)R中ode函數(shù)默認(rèn)算法lsoda的Julia算法）、AutoTsit5（Rosenbrock23（））和CVODE_BDF（），其中l(wèi)soda（）算法需要加載LSODA.jl包，CVODE_BDF（）算法需要加載Sundials.jl包。Julia中使用@time命令獲取程序的運(yùn)行時(shí)間。計(jì)算所用的計(jì)算機(jī)處理器為Intel?Core i7-8565U，內(nèi)存為8 GB，操作系統(tǒng)為Microsoft?Windows 10專業(yè)版。

2 結(jié)果與討論

2.1 程序編寫

由于所編寫的程序較復(fù)雜，不再列出。R的主體程序（即不包括參數(shù)賦值計(jì)算的語(yǔ)句）為14行，Julia的主體程序?yàn)?6行，R較Julia程序更簡(jiǎn)潔。R程序只經(jīng)過(guò)修改一些基本語(yǔ)法便可在Julia上運(yùn)行，反之亦然，兩者語(yǔ)法非常類似，可移植性較強(qiáng)。

2.2 模型運(yùn)行時(shí)間

對(duì)于每個(gè)程序，運(yùn)行5次，計(jì)算得出運(yùn)行時(shí)間的平均值，詳見(jiàn)表3。

表3 R和Julia中EQC動(dòng)力學(xué)模型的運(yùn)行時(shí)間

2.3 模型敏感性分析方法和運(yùn)行時(shí)間

模型的敏感性分析方法主要分為全局敏感性分析和局部敏感性分析兩個(gè)方面，具體可參考陳衛(wèi)平等［28］的總結(jié)。在R和Julia中，全局敏感性分析方法均有已發(fā)布的包可使用，用戶不必自行編寫敏感性分析程序，使用非常方便。R中模型全局敏感性分析包ODEsensitivity，提供了Morris法和Sobol法兩種敏感性分析方法；Julia中DiffEqSensitivity.jl包可用來(lái)進(jìn)行全局敏感性分析，包含了Morris、Sobol、eFAST和Regression 4種方法。另外，Julia中還提供了局部敏感性分析的函數(shù)，而R中暫未見(jiàn)相關(guān)功能。使用R和Julia中Morris和Sobol方法對(duì)表2中4個(gè)半衰期的敏感性進(jìn)行分析，每個(gè)參數(shù)的值浮動(dòng)±10%并設(shè)為均勻分布，計(jì)算模型敏感性分析的運(yùn)行時(shí)間。求解算法均選擇lsoda，Sobol法中模擬次數(shù)選擇1 000次，結(jié)果見(jiàn)表4。

表2 所選4種物質(zhì)的性質(zhì)參數(shù)

表4 R和Julia中EQC動(dòng)力學(xué)模型敏感性分析的運(yùn)行時(shí)間

2.4 軟件的比較

2.4.1 運(yùn)行時(shí)間

由表3可見(jiàn)，對(duì)于相同的模型，Julia的運(yùn)行時(shí)間至多為R的約1/10，尤其在復(fù)雜模型的lsoda算法中運(yùn)行時(shí)間約為R的1/100，計(jì)算效率非常高。R的主要應(yīng)用領(lǐng)域?yàn)閿?shù)據(jù)統(tǒng)計(jì)，在科學(xué)計(jì)算方面的計(jì)算效率并沒(méi)有進(jìn)行很好的優(yōu)化；Julia主要的應(yīng)用領(lǐng)域?yàn)榭茖W(xué)計(jì)算，從本研究的結(jié)果可見(jiàn)，其在多介質(zhì)動(dòng)力學(xué)模型計(jì)算中的計(jì)算效率是R的至少10倍。在敏感性分析方面，由表4可見(jiàn)，敏感性分析比模型求解需要更長(zhǎng)的運(yùn)行時(shí)間，Julia的運(yùn)行效率同樣是R的至少10倍。另外，不同物質(zhì)在同一模型中的運(yùn)行時(shí)間差別很小，這表明，化學(xué)物質(zhì)性質(zhì)不影響模型的運(yùn)行時(shí)間。

2.4.2 算法

目前R中ode函數(shù)提供了包括lsoda、bdf和radau等算法在內(nèi)的17種算法，可計(jì)算剛性和非剛性模型，默認(rèn)算法lsoda可自動(dòng)選擇剛性或非剛性求解器。多介質(zhì)動(dòng)力學(xué)模型中不同介質(zhì)間化學(xué)物質(zhì)含量往往相差很大，大多數(shù)情況下模型為剛性模型，這時(shí)不能使用rk、rk4或euler等求解非剛性模型的算法進(jìn)行計(jì)算，故R中適合多介質(zhì)動(dòng)力學(xué)模型求解的算法不多。本研究中選擇了效率最高的3種算法，對(duì)比來(lái)看，lsoda和bdf算法的運(yùn)行時(shí)間最短，計(jì)算效率最高，是R中求解多介質(zhì)動(dòng)力學(xué)模型的首選算法，但兩者求解復(fù)雜模型的運(yùn)行時(shí)間均很長(zhǎng)。

Julia中求解微分方程模型的算法非常多，可用來(lái)進(jìn)行多介質(zhì)動(dòng)力學(xué)模型求解的算法超過(guò)10種，為用戶提供了豐富的選擇。除了本研究中使用的lsoda（）、AutoTsit5（Rosenbrock23（））和CVODE_BDF（） 3種算法外，還可以嘗試QNDF（）、TRBDF2（）和RadauIIA（）等多種算法。由表3可見(jiàn)，lsoda（）算法的運(yùn)行時(shí)間最少，均不超過(guò)5 s，是Julia中是求解多介質(zhì)動(dòng)力學(xué)模型的首選算法。在敏感性分析方法中，Julia同樣提供了豐富的選擇，而R中目前僅提供了兩種方法。

2.4.3 便捷性

由以上討論可知，從運(yùn)行時(shí)間和算法方面比較，Julia相對(duì)R有很大優(yōu)勢(shì)。然而，Julia發(fā)展時(shí)間較短，很多軟件暫未提供Julia的接口，其應(yīng)用生態(tài)暫沒(méi)有R完善；Julia的中文資料仍然較少，中文社區(qū)也相對(duì)較??；目前國(guó)內(nèi)Julia服務(wù)器很少，連接不穩(wěn)定，Julia的安裝相對(duì)R更繁瑣。綜上，Julia的應(yīng)用暫不如R便捷。

為了結(jié)合R的便捷性和Julia的計(jì)算優(yōu)勢(shì)，RACKAUCKAS［29］在2020年8月首次發(fā)布了R中的diffeqr包，可令用戶在R中直接調(diào)用Julia的微分方程計(jì)算函數(shù)，R中的JuliaCall包也可直接調(diào)用Julia中的其他函數(shù)。但目前diffeqr和JuliaCall包暫不完善，在針對(duì)復(fù)雜模型進(jìn)行快速高效的模型計(jì)算和分析時(shí)，Julia仍然是較好的選擇。

2.4.4 與其他軟件的比較

結(jié)合以上結(jié)果和參考文獻(xiàn)［18］，對(duì)多介質(zhì)動(dòng)力學(xué)模型求解中常用的軟件R、Julia、Matlab、SciPy、Mathematica、C++、Anylogic和VB，從是否免費(fèi)、算法選擇、計(jì)算效率、模型分析工具、便捷性等方面進(jìn)行比較，詳見(jiàn)表5。綜合來(lái)看，Julia在多介質(zhì)動(dòng)力學(xué)模型求解中優(yōu)勢(shì)較大，對(duì)于復(fù)雜的模型，建議首先選用Julia；但由于Julia的便捷性不如R，對(duì)于簡(jiǎn)單的多介質(zhì)動(dòng)力學(xué)模型，R也是可以滿足需求的，這時(shí)使用R更方便。

表5 多介質(zhì)動(dòng)力學(xué)模型求解中的常用軟件比較

3 結(jié)論

a）分別使用R和Julia，編寫了4種化學(xué)物質(zhì)在EQC動(dòng)力學(xué)模型中的求解程序和半衰期的敏感性分析程序，得到模型求解程序和敏感性分析程序的運(yùn)行時(shí)間。結(jié)果表明，Julia的計(jì)算效率是R的至少10倍，在一些情況下Julia的計(jì)算效率甚至是R的近100倍?；瘜W(xué)物質(zhì)性質(zhì)不影響模型的運(yùn)行時(shí)間。lsoda是Julia中求解多介質(zhì)動(dòng)力學(xué)模型的首選算法，而lsoda和bdf是R中的首選算法。

b）對(duì)多介質(zhì)動(dòng)力學(xué)模型求解中常用的軟件，從是否免費(fèi)、算法選擇、計(jì)算效率、模型分析工具、便捷性等方面進(jìn)行了比較。對(duì)于簡(jiǎn)單模型，建議使用R來(lái)建模、求解和分析；對(duì)于復(fù)雜模型，則建議使用Julia。

c）目前我國(guó)對(duì)多介質(zhì)動(dòng)力學(xué)模型的研發(fā)投入較少，缺乏適用于我國(guó)本土環(huán)境的多介質(zhì)動(dòng)力學(xué)模型，結(jié)合本研究的結(jié)果，建議考慮使用開(kāi)源軟件Julia和R來(lái)進(jìn)行我國(guó)多介質(zhì)動(dòng)力學(xué)模型的研發(fā)和應(yīng)用工作。