沒(méi)有它們，給你引力波你也不認(rèn)識(shí)！ | 賽先生天文
● ● ●
封面圖：目前LIGO-Virgo官方證認(rèn)的50個(gè)引力波信號(hào)，包含了GWTC-1和GWTC-2，圖片來(lái)源：LIGO-Virgo | Frank Elavsky, Aaron Geller | Northwestern
撰文 | 吳仕超、趙天宇（北京師范大學(xué)）
責(zé)編 | 韓越揚(yáng)、呂浩然
通常來(lái)講，飛機(jī)和雷達(dá)之間是“貓和老鼠”般的競(jìng)爭(zhēng)關(guān)系：雷達(dá)發(fā)射電磁波，并利用回波尋找飛機(jī)，而回波會(huì)因?yàn)樾盘?hào)衰減而變得微弱。于是，雷達(dá)探測(cè)飛機(jī)則會(huì)變成一個(gè)典型的弱信號(hào)、強(qiáng)噪聲的問(wèn)題。人們由此發(fā)展出了相應(yīng)的匹配濾波（matched filtering）數(shù)據(jù)處理方法，以提高信號(hào)的信噪比[1]，使之更容易被探測(cè)到。雷達(dá)接收到的回波和其發(fā)射的電磁波具有基本相同的波形，所以雷達(dá)的匹配濾波可以直接實(shí)現(xiàn)。
而引力波的探測(cè)也存在相似的地方：從數(shù)據(jù)分析的角度來(lái)講，引力波探測(cè)也是典型的弱信號(hào)、強(qiáng)噪聲問(wèn)題。相應(yīng)地，人們自然也想到了通過(guò)使用匹配濾波的方法來(lái)增強(qiáng)引力波的探測(cè)能力。
為了實(shí)現(xiàn)匹配濾波，我們需要預(yù)先知道信號(hào)的形狀，也就是波形（waveform）。由于引力波探測(cè)器是被動(dòng)接收引力波信號(hào)，因此它的波形不能預(yù)先得知。對(duì)此，人們的應(yīng)對(duì)辦法是從理論上，針對(duì)各種引力波波源構(gòu)造引力波的波形庫(kù)（template bank），也就是所謂的引力波模板（template）。
因?yàn)閷?shí)際的引力波信號(hào)，波形受波源參數(shù)影響且參數(shù)未知，所以引力波模板庫(kù)通常需要覆蓋整個(gè)可能的波源參數(shù)空間，一個(gè)模板庫(kù)通常包含數(shù)百萬(wàn)個(gè)引力波模板。不準(zhǔn)確的引力波模板會(huì)與實(shí)際引力波信號(hào)波形不吻合，損失匹配濾波信噪比，從而導(dǎo)致引力波信號(hào)丟失。所以引力波波源建模至關(guān)重要。
圖1：GW150914在引力波探測(cè)器中的波形及軌道演化示意圖，圖片來(lái)源：文獻(xiàn)[2]。
在引力波被實(shí)際探測(cè)到之前，雙星并合系統(tǒng)是科學(xué)家預(yù)測(cè)的最有可能的波源，同時(shí)也是引力波天文學(xué)研究中最重要的引力波波源之一。因此，它也是人們首先想到要構(gòu)建引力波模板的系統(tǒng)。
雙星系統(tǒng)中的雙黑洞系統(tǒng)不涉及復(fù)雜的潮汐相互作用，是理論描述上最“干凈”的引力波波源系統(tǒng)。為了獲得雙黑洞系統(tǒng)的引力波波形，人們需要也只需要求解愛(ài)因斯坦場(chǎng)方程即可。但愛(ài)因斯坦場(chǎng)方程的復(fù)雜性是出了名的！一兩千頁(yè)的專著《愛(ài)因斯坦場(chǎng)方程的精確解》不干別的，只為了列舉到現(xiàn)在為止已知的愛(ài)因斯坦場(chǎng)方程的解。
圖2：非常“簡(jiǎn)單”的愛(ài)因斯坦場(chǎng)方程
更“可惜”的是，一兩千頁(yè)中與天文系統(tǒng)有對(duì)應(yīng)關(guān)系的解并不多。換句話說(shuō)，愛(ài)因斯坦場(chǎng)方程的解析求解非常困難。愛(ài)因斯坦本人也深諳其引力場(chǎng)方程的困難程度，他本人以及后來(lái)非常多的科學(xué)家們發(fā)展出弱場(chǎng)、低速、微擾等近似愛(ài)因斯坦場(chǎng)方程的計(jì)算分析辦法。但對(duì)于引力波輻射最強(qiáng)的雙星并合過(guò)程，雙黑洞的相對(duì)速度在臨近并合前會(huì)越來(lái)越快，最后接近光速（如圖1），引力場(chǎng)強(qiáng)超過(guò)地球附近場(chǎng)強(qiáng)6個(gè)量級(jí)，廣義相對(duì)論效應(yīng)強(qiáng)到一切近似計(jì)算方法都失效。
于是，人們想到了數(shù)值方法求解愛(ài)因斯坦場(chǎng)方程，也就是通常人們說(shuō)到的數(shù)值相對(duì)論（numerical relativity）?；谶@個(gè)設(shè)想，基普·索恩（Kip Thorne）早年畫(huà)出了雙黑洞并合引力波波形的示意圖（如圖3）。
圖3：索恩設(shè)想的雙黑洞引力波波形（下方坐標(biāo)軸內(nèi)），圖片來(lái)源：Kip Thorne。
既然都“淪落”到數(shù)值求解愛(ài)因斯坦場(chǎng)方程了，無(wú)非就是寫(xiě)寫(xiě)程序，花費(fèi)些電力資源，喝喝咖啡等著波形出來(lái)就可以了。這也是人們最早對(duì)數(shù)值相對(duì)論的預(yù)期。帶著這個(gè)想法，Susan G. Hahn和Richard W. Lindquist最早在1964年就開(kāi)始嘗試數(shù)值求解雙黑洞的愛(ài)因斯坦場(chǎng)方程。
Hahn是計(jì)算數(shù)學(xué)大家Peter Lax的學(xué)生，當(dāng)時(shí)他在IBM上班。Lindquist是有名的廣義相對(duì)論學(xué)家。等他們寫(xiě)好程序，啟動(dòng)計(jì)算機(jī)，坐下來(lái)準(zhǔn)備喝咖啡的時(shí)候，預(yù)期以外的事情發(fā)生了：程序運(yùn)行十多步就崩潰了, 輸出了一大堆NaN（編者注：NaN在計(jì)算機(jī)科學(xué)中代表一堆意義不明的值），沒(méi)有得到任何有物理意義的結(jié)果。兩人瞪直眼睛，為什么？這件事開(kāi)啟了數(shù)值相對(duì)論學(xué)家們長(zhǎng)達(dá)近半個(gè)世紀(jì)的苦思。這個(gè)問(wèn)題被人們稱之為數(shù)值相對(duì)論的穩(wěn)定性問(wèn)題。
廣義協(xié)變性原理是廣義相對(duì)論的理論基礎(chǔ)之一。該原理從著名的馬赫原理發(fā)展而來(lái)，愛(ài)因斯坦也以之為傲。該原理使得愛(ài)因斯坦場(chǎng)方程以張量形式出現(xiàn)，與數(shù)值計(jì)算所需的偏微分方程僅一紙之隔，人們需要選擇坐標(biāo)系（圖4）把愛(ài)因斯坦場(chǎng)方程寫(xiě)成偏微分方程的形式。
20世紀(jì)70年代，普林斯頓大學(xué)的Kenneth Eppley和德克薩斯大學(xué)奧斯汀分校的Larry Smarr對(duì)于愛(ài)因斯坦場(chǎng)方程數(shù)值計(jì)算的坐標(biāo)選擇做了深入研究，發(fā)現(xiàn)坐標(biāo)選擇會(huì)嚴(yán)重影響數(shù)值計(jì)算的穩(wěn)定性。他們針對(duì)一些特殊物理情形, 充分發(fā)掘其內(nèi)在對(duì)稱性, 發(fā)展了特殊、精致的坐標(biāo)以供數(shù)值計(jì)算使用, 取得了不錯(cuò)的效果。但缺點(diǎn)是這些坐標(biāo)選擇不具有通用性, 無(wú)法對(duì)非對(duì)稱性時(shí)空做一般性推廣。他們還提出了自動(dòng)搜索對(duì)稱性的坐標(biāo)條件等, 但發(fā)現(xiàn)還是不能保證計(jì)算的穩(wěn)定性。
20世紀(jì)80年代，Tsvi Piran對(duì)一些簡(jiǎn)單的模型化問(wèn)題給出了初步的數(shù)值計(jì)算結(jié)果。這種利用時(shí)空對(duì)稱性選擇坐標(biāo)系的方式在引力/規(guī)范對(duì)偶數(shù)值相對(duì)論問(wèn)題中也被經(jīng)常使用。但可惜的是，我們面對(duì)的真實(shí)引力波波源系統(tǒng)往往什么對(duì)稱性都沒(méi)有，比如雙黑洞并合系統(tǒng)。
圖4：Kruskal-Szekeres坐標(biāo)下的史瓦西黑洞，圖片來(lái)源：文獻(xiàn)[3]。
直到20世紀(jì)90年代初，LIGO的初步硬件（initial LIGO）建設(shè)正式啟動(dòng)[4]。眼看引力波探測(cè)的數(shù)據(jù)處理很快就要開(kāi)始，人們對(duì)引力波模板的需求迫在眉睫。人們急紅了眼，猜測(cè)數(shù)值相對(duì)論所遇到的不穩(wěn)定性問(wèn)題會(huì)不會(huì)是由于計(jì)算所用網(wǎng)格不夠精細(xì)導(dǎo)致的？為了把計(jì)算網(wǎng)格鋪得足夠細(xì)，就需要用非常多的計(jì)算節(jié)點(diǎn)并行計(jì)算。為了把復(fù)雜的愛(ài)因斯坦場(chǎng)方程和復(fù)雜的并行計(jì)算綜合到一起，就需要許多人的協(xié)作。
綜上，美國(guó)多個(gè)研究所和高校聯(lián)合組成“雙黑洞大挑戰(zhàn)聯(lián)盟”（The Binary Black Hole Grand Challenge Alliance，圖5）, 計(jì)劃構(gòu)建大型數(shù)值相對(duì)論軟件以實(shí)現(xiàn)對(duì)雙黑洞問(wèn)題的數(shù)值模擬。一款名為Cactus的軟件[5]就是在這個(gè)背景下誕生的，而如今數(shù)值相對(duì)論領(lǐng)域最成熟的開(kāi)源軟件Einstein Toolkit[6]就是基于Cactus發(fā)展起來(lái)的。但該聯(lián)合團(tuán)隊(duì)的研究結(jié)果卻表明，數(shù)值相對(duì)論的不穩(wěn)定性問(wèn)題并非只是計(jì)算分辨率的問(wèn)題。雙黑洞大挑戰(zhàn)聯(lián)盟后來(lái)也就不了了之了。
圖5：上世紀(jì)90年代雙黑洞大挑戰(zhàn)聯(lián)盟的介紹網(wǎng)頁(yè)（http://www.crpc.rice.edu/CRPC/demos/BlackHole/）
數(shù)值計(jì)算不穩(wěn)定是不是因?yàn)樗?jì)算的偏微分方程本身就不穩(wěn)定呢？針對(duì)這個(gè)想法，若干的數(shù)值相對(duì)論工作者從偏微分方程的角度出發(fā)分析計(jì)算所用偏微分方程的雙曲性質(zhì)。人們發(fā)現(xiàn)，雖然各種各樣的用來(lái)數(shù)值計(jì)算的偏微分方程形式都來(lái)自于原來(lái)的愛(ài)因斯坦場(chǎng)方程，但各自的偏微分方程性質(zhì)卻千差萬(wàn)別。
基于對(duì)稱雙曲性質(zhì)的要求，加州理工學(xué)院和康奈爾大學(xué)聯(lián)合數(shù)值相對(duì)論小組提出了廣義調(diào)和坐標(biāo)形式的計(jì)算用愛(ài)因斯坦偏微分方程。這些方程形式在2003年之前就被提出來(lái)了，但數(shù)值相對(duì)論的不穩(wěn)定問(wèn)題依然進(jìn)展緩慢。
Bernd Brügmann教授1996年博士畢業(yè)，論文是圈量子引力。圈量子引力界的著名物理學(xué)家Martin Bojowald和Thomas Thiemann都是他的師弟。酷愛(ài)中國(guó)圍棋的Brügmann教授喜歡通過(guò)簡(jiǎn)單的組合得到復(fù)雜事物，這促成他從圈量子引力到數(shù)值相對(duì)論的研究興趣轉(zhuǎn)變。
不同于雙黑洞大挑戰(zhàn)聯(lián)盟關(guān)于計(jì)算網(wǎng)格問(wèn)題的想法，Brügmann教授把雙黑洞問(wèn)題當(dāng)作多尺度物理問(wèn)題來(lái)看待，他最早把科學(xué)計(jì)算領(lǐng)域的自適應(yīng)網(wǎng)格細(xì)化（adaptive mesh refinement）技術(shù)引入到數(shù)值相對(duì)論（圖6），建立了現(xiàn)在稱之為BAM的數(shù)值相對(duì)論軟件。Brügmann教授喜歡合作，但只限于十人以下的合作。這就是為什么他在德國(guó)馬普引力所工作期間拒絕了把BAM引入到Cactus當(dāng)中。所以Cactus的自適應(yīng)網(wǎng)格細(xì)化大科學(xué)計(jì)算平臺(tái)部分是2003年之后才由Eric Schnetter發(fā)展起來(lái)的，現(xiàn)在被稱為Carpet。
1999年，Brügmann教授基于他發(fā)展的自適應(yīng)網(wǎng)格細(xì)化技術(shù)把雙黑洞的數(shù)值計(jì)算時(shí)間推進(jìn)到了35個(gè)時(shí)間單位（數(shù)值計(jì)算把雙黑洞問(wèn)題無(wú)量綱化后的時(shí)間單位）。以現(xiàn)在計(jì)算雙黑洞約1000個(gè)時(shí)間單位的標(biāo)準(zhǔn)來(lái)看，35個(gè)時(shí)間單位微不足道，但在當(dāng)時(shí)這已經(jīng)是不小的進(jìn)展了。到2000年，Brügmann教授和合作者一起把計(jì)算時(shí)間推進(jìn)到了50個(gè)時(shí)間單位。
圖6：雙黑洞計(jì)算中的自適應(yīng)網(wǎng)格細(xì)化，越靠近黑洞的區(qū)域網(wǎng)格越密集。圖片來(lái)源：文獻(xiàn)[7]
2000年對(duì)引力波發(fā)展是個(gè)特殊的年份，一方面LIGO初步硬件基本建成，眼看就可以開(kāi)展數(shù)據(jù)采集了；另一方面，Kip Thorne這一年60歲大壽。對(duì)引力波探測(cè)的極度樂(lè)觀和對(duì)數(shù)值相對(duì)論發(fā)展的極度悲觀使得索恩在其生日宴會(huì)上說(shuō)“很可能引力波的成功探測(cè)比數(shù)值相對(duì)論的突破還要早日實(shí)現(xiàn)”。當(dāng)然，后來(lái)的發(fā)展和索恩的這句話正好相反。
2000年以后，數(shù)值相對(duì)論繼續(xù)發(fā)展。Brügmann教授與其合作者相繼提出了黑洞穿刺的處理方式和針對(duì)BSSN計(jì)算方程形式的特殊坐標(biāo)選擇條件。這些技術(shù)把雙黑洞計(jì)算時(shí)間推進(jìn)到150個(gè)時(shí)間單位。也許做事精細(xì)是德國(guó)人的特點(diǎn)，Brügmann教授一直“固執(zhí)”地認(rèn)為黑洞穿刺點(diǎn)存在坐標(biāo)奇點(diǎn)，處理坐標(biāo)奇點(diǎn)的特殊函數(shù)使得穿刺點(diǎn)不能動(dòng)。
在人們還在順著B(niǎo)rügmann教授的思路往下思考時(shí)，由加州理工學(xué)院和康奈爾大學(xué)聯(lián)合的數(shù)值相對(duì)論小組正在廣義調(diào)和坐標(biāo)方程形式和譜方法兩個(gè)問(wèn)題上苦苦掙扎。數(shù)值相對(duì)論界有名的“獨(dú)狼”——Frans Pretorius于2003年從Matthew Choptuik處學(xué)成畢業(yè)到加州理工學(xué)院進(jìn)行博士后研究。
獨(dú)狼就是獨(dú)狼，一個(gè)人把調(diào)和坐標(biāo)方程形式套用到他在博士期間發(fā)展的自適應(yīng)網(wǎng)格細(xì)化軟件平臺(tái)上。2005年秋，數(shù)值相對(duì)論學(xué)家們從睡夢(mèng)中被人叫醒，“雙黑洞計(jì)算穩(wěn)定性問(wèn)題解決了！”。人們都覺(jué)得難以置信，Pretorius卻報(bào)告說(shuō)：“我就這么一做，沒(méi)問(wèn)題??！”。
當(dāng)人們還在思考到底發(fā)生了什么的時(shí)候，2006年，美國(guó)德克薩斯大學(xué)布朗斯維爾分校的數(shù)值相對(duì)論小組（現(xiàn)在的羅切斯特理工學(xué)院數(shù)值相對(duì)論小組）放棄Brügmann教授固定黑洞穿刺點(diǎn)的想法獲得了成功。同年美國(guó)國(guó)家航空航天局NASA的數(shù)值相對(duì)論小組不小心忘記了固定穿刺的黑洞，計(jì)算出與德克薩斯大學(xué)布朗斯維爾分校類似的結(jié)果。所以他們兩家同時(shí)宣布，“我們也可以穩(wěn)定計(jì)算雙黑洞了！”。和Pretorius一起，他們?nèi)齻€(gè)小組計(jì)算所得的引力波波形對(duì)比結(jié)果如圖7所示。明眼人皆可看出，異常吻合！
圖7：最早實(shí)現(xiàn)雙黑洞穩(wěn)定計(jì)算的三個(gè)數(shù)值相對(duì)論小組所得引力波波形對(duì)比。圖片來(lái)源：文獻(xiàn)[8]
后來(lái)人們認(rèn)識(shí)到，方程形式、坐標(biāo)條件、網(wǎng)格精細(xì)程度、邊界條件、黑洞奇點(diǎn)處理方式和數(shù)值計(jì)算方法等都對(duì)愛(ài)因斯坦場(chǎng)方程數(shù)值計(jì)算的穩(wěn)定性起到了關(guān)鍵性作用。而且這些問(wèn)題的處理方式需要合適地搭配才能保證數(shù)值相對(duì)論計(jì)算的穩(wěn)定性。
直到現(xiàn)在，也沒(méi)有嚴(yán)格的數(shù)學(xué)理論給出數(shù)值相對(duì)論穩(wěn)定計(jì)算的充分必要條件，更別提數(shù)值相對(duì)論的計(jì)算收斂性理論。所以有人會(huì)說(shuō)數(shù)學(xué)意義的數(shù)值相對(duì)論還沒(méi)有被建立起來(lái)，甚至可以說(shuō)還沒(méi)有開(kāi)始。但無(wú)論如何，世界上已經(jīng)有十多個(gè)數(shù)值相對(duì)論小組能夠穩(wěn)定計(jì)算雙黑洞和其他雙星問(wèn)題，并給出相應(yīng)的引力波波形。
實(shí)際上，從數(shù)值相對(duì)論的“能算”到引力波波形模板的建立，中間還隔著遍歷參數(shù)空間的問(wèn)題，僅帶有軌道進(jìn)動(dòng)的圓軌道雙黑洞模板由15個(gè)參數(shù)來(lái)描述，如果是帶有軌道偏心率或者包含中子星的模板則需要更多參數(shù)來(lái)描述，人們把這個(gè)問(wèn)題稱為數(shù)值相對(duì)論問(wèn)題中的維數(shù)災(zāi)難。
同時(shí)，圖3中索恩理想化的三段相接波形的想法在實(shí)際數(shù)據(jù)處理中也會(huì)帶來(lái)若干問(wèn)題。現(xiàn)在，人們對(duì)這個(gè)問(wèn)題的處理方式是，在深度理解數(shù)值相對(duì)論結(jié)果的基礎(chǔ)上，建立半解析波形模型，類似一種“模糊”處理，將精度降低。目前LIGO-Virgo科學(xué)合作組織主要使用的模板EOBNR和IMRPhenom系列都是這樣的波形模型。
由于數(shù)值相對(duì)論波形的計(jì)算成本過(guò)于高昂（一個(gè)波形需要約3萬(wàn)核心小時(shí)的計(jì)算資源，一次引力波信號(hào)搜索需要約1000萬(wàn)個(gè)波形，合計(jì)約10的11次方核心小時(shí)的計(jì)算資源），很難直接用于實(shí)際數(shù)據(jù)分析，生成速度更快且精度尚可的半解析模型就成為了目前引力波數(shù)據(jù)分析的首選。在建立這些模型的過(guò)程中，數(shù)值相對(duì)論的波形會(huì)被當(dāng)做標(biāo)準(zhǔn)波形來(lái)作比對(duì)。更具體的建模細(xì)節(jié)參見(jiàn)文獻(xiàn)[9]。
從數(shù)值相對(duì)論的計(jì)算、“求穩(wěn)”到“模糊”處理，目前，人們還在繼續(xù)對(duì)雙星并合系統(tǒng)引力波波形做更深入的理解，力圖建立完備的雙星并合系統(tǒng)引力波模型和模板庫(kù)，為未來(lái)的空間引力波探測(cè)和下一代的地面引力波探測(cè)提供堅(jiān)實(shí)的理論支持。
除了雙星系統(tǒng)，我們還可能探測(cè)到其它引力波波源，也許廣義相對(duì)論存在適用邊界，我們還需要其它的引力波波源和引力理論所描述的波形模板。數(shù)值相對(duì)論和波形分析建模的故事還會(huì)伴隨引力波天文學(xué)很長(zhǎng)一段時(shí)間。
作者簡(jiǎn)介：
· 吳仕超，北京師范大學(xué)天文系研究助理，KAGRA科學(xué)合作組織成員。研究方向是引力波波形模板構(gòu)建和引力波數(shù)據(jù)分析。
· 趙天宇，北京師范大學(xué)天文系博士研究生。研究方向是數(shù)值相對(duì)論。
參考資料：
1.Matched filter https://en.wikipedia.org/wiki/Matched_filter
2.Abbott B P, Abbott R, Abbott T D, et al. GW150914: Implications for the stochastic gravitational-wave background from binary black holes[J]. Physical review letters, 2016, 116(13): 131102.
3.Misner C W, Thorne K S, Wheeler J A. Gravitation [M]. Macmillan, 1973.
4.Abramovici A, Althouse W E, Drever R W P, et al. LIGO: The laser interferometer gravitational-wave observatory [J]. Science, 1992, 256(5055): 325-333.
5.Cactus http://cactuscode.org/
6.Einstein Toolkit http://einsteintoolkit.org/
7.Brügmann B. Adaptive mesh and geodesically sliced Schwarzschild spacetime in 3+ 1 dimensions[J]. Physical Review D, 1996, 54(12): 7361.
8.Baker J G, Campanelli M, Pretorius F, et al. Comparisons of binary black hole merger waveforms[J]. Classical and Quantum Gravity, 2007, 24(12): S25.
9.蔡榮根, 曹周鍵, 韓文標(biāo). 并合雙星系統(tǒng)的引力波理論模型[J]. 科學(xué)通報(bào), 2016 (14): 1525-1535.
10.Abbott R, Abbott T D, Abraham S, et al. GW190412: Observation of a binary-black-hole coalescence with asymmetric masses [J]. arXiv preprint arXiv:2004.08342, 2020.
11.Abbott R, Abbott T D, Abraham S, et al. GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object [J]. The Astrophysical Journal Letters, 2020, 896(2): L44.
12.Cao Z, Han W B. Waveform model for an eccentric binary black hole based on the effective-one-body-numerical-relativity formalism [J]. Physical Review D, 2017, 96(4): 044028.
13.曹周鍵, 都志輝. 數(shù)值相對(duì)論與引力波天文學(xué)[J]. 中國(guó)科學(xué): 物理學(xué), 力學(xué), 天文學(xué), 2017, 47(001): 55-72.
14.曹周鍵，愛(ài)因斯坦方程和數(shù)值相對(duì)論，《數(shù)學(xué)與人文》叢書(shū)第二十五輯，68頁(yè)，高等教育出版社 (2018).
15.Baumgarte T W, Shapiro S L. Numerical relativity: solving Einstein's equations on the computer [M]. Cambridge University Press, 2010.
16.Cai R G, Cao Z, Guo Z K, et al. The gravitational-wave physics [J]. National Science Review, 2017, 4(5): 687-706.
制版編輯 | Livan
