集合卡爾曼濾波同化估算地表水熱通量的研究

劉翔舸，王鵬新，唐伯惠，黃健熙，陶 欣

(1. 中國(guó)農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院 北京 海淀區(qū) 100083；2. 中國(guó)科學(xué)院地理科學(xué)與資源研究所 北京 朝陽區(qū) 100101；3. 美國(guó)馬里蘭大學(xué)帕克分校地理科學(xué)系 美國(guó) 馬里蘭帕克大學(xué)城 20742)

地表水熱通量在地面表層范圍內(nèi)的大氣、水文和生態(tài)等諸多領(lǐng)域研究中具有重要的地位，其影響著地表水汽和能量的輸送過程。在農(nóng)業(yè)大環(huán)境中，地表水熱通量的估計(jì)對(duì)于農(nóng)業(yè)管理、作物估產(chǎn)以及農(nóng)業(yè)生態(tài)系統(tǒng)的研究具有十分重要的意義。目前，在能量平衡原理、空氣動(dòng)力學(xué)原理、土壤-植被-大氣傳輸理論的不斷發(fā)展下，已經(jīng)建立了多種估算地表水熱通量的模型。如研究者利用建立相關(guān)的能量和物質(zhì)交換過程方法的單層或多層模型，進(jìn)行水熱通量的估算研究[1-2]。

另一方面，陸面數(shù)據(jù)同化方法的發(fā)展使數(shù)據(jù)同化技術(shù)在陸面過程研究中得到較多的應(yīng)用[3]。其中，數(shù)據(jù)同化技術(shù)在實(shí)際應(yīng)用中可同時(shí)利用多源觀測(cè)數(shù)據(jù)結(jié)合動(dòng)態(tài)模型的多層次架構(gòu)的特點(diǎn)，能夠?qū)⒛Ｐ团c觀測(cè)通過數(shù)據(jù)同化算法有機(jī)結(jié)合在一起，共同作用于對(duì)目標(biāo)參數(shù)的估算反演中，從而獲得高分辨率的估算結(jié)果，使最終估算的目標(biāo)數(shù)據(jù)質(zhì)量得到提升。因此，數(shù)據(jù)同化技術(shù)在估算地表水熱通量的相關(guān)研究也不同程度地開展起來，文獻(xiàn)[4]應(yīng)用變分同化方法，將遙感反演的地表溫度數(shù)據(jù)同化入陸面能量平衡與地表溫度為約束條件的動(dòng)態(tài)方程組中，得到了潛熱和地表蒸散參數(shù)的有效估算。文獻(xiàn)[5]應(yīng)用擴(kuò)展卡爾曼濾波方法，并通過設(shè)計(jì)方案中考慮觀測(cè)和模型誤差的不確定性，將衛(wèi)星反演溫度誤差與實(shí)測(cè)表層土壤溫度進(jìn)行擬合，以擬合作為觀測(cè)誤差同化入地表陸面模型中，估計(jì)下層土壤溫度與地表蒸散。文獻(xiàn)[6]以能量平衡方程的方法提出了一種近地表陸面模式，并應(yīng)用變分同化方法與遙感反演的地表溫度耦合，實(shí)現(xiàn)對(duì)地表熱通量的估算。文獻(xiàn)[7]應(yīng)用弱約束的四維變分方法，對(duì)簡(jiǎn)單陸面模式進(jìn)行了數(shù)值試驗(yàn)，將地面站實(shí)測(cè)地表溫度與陸面模式等關(guān)鍵參數(shù)共同建立代價(jià)函數(shù)，實(shí)現(xiàn)地表水熱通量的估計(jì)。文獻(xiàn)[8]采用通用陸面模式結(jié)合EnKF對(duì)地表水熱通量進(jìn)行估算。文獻(xiàn)[9]綜合考慮了能量傳遞過程中土壤表層與植被冠層溫度之間的差異，并應(yīng)用變分同化將地表溫度同化入一個(gè)雙層模型，估算了地表水熱通量數(shù)據(jù)。

由于多數(shù)地表水熱通量的研究工作是基于國(guó)外通量觀測(cè)站數(shù)據(jù)的支持，對(duì)國(guó)內(nèi)研究的參考價(jià)值有限。因此本文嘗試?yán)脟?guó)內(nèi)試驗(yàn)站點(diǎn)的觀測(cè)數(shù)據(jù)，通過搭建一個(gè)以同化地表溫度改進(jìn)陸面過程模型中，地表水熱通量估計(jì)的數(shù)據(jù)同化系統(tǒng)，為數(shù)據(jù)同化方法在農(nóng)田環(huán)境中的地表水熱通量估計(jì)的進(jìn)一步研究提供借鑒與基礎(chǔ)。

1 研究站點(diǎn)概況與數(shù)據(jù)預(yù)處理

本文選用的試驗(yàn)數(shù)據(jù)來自山東禹城農(nóng)業(yè)綜合試驗(yàn)站(36.83N ，116.57E )，如圖1所示，站點(diǎn)位于山東省禹城市，是中國(guó)科學(xué)院首批的開放式試驗(yàn)站，并于1999年由國(guó)家科技部設(shè)立為首批國(guó)家重點(diǎn)試驗(yàn)站。該站地理環(huán)境處于黃淮海平原，地處暖溫帶半濕潤(rùn)季風(fēng)氣候區(qū)。年平均氣溫13.1℃，年均降水量582 mm，光熱資源豐富，雨熱同期，是主要的農(nóng)業(yè)生產(chǎn)區(qū)，站點(diǎn)農(nóng)作物主要是小麥、玉米等，該站點(diǎn)同時(shí)具有農(nóng)業(yè)生產(chǎn)的有利環(huán)境和農(nóng)業(yè)作物類型的代表性。

禹城農(nóng)業(yè)綜合試驗(yàn)站通過不斷發(fā)展，陸續(xù)建立了多種類型儀器的監(jiān)測(cè)與研究系統(tǒng)，可具備同時(shí)監(jiān)測(cè)氣象數(shù)據(jù)、土壤數(shù)據(jù)、植物信息等能力。為農(nóng)業(yè)的持續(xù)發(fā)展和陸面過程中的土壤-植物-大氣間的能量循環(huán)監(jiān)測(cè)研究，提供完備的科技手段和觀測(cè)數(shù)據(jù)支持。

2 EnKF同化系統(tǒng)的構(gòu)建

圖2所示為本文的研究數(shù)據(jù)同化系統(tǒng)的具體流程圖。從圖中可以大致看出，數(shù)據(jù)同化系統(tǒng)的主要結(jié)構(gòu)可分為：外部觀測(cè)、同化算法、動(dòng)態(tài)模型3個(gè)主要部分組成。數(shù)值試驗(yàn)研究中的外部觀測(cè)數(shù)據(jù)與氣象驅(qū)動(dòng)數(shù)據(jù)由上述國(guó)內(nèi)禹城試驗(yàn)站點(diǎn)的實(shí)測(cè)數(shù)據(jù)提供，下面將對(duì)數(shù)據(jù)同化系統(tǒng)的同化算法和陸面過程模型分別進(jìn)行描述。

圖2 同化系統(tǒng)流程圖

2.1 EnKF同化方法

為了有效地應(yīng)對(duì)非線性問題，從而合理得到誤差的估計(jì)，一種集合的卡爾曼濾波(ensemble Kalman filter，EnKF)發(fā)展起來[10-18]。EnKF運(yùn)用Monte-Carlo方法與集合的思想，主要解決了一般同化方法在實(shí)際應(yīng)用中背景誤差協(xié)方差矩陣估計(jì)和預(yù)報(bào)困難的問題；解決非線性系統(tǒng)的近似問題；避免使用變分模式下的伴隨模式；較易考慮模式誤差問題；并可以實(shí)現(xiàn)并行化計(jì)算[19]。因此，本文選擇EnKF方法作為數(shù)據(jù)同化系統(tǒng)中的優(yōu)化算法。

EnKF算法作為應(yīng)用最廣泛的非線性濾波同化方法，最初由文獻(xiàn)[12]提出并應(yīng)用于海洋數(shù)據(jù)同化領(lǐng)域。但是隨著近年來陸面數(shù)據(jù)同化研究的不斷發(fā)展，EnKF方法也逐步應(yīng)用于陸面過程研究中，其EnKF算法由KF算法推衍而來，主要結(jié)構(gòu)仍分為預(yù)報(bào)階段和分析階段兩個(gè)組成部分不變，具體算法描述如下：將動(dòng)態(tài)模型M在每個(gè)隨機(jī)變量第i時(shí)刻的模型的預(yù)報(bào)值xi作為背景場(chǎng)的狀態(tài)集合向量bX ，

2.2 動(dòng)態(tài)模型描述

式中，Rn是地表凈輻射；H和LE分別是顯熱通量和潛熱通量；P是地表的熱慣量； 為每日頻率(1/86 400 s)；Tsi表示前一天第i時(shí)刻的地表溫度值；表示土壤溫度；DN表示一天內(nèi)的時(shí)間步長(zhǎng)。式(11)為該模式的主體部分，其中熱通量數(shù)據(jù)(潛熱、顯熱與凈輻射)依賴于地表溫度Ts的狀態(tài)估計(jì)。因此，準(zhǔn)確估算與優(yōu)化當(dāng)前時(shí)刻的地表溫度Ts，對(duì)下一刻的地表水熱通量數(shù)據(jù)的估計(jì)至關(guān)重要。所以，可認(rèn)為通過改變地表溫度的預(yù)報(bào)，進(jìn)而可改變地表水熱通量的估計(jì)精度。模型中的參數(shù)具體求解與設(shè)定參考文獻(xiàn)[6-7, 20]。

2.3 同化系統(tǒng)中誤差的確定

在數(shù)據(jù)同化系統(tǒng)的研究中，由于系統(tǒng)中模型表述狀態(tài)變量的復(fù)雜性，在實(shí)際應(yīng)用中同化系統(tǒng)的誤差估計(jì)與處理是受到廣泛關(guān)注的問題。數(shù)據(jù)同化的目的是在動(dòng)態(tài)模型的動(dòng)力框架內(nèi)，通過綜合不同來源的不同時(shí)空分辨率的常規(guī)和非常規(guī)觀測(cè)，將動(dòng)態(tài)模型和各種觀測(cè)信息集合成為不斷依靠觀測(cè)而調(diào)整優(yōu)化模型的軌跡，減小模型預(yù)報(bào)誤差而獲得更高精度的預(yù)報(bào)值。因此，誤差的研究應(yīng)始終貫穿于數(shù)據(jù)同化的研究之中[21]。

通常情況認(rèn)為動(dòng)態(tài)模型在運(yùn)行一段時(shí)間后將自動(dòng)達(dá)到一種穩(wěn)定和平衡的狀態(tài)，此時(shí)的模型模擬結(jié)果在一定程度上具有模型自身誤差的代表性。本文根據(jù)該方法，利用運(yùn)行一段時(shí)間后的模型模擬結(jié)果與地面試驗(yàn)站點(diǎn)的實(shí)測(cè)值進(jìn)行比較，可獲得統(tǒng)計(jì)上的模型誤差。

如圖3所示，利用禹城試驗(yàn)站點(diǎn)提供的氣象數(shù)據(jù)使模型從第91天開始，持續(xù)運(yùn)行90天的模擬結(jié)果。圖中所示從第100天左右開始，模型模擬的LST結(jié)果高于地表觀測(cè)的LST結(jié)果，之后在一段時(shí)間內(nèi)(第100天～160天)的觀測(cè)值與模擬值的差別較大，此外，圖中顯示在第165天左右，試驗(yàn)站點(diǎn)實(shí)測(cè)的地表溫度數(shù)值有部分超過320 K，認(rèn)為該部分?jǐn)?shù)據(jù)存在較大的測(cè)量誤差，在以后計(jì)算中以320 K為閾值，對(duì)觀測(cè)站點(diǎn)實(shí)測(cè)的地表溫度數(shù)據(jù)進(jìn)行質(zhì)量控制。

為了研究并得到模型的誤差，認(rèn)為模型自身運(yùn)行10天后的模擬結(jié)果即可達(dá)到穩(wěn)定。將模型在一段時(shí)間內(nèi)的(第100天～160天)對(duì)應(yīng)試驗(yàn)站點(diǎn)的觀測(cè)值與模型獨(dú)立運(yùn)行的模擬值進(jìn)行以天為單位的差值平均，得到該時(shí)段一天內(nèi)的模型誤差變化趨勢(shì)，如圖4所示。圖中的結(jié)果表示，地表溫度的誤差主要體現(xiàn)在白天階段，并且白天的誤差相比夜間變化的波動(dòng)較大。綜合考慮，模型地表溫度平均高于地面實(shí)測(cè)地表溫度約為1.5 K，本文在試驗(yàn)中假設(shè)模型運(yùn)行時(shí)的模型誤差為定值，不隨其他變量相關(guān)，通過上述統(tǒng)計(jì)方式得出平均每天模型誤差的均值，作為同化系統(tǒng)中模型誤差的初始設(shè)定。

圖3 山東禹城站點(diǎn)模擬地表溫度結(jié)果與觀測(cè)結(jié)果比較圖

圖4 一天內(nèi)模擬值與觀測(cè)值的差值變化

3 結(jié)果分析與討論

3.1 集合大小的確定

EnKF集合的多少?zèng)Q定最終整個(gè)數(shù)據(jù)同化系統(tǒng)的效率和精度。集合數(shù)目過多影響系統(tǒng)執(zhí)行效率，過少卻影響最終的同化精度。因此在本文試驗(yàn)中，考慮利用調(diào)整集合m的大小，通過系統(tǒng)運(yùn)行結(jié)果的均方根誤差(root mean square error，RMSE)變化，試確定適合本同化系統(tǒng)的集合數(shù)目m。

在確定集合大小的試驗(yàn)中，同化步長(zhǎng)，模型誤差、背景場(chǎng)誤差、模型驅(qū)動(dòng)參數(shù)等均取值相同，只調(diào)整集合m的大小，將得到的同化后地表溫度結(jié)果與試驗(yàn)站觀測(cè)的地表溫度對(duì)比，利用誤差RMSE變化分析EnKF的集合數(shù)目的影響，如圖5所示。

由圖5可以看出，RMSE變化趨勢(shì)隨著集合數(shù)目的增大顯著減小，當(dāng)集合大小為30時(shí)RMSE最小，之后逐漸趨于穩(wěn)定。因此取m=30。

3.2 同化結(jié)果

利用已設(shè)定的數(shù)據(jù)同化系統(tǒng)，對(duì)山東禹城2009年4～6月內(nèi)進(jìn)行EnKF的數(shù)據(jù)同化試驗(yàn)。考慮數(shù)據(jù)的連續(xù)性和演示效果的清晰性，選擇試驗(yàn)站點(diǎn)的第91天開始之后30天的數(shù)據(jù)同化結(jié)果進(jìn)行展示。

圖5 集合大小對(duì)同化系統(tǒng)的影響

圖6所示為模型模擬的LST和同化后的LST的趨勢(shì)變化，兩者總體趨勢(shì)與實(shí)測(cè)的LST變化相同，但是從圖中可以明顯看出，同化后的LST趨勢(shì)更接近實(shí)測(cè)的LST曲線，說明相較于模型模擬的LST，同化后的LST精度得到了有效的改善。其中，模型模擬LST與地面實(shí)測(cè)LST的RMSE為4.85 K，經(jīng)過同化后LST的RMSE則下降到1.57 K。證明在同化過程中對(duì)于LST估計(jì)精度的改善較為明顯。因此，伴隨在整個(gè)同化的過程中，改善的LST將使模型最終估計(jì)的地表水熱通量的精度同時(shí)得到提高。

圖6 山東禹城站點(diǎn)地表溫度、地表水熱通量的同化結(jié)果

另一方面，由于試驗(yàn)中缺少對(duì)地面觀測(cè)通量數(shù)據(jù)的支持，為了驗(yàn)證數(shù)據(jù)同化后的通量結(jié)果，利用MODIS的數(shù)據(jù)產(chǎn)品進(jìn)行輔助驗(yàn)證。將同化后得到的水熱通量中的潛熱通量與MODIS的ET產(chǎn)品(MOD16A2)進(jìn)行對(duì)比分析。MODIS的ET產(chǎn)品分為8天和月的ET/LE的全球產(chǎn)品，分辨率為1 km′1 km。驗(yàn)證試驗(yàn)選擇MODIS對(duì)應(yīng)禹城站點(diǎn)的8天ET產(chǎn)品與陸面模型獨(dú)自模擬的ET(未同化)和同化系統(tǒng)估算的ET進(jìn)行對(duì)比驗(yàn)證，驗(yàn)證試驗(yàn)將分別從ET變化趨勢(shì)和RMSE的結(jié)果兩方面做進(jìn)一步的判斷與分析。

MOD16的ET/LE產(chǎn)品采用8天累加的ET值。文中的數(shù)據(jù)均為30 m ins，因此需要對(duì)數(shù)據(jù)進(jìn)行時(shí)間尺度上的匹配。根據(jù)MODIS產(chǎn)品的定義，本文將模型模擬的ET與同化結(jié)果的ET，在8天內(nèi)的數(shù)據(jù)累加生成8天的數(shù)據(jù)產(chǎn)品(mm/8 d)，與MOD16產(chǎn)品對(duì)比結(jié)果。

圖7所示為2009年第89天～177天(4～6月)間的8天的ET結(jié)果隨時(shí)間的變化的趨勢(shì)情況，圖中的3個(gè)方法得到的ET的總體變化趨勢(shì)相同。由圖中可以明顯看出，模型模擬得到的ET值總體偏高(與MOD16的ET之間的RMSE=4.18 mm)。

圖7 同化后ET結(jié)果對(duì)比圖

而經(jīng)過同化后的ET整體趨勢(shì)與模型相比略低，但ET值仍高于MODIS的ET產(chǎn)品(與MOD16的ET之間的RMSE=2.99 mm)。綜合LST與ET結(jié)果的RMSE對(duì)比試驗(yàn)，認(rèn)為經(jīng)過數(shù)據(jù)同化后的模型模擬結(jié)果(LST、H和LE)均得到了有效改善。因此，本文構(gòu)建的利用地表溫度改進(jìn)地表水熱通量估計(jì)數(shù)據(jù)同化系統(tǒng)具有可應(yīng)用性與適用性。

4 結(jié) 論

本文利用山東禹城試驗(yàn)站點(diǎn)數(shù)據(jù)通過EnKF同化方法，改進(jìn)模型地表水熱通量估算精度的嘗試是成功的，具有較高的可重復(fù)性與穩(wěn)定性，并具備一定的適用性。在EnKF的數(shù)據(jù)同化系統(tǒng)中，探討了系統(tǒng)中陸面模型誤差的設(shè)定問題，認(rèn)為一定時(shí)間內(nèi)模型誤差的變化可用統(tǒng)計(jì)的方法標(biāo)定，進(jìn)而運(yùn)用于整個(gè)數(shù)據(jù)同化過程。同時(shí)，探討了EnKF算法中集合大小對(duì)同化結(jié)果的影響，通過對(duì)集合大小變化的誤差協(xié)方差RMSE的比較分析，確定已建立同化系統(tǒng)的集合數(shù)目。最后，利用該同化系統(tǒng)同化地表溫度數(shù)據(jù)，提高了模型的地表水熱通量的最終估計(jì)精度，認(rèn)為該同化系統(tǒng)方案可有效地改善陸面過程模型對(duì)地表水熱通量的估計(jì)。該研究方法和結(jié)果在農(nóng)田地表水熱通量的估算方面具有參考價(jià)值。

感謝中國(guó)科學(xué)院地理科學(xué)與資源研究所提供山東禹城試驗(yàn)站觀測(cè)數(shù)據(jù)。

[1] 劉紹民, 李小文, 施生錦, 等. 大尺度地表水熱通量的觀測(cè)、分析與應(yīng)用[J]. 地球科學(xué)進(jìn)展, 2010, 25(11): 1113-1127.LIU Shao-m in, LI Xiao-wen, SHI Sheng-jin, et al.Measurement, analysis and application of surface energy and water vapor fluxes at large scale[J]. Advances in Earth Science, 2010, 25(11): 1113-1127.

[2] 辛?xí)灾? 田國(guó)良, 柳欽火. 地表蒸散定量遙感的研究進(jìn)展[J]. 遙感學(xué)報(bào), 2003(3): 233-240.

XIN Xiao-zhou, TIAN Guo-liang, LIU Qin-huo. A review of researches on remote sensing of land surface evapotranspiration[J]. Journal of Remote Sensing, 2003(3):233-240.

[3] LIANG S. Advances in land remote sensing, system,modeling, inversion and application[M]. New York, USA:Springer, 2008.

[4] BONI G, ENTEKHABI D, CASTELLI F. Land data assimilation w ith satellite measurements for the estimation of surface energy balance components and surface control on evaporation[J]. Water Resources Research, 2001, 37(6):1713-1722.

[5] KUMAR P, KALEITA A L. Assimilation of near-surface temperature using extended Kalman filter[J]. Advances in Water Resources, 2003, 26(1): 79-93.

[6] CAPARRINI F, CASTELLI F, ENTEKHABI D. Estimation of surface turbulent fluxes through assimilation of radiometric surface temperature sequences[J]. Journal of Hydrometeorology, 2004, 5(1): 145-159.

[7] JUN Q, LIANG S L, LIU R G, et al. A weakconstraint-based data assimilation scheme for estimating surface turbulent fluxes[J]. Geoscience and Remote Sensing Letters, 2007, 4(4): 649-653.

[8] XU T R, LIANG S L, LIU S M. Estimating turbulent fluxes through assimilation of geostationary operational environmental satellites data using ensemble Kalman filter[J]. Journal of Geophysical Research-Atmospheres,2011: 116(D09): 109-125.

[9] BATENI S M, LIANG S. Estimating surface energy fluxes using a dual-source data assimilation approach adjoined to the heat diffusion equation[J]. Journal of Geophysical Research-Atmospheres, 2012, 117(D17): 118-132.

[10] EVENSEN G. Sequential data assimilation w ith a nonlinear quasi-geostrophic model using monte-carlo methods to forecast error statistics[J]. Journal of Geophysical Research-Oceans, 1994, 99(C5): 10143-10162.

[11] EVENSEN G. Advanced data assimilation for strongly nonlinear dynamics[J]. Monthly Weather Review, 1997,125(6): 1342-1354.

[12] EVENSEN G. The ensemble Kalman filter: theoretical formulation and practical implementation[J]. Ocean Dynamics, 2003, 53(4): 343-367.

[13] BURGERS G, VAN LEEUWEN P J, EVENSEN G.Analysis scheme in the ensemble Kalman filter[J].1998(126): 1719-1724.

[14] BISHOP C H, TOTH Z. Ensemble transformation and adaptive observations[J]. Journal of the Atmospheric Sciences, 1999, 56(11): 1748-1765.

[15] BOCQUET M. Ensemble Kalman filtering w ithout the intrinsic need for inflation[J]. Nonlinear Processes in Geophysics, 2011, 18(5): 735-750.

[16] WANG X G, BISHOP C H, JULIER S J. Which is better,an ensemble of positive-negative pairs or a centered spherical simplex ensemble?[J]. Monthly Weather Review,2004, 132(7): 1590-1605.

[17] WHITAKER J S, HAM ILL T M. Ensemble data assim ilation w ithout perturbed observations[J]. Monthly Weather Review, 2006, 134(6): 1722-1722.

[18] WANG X G, HAM ILL T M, WHITAKER J S, et al. A comparison of the hybrid and EnSRF analysis schemes in the presence of model errors due to unresolved scales[J].Monthly Weather Review, 2009, 137(10): 3219-3232.

[19] 黃健熙, 武思杰, 劉興權(quán), 等. 基于遙感信息與作物模型集合卡爾曼濾波同化的區(qū)域冬小麥產(chǎn)量預(yù)測(cè)[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2012, 28(4): 142-148.

HUANG Jian-xi, WU Si-jie, LIU Xing-quan, et al.Regional w inter wheat yield forecasting based on assim ilation of remote sensing data and crop grow th model w ith ensemble Kalman method[J]. Transactions of the Chinese Society of Agricultural Engineering, 2012, 28(4):142-148.

[20] CAPARRINI F, CASTELLI F, ENTEKHAB D. Variational estimation of soil and vegetation turbulent transfer and heat flux parameters from sequences of multisensor imagery[J].Water Resources Research, 2004, 40(12): 5-20.

[21] 擺玉龍, 李新, 韓旭軍. 陸面數(shù)據(jù)同化系統(tǒng)誤差問題研究綜述[J]. 地球科學(xué)進(jìn)展, 2011, 26(8): 795-804.

BAI Yu-long, LI Xin, HAN Xu-jun. A review of error problems for land data assim ilation systems[J]. Advances in Earth Science, 2011, 26(8): 795-804.

編 輯 黃 莘