基于不同失穩判據的黃土邊坡穩定性數值計算分析

王昕宇 鐘佩文 張慧莉 陳卓 田堪良

摘 要：安全系數是邊坡穩定性評價的一項重要指標，依據不同失穩判據計算得到的安全系數存在差異，會直接影響邊坡穩定性狀態的評價。以萬花山治溝造地工程臺階型開挖邊坡為對象，運用ＦＬＡＣ３Ｄ 有限差分軟件，采用強度折減法，選擇計算不收斂判據、塑性區貫通判據、位移突變判據３ 種失穩判據，對黃土邊坡安全系數進行數值計算。對比分析３ 種判據所得的分析結果及優缺點表明：３ 種失穩判據得出的安全系數大小排序為位移突變判據＞計算不收斂判據＞塑性區貫通判據。計算不收斂判據使用簡便，求解快速，但不能展示折減過程坡體塑性區的變化情況；塑性區貫通判據不宜單獨作為邊坡失穩判據，需要結合不平衡力綜合評估，但可以展示坡體剪切和拉伸塑性區；位移突變判據求解過程復雜，但能真實地反映坡體位移變化導致的臨界破壞狀態。建議類似工程優先選用求解過程快速的計算不收斂判據，然后結合塑性區貫通情況與坡體位移變化狀態對邊坡穩定性進行綜合評價。

關鍵詞：黃土邊坡；失穩判據；安全系數；強度折減法；穩定性分析

中圖分類號：ＴＶ２２３．６；ＴＵ４５７ 文獻標志碼：Ａ ｄｏｉ：１０．３９６９／ ｊ．ｉｓｓｎ．１０００－１３７９．２０２４．０３．０２６

引用格式：王昕宇，鐘佩文，張慧莉，等．基于不同失穩判據的黃土邊坡穩定性數值計算分析［Ｊ］．人民黃河，２０２４，４６（３）：１４３－１４７，１５５．

強度折減法提出后在邊坡穩定性分析方面得到了廣泛應用［１－３］ ?；趶姸日蹨p法得到的安全系數往往依賴于失穩判據的選擇［４］ 。邊坡失穩判據主要有計算不收斂判據、塑性區貫通判據、位移突變判據３種［５］ 。目前對于邊坡失穩判據尚未達成共識，導致不同研究人員針對同一個邊坡得出的安全系數有很大差異［６－８］ 。邊坡穩定性研究需要結合實際工況進行分析，針對實際工況選用何種失穩判據得出的安全系數更為可靠，一直是學界研究的熱點。王其寬等［９］ 基于３ 種不同失穩判據，針對露天礦山堆載邊坡的特殊情況研究得出，計算不收斂判據對應的安全系數較大，位移突變判據對應的安全系數較小，需要結合３ 種判據所得結果綜合評價露天礦山邊坡的穩定性。陳林杰等［１０］ 通過對地鐵車站隧道交叉段的土體進行穩定性研究，得出采用位移突變判據與塑性區貫通判據相結合的方法，可以較準確地得到工況的安全系數。王鵬等［１１］ 采用３ 種失穩判據對貴州仁遵某在建公路橋梁岸坡進行穩定性分析，收斂性判據得到的安全系數最大，塑性區貫通判據次之，位移突變性判據所得結果最小，宜以位移突變性判據進行穩定性分析。郝進鋒等［１２］ 研究單洞室穩定性，得出單洞室塑性區圍繞開挖面存在，計算不收斂判據所得安全系數比位移突變判據的大，塑性區貫通判據所得安全系數比位移突變判據的小，宜采用位移突變判據求解安全系數。邵冠慧等［１３］ 總結分析國內外對于隧道穩定性的判別方法后，得出對于不同的地質條件，失穩判據的選擇有所不同。

綜合上述研究發現，在失穩判據選擇方面，針對不同的實際工況，學者進行了大量的研究，然而，目前采用３ 種失穩判據對黃土邊坡進行穩定性分析的研究還比較少。近年來，發軔于延安的治溝造地工程，是針對黃土高原區特殊的丘陵溝壑地貌啟動的一項治溝保生態、造地惠民生的系統工程［１４］ 。但在工程實施過程中，會形成大量裸露的黃土邊坡，這些邊坡在遇到工程切坡以及坡頂荷載超載等情況時會發生失穩，從而導致滑坡的發生［１５］ 。因此，針對黃土邊坡，選擇何種失穩判據進行穩定分析就顯得尤為重要。

ＦＬＡＣ３Ｄ 有限差分軟件以單元格為基本單元展示模擬材料的剪切、拉伸破壞情況，在模擬巖體材料的受力破壞和變形方面應用廣泛［１６］ 。因此，本研究選用ＦＬＡＣ３Ｄ 有限差分軟件，根據陜北黃土高原區萬花山治溝造地工程中臺階型開挖邊坡的實際情況建立幾何模型，采用３ 種失穩判據，逐一對黃土邊坡穩定性進行數值計算，對比分析３ 種判據評價結果的差異及優缺點，為實際工況中黃土邊坡穩定性評價提供技術方法。

１ 參數選取與模型建立

１．１ 參數選取

試驗所用土樣采自陜西省延安市萬花山，是當地治溝造地工程中涉及的主要土壤，模型建立所需邊坡土樣的物理力學參數取自鐘佩文等［１７］ ２０１６ 年５ 月１日至７ 月１０ 日試驗結果，見表１，土樣屬濕陷性（Ｑ３ ）黃土。

１．２ 模型建立

以萬花山治溝造地工程中兩級臺階的邊坡尺寸為依據，建立黃土邊坡數值計算模型，如圖１ 所示。單極坡高４ ｍ，邊坡高８ ｍ，單極坡比１ ∶ ０．５，兩級中間為２ ｍ寬的平臺?；冢停铮瑁颍茫铮酰欤铮恚?強度理論，利用ＦＬＡＣ３Ｄ 軟件進行數值分析。模型的底部和四周均采用固定邊界來模擬半無限場地，其中左右邊界固定ｘ方向，前后邊界固定ｙ 方向，底部固定ｚ 方向。

２ 強度折減法

２．１ 計算原理

利用ＦＬＡＣ３Ｄ 有限差分軟件，采用強度折減法計算邊坡安全系數。具體計算方法為：利用折減系數Ｆ ，對土體的黏聚力ｃ 和內摩擦角φ 進行折減，當達到失穩判據的條件時，邊坡處于臨界失穩狀態，此時原始黏聚力與折減后臨界破壞時黏聚力的比值即為安全系數。強度折減法基本公式為［１８］

ｃＦ ＝ ｃ０ ／ Ｆ （１）

ｔａｎ φＦ ＝ ｔａｎ φ０ ／ Ｆ （２）

式中： Ｆ 為折減系數， ｃ０ 為折減前的黏聚力， φ０ 為折減前的內摩擦角， ｃＦ 、φＦ 分別為邊坡到達臨界破壞狀態時的黏聚力和內摩擦角。

２．２ 邊坡失穩判據

１）計算不收斂判據。邊坡失穩破壞時，滑動區土體發生滑動并無限發展，這種狀態下，滑移面上的塑性位移會發生突變，在力的收斂標準和位移的收斂標準下，程序的求解結果均不收斂，因此可以把數值計算結果是否收斂作為判別邊坡失穩的依據。

２）塑性區貫通判據。采用強度折減法計算安全系數時，邊坡內部土體受應力作用逐漸發生塑性破壞，塑性區從坡腳發展到坡頂的過程，坡體內部會產生不同程度的塑性變形，塑性區貫通時，邊坡處于臨界失穩狀態。

３）位移突變判據。邊坡內部土體受應力作用發生剪切破壞時，邊坡部分土體會發生略微移動，邊坡達到臨界失穩狀態時，部分土體的位移量會突增，通過在邊坡內部設置位移監測點，觀察位移變化，從而判斷邊坡是否失穩。

３ 自然狀態邊坡穩定性分析

３．１ 計算不收斂判據結果分析

調用ｓｌｏｖｅ ｆｏｓ 命令使用內置模塊對邊坡模型進行分析，采用強度折減法求解安全系數，得出的剪切應變增量云圖見圖２。發生剪切應變的區域主要在上級臺階，剪切應變增量為０．００３ ０～０．０１２ ３，應變帶呈圓弧狀從坡腳逐漸向坡頂方向發展，最大剪切應變增量發生在坡腳處，為０．０１２ ３，產生局部塑性變形，從而導致上級邊坡存有潛在塑性貫通區域，此區域受土體自重的影響，有發展為滑移面的可能。依據計算不收斂判據，此時的折減系數即安全系數，數值為１．９５３ １。

３．２ 塑性區貫通判據結果分析

采用強度折剪法計算邊坡的安全系數時可以展示剪切塑性區和拉伸塑性區的狀態。通過ＦＬＡＣ３Ｄ 軟件得到的折減系數在１．７００ ０～２．０００ ０ 范圍內坡體塑性區貫通情況如圖３ 所示。對于剪切塑性區，軟件輸出結果中以ｓｈｅａｒ－ｎ 標識的區域土體正處于破壞狀態，以ｓｈｅａｒ－ｐ 標識的區域土體在折減過程中曾進入過破壞狀態，但此時已退出；對于拉伸塑性區，軟件輸出結果中相應土體狀態則分別以ｔｅｎｓｉｏｎ－ｎ 與ｔｅｎｓｉｏｎ－ｐ 標識。當處于破壞狀態的區域貫通時，貫通區上方的土體單元會因自重作用而形成滑動帶，此時邊坡處于臨界失穩狀態。

折減系數為１．７００ ０ 時上級臺階的坡腳處出現剪切塑性區，折減系數增大至１．８００ ０ 的過程中，剪切塑性區不斷向坡頂發展，剪切塑性區下方為退出破壞狀態區域、上方為處于破壞狀態區域，這表明隨著折減系數的增大，剪切應力主要集中在該區域上部，從而產生土體局部塑性變形并向上發展。折減系數增大至１．８５０ ０時，圖３（ｃ）上級臺階上部出現拉伸塑性區ｔｅｎ?ｓｉｏｎ－ｎ，表明此區域正出現拉應力集中，邊坡塑性區即將貫通。折減系數為１．８９５ ０ 時，圖３（ｆ）拉伸塑性區ｔｅｎｓｉｏｎ－ｎ 恰好貫通，該判據下臨界折減系數１．８９５ ０即安全系數。隨著折減系數的繼續增大，貫通區域擴大，當折減系數增大至１．９５０ ０ 時，邊坡靠近臨空面部分開始出現拉伸塑性區ｔｅｎｓｉｏｎ－ｐ，表明該區域曾處于破壞狀態。這說明隨著塑性貫通區域的繼續擴大，邊坡靠近臨空面部分受貫通區域拉伸破壞力以及土體自重的影響開始出現滑坡前兆，當滑坡發生時，這部分土體脫離邊坡整體，所受拉伸破壞力消失。通過觀察不同折減系數的塑性區貫通情況發現，塑性區主要在上級臺階左部，下級臺階無塑性區，穩定性好。實際情況中，底部不會出現大量剪切塑性區，原因是底部為固定約束，存在邊界效應［１９］ 。

３．３ 位移突變判據結果分析

通過ＦＬＡＣ３Ｄ 軟件中ＦＩＳＨ 語言工具在邊坡關鍵位置設置監測點，繪制不同折減系數的位移變化曲線，從而找到位移突變點，得到邊坡安全系數。

依據位移突變判據，坡體ｘ 方向位移云圖及潛在滑移面如圖４ 所示，滑移線左側靠近臨空面部分有不同程度ｘ 方向的位移，滑移線右側位移小于０．０１ ｍｍ。依據圖４ 邊坡ｘ 方向位移變化情況，在邊坡上級臺階及下級臺階設置監測點Ｐ０１～Ｐ１２（部分監測點布設位置見圖４），可以得到不同折減系數的監測點位移（部分結果見表２，負號代表ｘ 軸負方向，即臨空面方向），可繪制位移變化曲線。

通過觀察監測點的位移變化情況發現：折減系數大于１．９５３ １ 時，Ｐ０５、Ｐ０６、Ｐ１１、Ｐ１２ 監測點的位移有明顯變化，折減系數增至１．９６８ ８ 時，Ｐ０５ 監測點位移從－１２．９３４ ０ ｍｍ 增至－８９．２２０ ０ ｍｍ，Ｐ０６ 監測點位移從－１１．８８８ ０ ｍｍ 增至－８４．２８３ ０ ｍｍ，Ｐ１１ 監測點位移從－１２．００９ ０ ｍｍ 增至８２．３２６ ０ ｍｍ，Ｐ１２ 監測點位移從－１１．９２９ ０ ｍｍ 增至－８３．４７８ ０ ｍｍ；其余監測點位移在整個過程中沒有發生明顯變化，說明所在部位穩定性較好。此方法得到的安全系數選?。保梗担?１ 和１．９５７ ０的中間值１．９５５ １。

３．４ ３ 種判據所得安全系數差異分析

計算不收斂判據、塑性區貫通判據、位移突變判據得到的安全系數分別為１．９５３ １、１．８９５ ０、１．９５５ １。位移突變判據與計算不收斂判據結果接近，差值僅為０．００２ ０，塑性區貫通判據結果偏小，與位移突變判據和計算不收斂判據差值分別為０．０６０ １、０．０５８ １。

當不平衡力計算結果趨于一個非零值時可認為發生塑性變形。因此，使用ＦＬＡＣ３Ｄ 軟件對折減過程中邊坡模型不平衡力進行計算，分析邊坡失穩過程中塑性區的變化狀態，如圖５ 所示。

由圖５（ａ）和圖５（ｂ）可以看出，隨計算時步逐漸增加，不平衡力逐漸減小最后趨于零，這表明在折減系數為１．８９５ ０ 和１．９００ ０ 時，邊坡處于穩定狀態。但在塑性區貫通判據下，折減系數為１．９００ ０ 時（大于安全系數１．８９５ ０），邊坡已經失穩，此時塑性區貫通判據所得結果不準確。原因是塑性區貫通后，其附近的土體受邊界條件限制而阻止塑性區繼續發展，邊坡未達到臨界失穩狀態。由圖５（ｃ） 可知折減系數為１．９５０ ０時，不平衡力隨著計算時步增長逐漸減小，最后趨于一個非零值，說明邊坡已經發生塑性變形，且依據塑性區貫通判據，邊坡同樣為失穩狀態，此時塑性區貫通判據所得結果準確。結合圖３（ｈ）塑性區貫通情況，在邊坡模型上級臺階處，上方遍布正處于破壞狀態的拉伸塑性區ｔｅｎｓｉｏｎ－ｎ，而靠近臺階邊緣處充滿已經退出破壞狀態的拉伸塑性區ｔｅｎｓｉｏｎ－ｐ，此區域曾進入過破壞狀態，受力的影響已經變為滑動帶，邊坡失穩，圖５（ｃ）平衡力在后期略微起伏并趨于一個非零值，此時邊坡才真正發生臨界失穩。折減系數為１．９８０ 時，隨著計算時步的增長，不平衡力后期產生明顯波動，表明其不滿足不平衡力求解要求?？梢?，塑性區貫通不宜單獨作為邊坡失穩的判據，其會出現判斷結果偏小的情況，需要結合不平衡力綜合分析邊坡安全系數。綜合分析得到安全系數為１．９５０ ０，與計算不收斂判據、塑性區貫通判據所得結果的差值分別為０．００３ １ 與０．００５ １，差值明顯縮小。

由３ 種判據計算結果可以看出：數值收斂判據直接調用軟件內置模塊對邊坡進行穩定性分析，使用簡便，計算快速，但不能展示隨著折減系數的增大，坡體塑性區的變化情況。塑性區貫通判據得到的安全系數偏小，不夠真實準確，原因是塑性區貫通是邊坡失穩破壞的必要條件，而不是充分條件［２０］ ，此判據只能模擬土的強度特性，不能模擬實際應力—應變關系［２１］ ，需要結合不同折減系數的不平衡力進行綜合分析，但通過圖３ 可得到隨著折減系數的增大整個邊坡潛在的破壞區域分布以及滑動面相對位置，在這方面此判據優于其他兩個判據［２２－２３］ 。位移突變判據得到的安全系數與計算不收斂判據所得結果接近，需要在邊坡關鍵位置設置監測點，過程復雜，但能夠根據坡體位移變化情況真實反映邊坡的臨界破壞狀態［２４］ 。

４ 結論

１）基于不同失穩判據對黃土邊坡穩定性進行分析，所得安全系數由大到小排序為：位移突變判據１．９５５ １，計算不收斂判據１．９５３ １，塑性區貫通判據１．８９５ ０。位移突變判據與計算不收斂判據結果接近，差值僅為０．００２ ０，塑性區貫通判據所得結果偏小，與位移突變判據和計算不收斂判據差值分別為０．０６０ １、０．０５８ １。

２）計算不收斂判據直接調用軟件內置模塊進行邊坡穩定性分析，使用簡便，計算快速，但不展示折減過程坡體塑性區的變化情況。塑性區貫通判據所得安全系數較小，原因是塑性區貫通時邊坡未必處于臨界失穩狀態，不宜單獨作為邊坡失穩判據，需要結合不平衡力綜合評估，但此判據可以展示隨著折減系數的增大坡體的剪切和拉伸塑性區變化情況。位移突變判據分析需要在邊坡關鍵位置設置監測點，求解過程復雜，但能分析坡體位移變化狀態。

３）在實際工況中，針對類似萬花山治溝造地工程開挖的臺階型黃土邊坡，建議優先選用求解簡便、快速的計算不收斂判據進行穩定性分析、得到安全系數，并且參考其他判據得到的塑性區貫通情況與坡體位移變化狀態，以準確分析黃土邊坡穩定性狀態。

參考文獻：

［１］ 史俊濤，孔思麗，任琪．有限元強度折減法邊坡失穩判據適用性分析［Ｊ］．人民黃河，２０１４，３６（２）：１３５－１３７，１４０．

［２］ 牛越先，馬?。临|邊坡穩定性影響因素的研究［Ｊ］．水土保持研究，２０１１，１８（４）：２７３－２７６．

［３］ ＮＡＥＩＪ Ｍ，ＧＨＡＳＥＭＩ Ｈ，ＧＨＡＦＡＲＩＡＮ Ｄ，ｅｔ ａｌ．Ｅｘｐｌｉｃｉｔ ＦｉｎｉｔｅＥｌｅｍｅｎｔ Ａｎａｌｙｓｉｓ ｏｆ Ｓｌｏｐｅ Ｓｔａｂｉｌｉｔｙ ｂｙ Ｓｔｒｅｎｇｔｈ Ｒｅｄｕｃｔｉｏｎ［Ｊ］．Ｇｅｏｍｅｃｈａｎｉｃｓ ａｎｄ Ｅｎｇｉｎｅｅｒｉｎｇ，２０２１，２６（２）：１３３－１４６．

［４］ 李永亮，周國勝，李永鵬．有限元強度折減法邊坡失穩判據的適用性研究［Ｊ］．水利與建筑工程學報，２０１８，１６（５）：１２５－１２９．

［５］ 陶宇，梁偉橋，謝衛兵．有限元強度折減法分析邊坡穩定性的判據研究［Ｊ］．工業建筑，２０１９，４９（２）：１０３－１０６．

６］ ＤＯＮＧ Ｊ，ＷＡＮＧ Ｃ，ＨＵＡＮＧ Ｚ，ｅｔ ａｌ．Ｄｙｎａｍｉｃ Ｒｅｓｐｏｎｓｅ Ｃｈａｒａｃ? ｔｅｒｉｓｔｉｃｓ ａｎｄ Ｉｎｓｔａｂｉｌｉｔｙ Ｃｒｉｔｅｒｉａ ｏｆ ａ Ｓｌｏｐｅ ｗｉｔｈ ａ ＭｉｄｄｌｅＬｏｃｋｅｄ Ｓｅｇｍｅｎｔ［Ｊ］． Ｓｏｉｌ Ｄｙｎａｍｉｃｓ ａｎｄ Ｅａｒｔｈｑｕａｋｅ Ｅｎｇｉ? ｎｅｅｒｉｎｇ，２０２１，１５０：１０６８９９．

［７］ 賀俊，史俊濤，黃璐，等．強度折減法中非均質層狀土坡失穩判據適用性［Ｊ］．人民黃河，２０１４，３６（１２）：１２２－１２５．

［８］ 安光明，劉大安，鄭書彥，等．涪陵區植物油廠滑坡有限元強度折減分析［Ｊ］．水土保持研究，２００９，１６（１）：１７５－１７８．

［９］ 王其寬，梁世春，宋曉雷．基于不同失穩判據的堆載邊坡穩定性分析［Ｊ］．現代礦業，２０２１，３７（６）：１１０－１１４．

［１０］ 陳林杰，梁波．基于強度折減法的地鐵車站交叉段穩定性分析［Ｊ］．地下空間與工程學報，２０１８，１４（１）：２６０－２６５．

［１１］ 王鵬，譚明建，謝尚銪，等．貴州仁遵高速公路某特大橋岸坡穩定性分析［Ｊ］．四川建材，２０２２，４８（８）：１６７－１７０．

［１２］ 郝進鋒，劉明賓．強度折減法在單洞室穩定性分析中的失穩判據比較［Ｊ］．城市建筑，２０２２，１９（２）：１０２－１０５．

［１３］ 邵冠慧，李曉昭，趙曉豹，等．不同變形失穩模式下隧道圍巖穩定性判據研究的進展［Ｊ］．防災減災工程學報，２００９，２９（５）：５７７－５８２．

［１４］ 王健，肖鵬．富縣：治溝造地“保生態”“惠民生”［Ｊ］．西部大開發，２０１９（７）：１０６－１０７．

［１５］ ＢＲＵＳＣＨＩ Ｖ Ｍ，ＢＯＮＡＣＨＥＡ Ｊ，ＲＥＭＯＮＤＯ Ｊ，ｅｔ ａｌ．Ａｎａｌｙｓｉｓ ｏｆＧｅｏｍｏｒｐｈｉｃ ＳｙｓｔｅｍｓＲｅｓｐｏｎｓｅ ｔｏ Ｎａｔｕｒａｌ ａｎｄ Ｈｕｍａｎ Ｄｒｉｖｅｒｓｉｎ Ｎｏｒｔｈｅｒｎ Ｓｐａｉｎ： Ｉｍｐｌｉｃａｔｉｏｎｓ ｆｏｒ Ｇｌｏｂａｌ ＧｅｏｍｏｒｐｈｉｃＣｈａｎｇｅ［Ｊ］．Ｇｅｏｍｏｒｐｈｏｌｏｇｙ，２０１３，１９６（７）：２６７－２７９．

［１６］ 趙岳然．基于ＦＬＡＣ３Ｄ 的余吾礦高抽巷合理層位研究及應用［Ｊ］．能源技術與管理，２０１９，４４（６）：３８－４０．

［１７］ 鐘佩文，張慧莉，田堪良，等．持續降雨入滲對黃土邊坡穩定性的影響［Ｊ］．人民黃河，２０１８，４０（１）：７６－８１．

［１８］ ＺＩＥＮＫＩＥＷＩＣＺ Ｏ Ｃ，ＨＵＭＰＨＥＳＯＮ Ｃ，ＬＥＷＩＳ Ｒ Ｊ Ｇ．Ａｓｓｏ?ｃｉａｔｅｄ ａｎｄ Ｎｏｎ?Ａｓｓｏｃｉａｔｅｄ Ｖｉｓｃｏ?Ｐｌａｓｔｉｃｉｔｙ ａｎｄ Ｐｌａｓｔｉｃｉｔｙ ｉｎＳｏｉｌ Ｍｅｃｈａｎｉｃｓ［Ｊ］．Ｇéｏｔｅｃｈｎｉｑｕｅ，１９７５，２５（４）：６７１－６８９．

［１９］ 張愛軍，莫海鴻．有限元強度折減法中邊坡失穩位移突變判據的改進［Ｊ］．巖土力學，２０１３，３４（增刊２）：３３２－３３７．

［２０］ 李垠，程丹，蘇凱．基于有限元強度折減法的邊坡失穩判據統一性研究［Ｊ］．大地測量與地球動力學，２０１６，３６（１）：６９－７４．

［２１］ 牛志力．軟基高填路堤沉降及邊坡強度折減穩定性分析［Ｊ］．筑路機械與施工機械化，２００９，２６（４）：５４－５６．

［２２］ 趙尚毅，鄭穎人，張玉芳．極限分析有限元法講座：Ⅱ有限元強度折減法中邊坡失穩的判據探討［Ｊ］．巖土力學，２００５，２６（２）：３３２－３３６．

［２３］ 趙尚毅，時衛民，鄭穎人．邊坡穩定性分析的有限元法［Ｊ］．地下空間，２００１，７（增刊１）：４５０－４５４，５８９．

［２４］ 廖少波．基于位移突變判據的巖質邊坡動力穩定性研究［Ｊ］．地下空間與工程學報，２０１７，１３（增刊２）：９２１－９２５．

【責任編輯 張華巖】