?

應用VSD動態酸化模型確定區域酸沉降控制目標研究

2012-12-20 09:08趙和春謝紹東北京大學環境科學與工程學院環境模擬與污染控制重點實驗室北京100871
中國環境科學 2012年3期
關鍵詞:保護率控制目標點位

趙和春,謝紹東 (北京大學環境科學與工程學院,環境模擬與污染控制重點實驗室,北京 100871)

應用VSD動態酸化模型確定區域酸沉降控制目標研究

趙和春,謝紹東*(北京大學環境科學與工程學院,環境模擬與污染控制重點實驗室,北京 100871)

為了控制區域酸沉降污染,需要制定科學的區域大氣酸沉降控制目標.本研究建立了應用VSD動態模型的多點位模擬和累積頻率分布曲線統計方法,通過模擬各酸沉降情景下某一目標年區域內土壤理化特性的變化確定其酸沉降控制目標.將此方法應用于廣州-東莞-惠州地區,在現場測量區域內25點位土壤特征的基礎上,應用VSD模型模擬各點位土壤特征對酸沉降的響應,再將模擬結果繪制成累積頻率分布曲線,據此確定該區域酸沉降控制目標.結果表明,單獨控制S沉降時,若使得該區域生態保護率達到80%,則短期和長期S沉降的控制目標分別為7.68~12g/(m2?a)和10.24~16g/(m2?a);若生態保護率為95%,短期和長期S沉降控制目標分別為5.12~8g/(m2?a)和7.68~12g/(m2?a).同時控制 S和 BC沉降時,若生態保護率為 80%,當 BC沉降為 6.4~12.8g/(m2?a)時,短期和長期 S的控制目標分別為 2.56~4g/(m2?a)和5.12~8g/(m2?a);當BC沉降為4.8~9.6g/(m2?a)時,S的控制目標為2.56~4g/(m2?a).若生態保護率為95%,當BC沉降為6.4~12.8g/(m2?a)時,短期和長期S的控制目標分別為0.64~1g/(m2?a)和5.12~8g/(m2?a);當BC沉降為4.8~9.6g/(m2?a)時,短期和長期S的控制目標分別為0.64~1g/(m2?a)和2.56~4g/(m2?a);當BC沉降量降至2~4g/(m2?a),則80%和95%生態保護率下的S控制目標均為0.64~1g/(m2?a).

酸沉降;區域控制目標;VSD;累積頻率分布

酸沉降仍是目前的主要環境問題,對土壤、森林、湖泊、動植物和建筑等產生不同程度的影響[1].近年來,動態模型廣泛應用于模擬酸沉降對生態系統的影響.1990年起歐洲影響合作中心(CCE)每2年應用動態模型模擬歐洲地區不同生態系統的恢復過程,為歐洲酸沉降控制提供科學依據[2-6].中國在 20世紀末就有學者應用動態模型研究西南和東南部分省份的酸沉降,分析其酸化狀態和趨勢[7-10],但這些研究多局限于單個站點或單個生態系統的模擬,應用范圍較小,對于區域酸沉降控制意義不大.因此,需要開發應用動態模型確定區域酸沉降控制目標的方法.本研究應用簡單動態模型(VSD)聯合模擬區域內多個點位,并應用累積頻率分布(CFD)的方法統計模擬結果,以確定區域酸沉降控制目標.將此方法應用于酸沉降嚴重的廣州-東莞-惠州地區,基于該地區25個點位的土壤采樣和現場調查,應用VSD模型模擬不同酸沉降控制情景下區域內各點位土壤性質的變化,進而確定區域酸沉降控制目標.

1 研究方法

1.1 VSD模型的基本原理

VSD模型是一個單層酸沉降動態模型,結構簡單,輸入參數相對較少[13].模型包括土壤中主要物理化學過程,如土壤礦物的化學風化、土壤陽離子交換、水鋁礦溶解、土壤有機酸溶解、植被吸收和氮固定等.由電荷平衡、化學平衡和質量平衡關系式構成,忽略SO42-的吸附以及Al的絡合反應形態,假設N完全硝化等.通過輸入模擬點位的土壤、植被、沉降等基本參數可得到土壤主要理化性質.其主要平衡關系如下[11].

1.1.1 電荷平衡 VSD模型中的電荷平衡方程如下式:

式中: BC=Ca+Mg+K;Org表示有機酸.由于假設完全硝化,因此[NH4+]=0

1.1.2 化學平衡 VSD模型的化學平衡主要包括水鋁礦溶解平衡、HCO3-的解離平衡、有機酸溶解平衡以及陽離子交換平衡.它們的平衡關系式如下:

水鋁礦平衡:式中: KAlox為解離常數;α為 pH-pAl關系斜率(α≤3).

HCO3

-解離平衡:

式中: K1為一級解離常數;KH為亨利常數;PCO2為土壤中二氧化碳分壓(101325 Pa).

有機酸溶解平衡:

式中: DOC為土壤中溶解的有機碳濃度(molC/m3);m 為有機碳中活性組分的濃度(mol/molC);Korg為解離常數,可由(5)式計算得到:

陽離子交換平衡主要由以下2個方程決定:

Gaines-Thomas方程:

Gapon方程:

式中: EX為離子X在土壤交換基上的平衡濃度; KAlBC和KHBC為Al-H以及H-BC交換反應的交換系數.

1.1.3 質量平衡 VSD模型中各離子的質量平衡方程如下式:

式中: Xtot為單位面積土壤中 X離子的總濃度, eq/m2; Xin為X離子的年總輸入量, eq/(m2·a),包括大氣沉降和植被利用等;Q為徑流量, m/a.

式中,

式中: θ為土壤含水率, m3/m3; z為土壤層厚度, m.

式中: EBC為 BC在土壤交換基上的平衡濃度;ρ為土壤容重,g/cm3;CEC 為離子交換容量, meq/kg.

1.2 累積頻率分布曲線的建立與酸沉降控制目標的確定

首先將研究區域劃分成一定大小的網格,在每個網格內選取具有代表性的點位,如代表性土壤或生態系統類型,隨后將這些點位的基本參數輸入VSD模型,以模擬不同酸沉降情景下和不同目標年時土壤主要理化性質的變化.假設區域內有n個模擬點位,動態模型模擬得到某一目標年下 n個點位的土壤某理化性質為 x1,x2,……, xn,x1<x2<……<xn,于是某一點位理化性質 xk的累積頻率Pk(0<Pk<1)定義為x1~xn中小于xk的x所占的比例,由此就可得到各點位的累積頻率P1,P2,……,Pn,以累積頻率P為x軸、理化性質x為 y軸作圖,就可得到各點位的該理化性質在某一目標年下的累積頻率分布曲線,如圖1所示.

圖1 累積頻率分布曲線Fig.1 Cumulative Frequency distribution curve

從圖1累積頻率分布曲線可看出,在土壤理化性質達到某閾值時模擬點位所占的比例,假設該理化性質閾值為xk,即當x<xk時土壤即處于酸化狀態,xk對應的累積頻率為 Pk%,也就是說有Pk%模擬點位的土壤的理化性質在閾值以下,即此時區域內有 Pk%的土壤將受到酸沉降的危害.因此,基于累積頻率分布曲線就可獲得某目標年內區域內土壤理化性質的恢復情況,進而根據設定的不同生態保護率確定區域酸沉降控制目標.

1.3 研究區域及輸入參數的確定

選取的研究區域為廣州-東莞-惠州一帶,如圖2所示.該區域為珠三角地區S沉降量最大的地區,區域內S沉降量多在12.8g/(m2?a)以上,部分地區可達到 20g/(m2?a).該區域鹽基陽離子沉降量也很高,多在 8g/(m2?a)以上,個別地區可達到16g/(m2?a)以上.區域內植被主要為闊葉林、針闊林、馬尾松和灌木等;土壤主要為自然酸性土壤,對酸沉降比較敏感.VSD模型的主要輸入參數如表1所示,包括土壤、植被、沉降等基本參數,如土壤礦物風化速率,植被對N和BC的吸收速率,主要離子的沉降速率及其他基本參數等.

圖2 研究區域及模擬點位分布Fig.2 The simulation region and sites, Guangzhou-Dongguan-Huizhou

表1 VSD動態模型的輸入參數Table 1 Input parameters of VSD model

將研究區域劃分成 0.2°×0.2°的網格,應用網格布點法采集土壤樣品,共采集25個點75個土壤樣品,用 X射線衍射法(XRD)測定土壤礦物組成.土壤風化速率應用PROFILE模型計算,其他土壤理化特性數據來自文獻調研和土壤普查[12-14],表2列出了部分點位的土壤風化速率和理化參數.

表2 區域內主要自然土壤土壤參數(部分點位)Table 2 The soil parameters at part of sites in the region

本研究收集了區域內不同植被類型的生產力以及優勢物種化學元素組成的資料[15-16],在此基礎上根據公式(11)計算得到區域內主要植被的氮和鹽基陽離子吸收速率,并根據已有研究結果總結出主要植被類型的干沉降因子[17-18],主要植被類型的相關參數見表3.

表3 區域內主要植被氮和鹽基陽離子吸收速率和干沉降因子Table 3 The uptake rates and dry deposition factor of major vegetables in the region

式中: Kt和Kb分別是干和枝的凈生產力; Xt和Xb分別為元素在干和枝中的含量.

各離子的總沉降量是干沉降量和濕沉降量的總和,應用式(12)計算得到:

式中: Xdep為某離子的總沉降量, kmol/(hm2?a);[X]為雨水中的組分濃度,mol/L;P為年降雨量,m/a; fDD為干沉降因子.雨水中各離子的濃度和降雨量來自區域內環境監測部門監測獲得的酸雨常規監測資料,見表4.

通過計算得到研究區域內當前 S沉降量為12.8~20g/(m2?a),BC沉降量為 8~16g/(m2?a).對于模擬點位各離子的歷史沉降數據,根據文獻[19]假設主要離子的歷史變化趨勢與 S沉降變化趨勢一致;未來沉降情景以 2010年為基準年和2020年為目標年,設定兩種控制情景:(1)單獨控制S沉降,設定不控制S沉降和基準年上分別削減 20%、40%、60%、80%S沉降幾種情景;(2)同時控制S和BC沉降,設定在基準年上分別削減20%、40%和 75%BC沉降時,分別同時削減20%、40%、60%、80%S和95%的S沉降.

表4 廣州-東莞-惠州主要離子濕沉降量年均值[kmol/(hm2?a)]Table 4 Average value of deposition of main ions in the Guangzhou-Dongguan-Huizhou[kmol/(hm2?a)]

2 結果與討論

2.1 模型校驗

以位于廣州從化流溪河地區的一個點位說明,其土壤為花崗巖赤紅壤,植被為闊葉林,是區域內代表性的土壤和植被類型.模型校驗以1900年為起點,認為該年份時尚未有人為污染的影響,現狀年為2002年,應用VSD模型中的模型校驗功能反復計算,主要的基本參數均來自于該地區已有研究實測值[20],通過不斷調整模型參數(表 5),計算出 2002年該地區土壤化學性質,結果列于表 6.比較表 6中模擬值與實測值可看出,模擬結果與當地實測值基本吻合,說明可應用 VSD模型和這些校驗后的參數模擬該區域內未來不同酸沉降下土壤理化性質的變化.

表5 廣州花崗巖赤紅壤VSD模型校驗參數Table 5 The calibration parameters of VSD model of red granite soil in Guangzhou

表6 廣州花崗巖赤紅壤土壤水性質模擬與實測對比(mmol/m3)Table 6 The simulation and measurement soil solution parameters of red granite soil in Guangzhou(mmol/m3)

2.2 區域酸沉降控制目標的確定

應用VSD動態模型模擬區域內25個點位的土壤理化特性,在(1)單獨控制S沉降和(2)同時控制S和BC沉降情景下模擬1900~2100年間研究區域內土壤鹽基飽和度(BS)的變化,再統計得到目標年下土壤 BS的累積頻率分布圖,據此確定區域酸沉降控制目標.

2.2.1 單獨控制S沉降 單獨控制S沉降時,得到模擬區域內2020年和2100年不同S沉降控制目標下土壤BS的累積頻率分布曲線,見圖3.

通常以BS<0.2作為土壤理化性質惡化的標志.由圖3可看出,2020年如維持當前沉降量不變,將有接近50%的點的BS<0.2,即區域內50%的土壤會受到酸沉降威脅;當削減20%S時,該比例將降到約35%;當削減比例提高到40%、60%和80%時,區域內將分別有20%、15%和5%土壤會受到酸沉降威脅.由圖3還可見,2100年維持當前沉降量不變時BS<0.2的點仍接近50%;當削減20%、40%、60%和80%的S沉降時,BS<0.2所占的比例分別下降到 20%、5%、0%和 0%.由此可見,削減S沉降對土壤性質的恢復效果比較明顯,且隨時間增長效果逐漸顯現,說明土壤理化性質的恢復是一個長期的響應過程.若要使得區域內80%以上的土壤不受到酸沉降的影響,短期內 S沉降的削減比例需達到40%,長期需達到20%,結合當前區域的S沉降量,在此生態保護率下該區域短期和長期的 S沉降量需控制為 7.68~12g/(m2?a)和 10.24~16g/(m2?a),該沉降量則可作為該區域內80%生態保護率下的S沉降目標負荷.若生態保護率提高到 95%,短期內需要削減80%的S沉降,長期需要達到40%,相應的S沉降量需要控制到5.12~8g/(m2?a)和7.68~12g/(m2?a),該沉降量便可作為該區域內 95%生態保護率下的S沉降目標負荷.該模擬結果與目前該區域的S沉降的臨界負荷的研究結果基本一致[21].

圖3 單獨控制S沉降時,2020年和2100年土壤BS累積頻率分布曲線Fig.3 The BS CFD curve in2020 and 2100 when S deposition is reduced alone

2.2.2 同時控制S和BC沉降 同時控制S和BC沉降時,應用VSD模型模擬得到圖4和圖5所示模擬區域在2020年和2100年不同S和BC控制方案下土壤 BS的累積頻率分布圖.模擬了在削減20%、40%以及75%BC沉降的同時削減S沉降的各種情景,由于難以獲得模擬區域BC的背景濃度數據,只好依據目前已有研究近似處理

[21-22],亦即當BC削減75%時可近似認為是其背景的濃度水平.

圖4 同時控制S和BC沉降時,2020年模擬區域內土壤BS累積頻率分布對比Fig.4 The BS CFD of the soil in the region in 2020 when S and BC deposition is reduced simultaneously

由圖4和圖5可看出,BC的削減可以明顯減弱削減酸性污染物的效果.特別是當BC的削減比例大于或等于S的削減比例時,模擬區域內幾乎所有點位的土壤BS均小于0.2,即整個區域面臨酸化的危險,說明該區域土壤本身的酸沉降承受能力較差,主要堿度來源于大氣BC的沉降.

圖5 同時控制S和BC沉降時,2100年模擬區域內土壤BS累積頻率分布對比Fig.5 The BS CFD of the soil in the region in 2100 whenS and BC deposition is reduced simultaneously

由圖4和圖5的累積頻率分布曲線可得到各種S和BC協同控制比例下模擬點位中BS<0.2的比例,進而得到不同生態保護率下酸沉降控制目標,結果見表7.

由表7可看出,隨著BC沉降量逐漸降低,相應S沉降控制目標也愈加嚴格.特別是當BC沉降量削減 75%時,相應 S沉降控制目標將降至1g/(m2?a)以下.因此,當區域內協同控制 S和 BC時,需根據保護率合理選取酸沉降控制目標,以達到最優控制效果.

3 結論

3.1 應用 VSD動態模型聯合模擬多點位生態系統理化特性,并應用累積頻率分布曲線的方法統計其模擬結果,可確定區域尺度酸沉降控制目標.通過計算各模擬點位在各酸沉降情景下某一目標年時土壤理化性質的累積頻率,并繪制累積頻率曲線,可得到當該理化性質達到某一閾值時區域內點位所占的比例,進而確定區域尺度的酸沉降控制目標.

3.2 將上述方法應用于廣州-東莞-惠州地區,模擬計算區域內 25個點位.模擬結果表明,在單獨削減S沉降時,隨著S沉降削減比例的增加區域內土壤的恢復情況變好.當生態保護率為 80%時,該區域短期和長期的S沉降量的控制目標分別為7.68~12g/(m2?a)和10.24~16g/(m2?a).若生態保護率提高到95%,相應的S沉降量的控制目標分別為5.12~8g/(m2?a)和7.68~12g/(m2?a).以上沉降控制目標可以作為該區域內80%和 95%生態保護率下的S沉降目標負荷.

3.3 控制S沉降的同時控制大氣BC沉降,會比較明顯的降低削減S的效果.當BC的削減比例大于或等于S沉降削減比例時,幾乎整個區域內均面臨酸沉降危害,說明區域內土壤本身承受酸沉降的能力較低,大氣 BC沉降是區域內重要的堿度源.當BC沉降為6.4~12.8g/(m2?a)時,80%生態保護率下短期和長期 S的控制目標為2.56~4g/(m2?a)和5.12~8g/(m2?a),95%生態保護率下短期和長期S的控制目標為0.64~1g/(m2.a)和5.12~8g/(m2?a);當 BC沉降為 4.8~9.6g/(m2?a)時, 80%生態保護率下 S的控制目標為 2.56~4g/ (m2?a), 95%生態保護率短期和長期 S的控制目標為0.64~1g/(m2?a)和2.56~4g/(m2?a).而當BC沉降量降至2~4g/(m2?a),則80%和95%生態保護率 下的S控制目標均為0.64~1g/(m2?a).

表7 廣州-東莞-惠州區域不同生態保護率下和不同鹽基陽離子沉降下硫沉降控制目標Table 7 The acid deposition control target under different ecosystem protection rate and deposition rate of BC in Guangzhou-Donguan-Huizhou

[1] 馮宗煒,曹洪法,周修萍,等.酸沉降對生態環境的影響及其生態恢復 [M]. 北京:中國環境科學出版社, 1999.

[2] Hettelingh J P, Downing R J, De De Smet P A M Mapping critical loads for Europe [R]. CCE Technical Report No.1, Bilthoven, the Netherlands, 1991.

[3] Slootweg J, Posch M, Hettelingh J P. Critical loads of nitrogen and dynamic modeling [R]. CCE Progress Report 2007, Bilthoven, the Netherlands, 2007.

[4] Ferrier R C, Wright R F, Cosby B J, et al. Application of the MAGIC model to the Norway spruce stand at Solling, Germany [J]. Ecological Modelling, 1995,83:77-84.

[5] Sverdulp H, Warfvinge P, Blake L, et al. Modelling recent and historic soil data from the Rothamsted experimented station, UK using SAFE [J]. Agriculture Ecosystems and Environment, 1995, 53:161-177.

[6] Langan S, Fransson L, Vanguelova E. Dynamic modeling of the response of UK forest soils to change in acid deposition using the SAFE model [J]. Science of the Total Environment, 2009,407: 5605-5619.

[7] Zhao D W, Seip H M. Assessing effects of acid deposition in southwestern China using the MAGIC model [J]. Water, Air and Soil Pollution, 1991,60:83-97.

[8] 趙殿五,張曉山,熊際翎.應用MAGIC模式確定酸沉降臨界負荷[J]. 中國環境科學, 1992,12(2):93-97.

[9] 謝紹東,郝吉明,周中平,等.柳州地區酸沉降臨界負荷的確定[J]. 環境科學, 1996,17(5):1-4.

[10] An J L, Huang M Y. Long-term soil acidification model (LTSAM) development and application for analyzing soil responses to acidic deposition [J]. Water, Air and Soil Pollution, 1999,110: 255-272.

[11] Alveteg M, Walse C, Sverdrup H. Evaluating simplifications used in regional applications of the SAFE and MAKEDEP models [J]. Ecological Modelling, 1998,107:265-277.

[12] 張遠東.珠江三角洲地區酸雨污染簡析 [J]. 環境科學研究, 1999,12(3):31-34.

[13] Posch M, Reinds G J. A very simple dynamic soil acidification model for scenarios analysis and target loads calculation [J]. Environmental Modelling and Software, 2009,24:329-340.

[14] 段 雷.中國酸沉降臨界負荷區劃研究 [D]. 北京:清華大學, 2000.

[15] 廣東省土壤普查辦公室.廣東土壤 [M]. 北京:科學出版社, 1993.

[16] 解憲麗,孫 波,周慧珍,等.中國土壤有機碳密度和儲量的估算與空間分布分析 [J]. 土壤學報, 2004,41(1):35-42.

[17] 楊 昆,管東生.珠江三角洲森林的生物量和生產力研究 [J].生態環境, 2006,15(1):84-88.

[18] 李燕燕.馬尾松-闊葉樹混交林生物量及礦質養分的研究 [D].福州:福建農林大學, 2000.

[19] 劉菊秀,溫達志,周國逸.廣東鶴山酸雨地區針葉林與闊葉林降水化學特征 [J]. 中國環境科學, 2000,20(3):198-202.

[20] 陳 勇.珠江三角洲城市森林群落對降水及土壤侵蝕的影響的研究 [D]. ???華南熱帶農業大學, 2004.

[21] Allen S L, Janja D H, Rudolf D H. Estimating historical anthropogenic global sulfur emission patterns for period 1850-1990 [J]. Atmospheric Environment, 1999,33:3435-3444.

[22] Larssen T, 何 毅,湯大鋼.中國酸沉降綜合影響觀測研究 [R]. 2003年研究報告.中國環境科學研究院,挪威水研究所, 2003.

[23] 張懿華.珠江三角洲酸沉降臨界負荷區劃及沉降控制選擇研究[D]. 北京:北京大學, 2009.

[24] Zhao Y, Duan L, Thorjorn L, et al. Simultaneous assessment of deposition effects of base cations, sulfur, and nitrogen using an extended critical load function for acidification [J]. Environmental Science and Technology, 2007,41(6):1815-1820.

Determining regional control targets of acid deposition using VSD model.

ZHAO He-chun, XIE Shao-dong*(State Key Joint Laboratory of Environmental Simulation and Pollution Control, College of Environmental Sciences and Engineering, Peking University, Beijing 100871, China). China Environmental Science, 2012,32(3):411~418

Setting scientific acid deposition control targets was necessary to control regional acid deposition pollution. A method to set such control targets by analyzing soil acidity at selected sites under different acid deposition scenarios in the target years, combining the VSD model and cumulative frequency distribution analysis, was proposed in this study. This method was applied to the Guangzhou-Dongguan-Huizhou region. The soil acidity in 25 sites was measured and simulated by using VSD model under different acid deposition scenario, and the results were described by using cumulative frequency distribution curve. The S deposition target under different scenarios was given in the results. If S deposition was controlled solely, the short-term and long-term S deposition control targets should be 7.68~12g/(m2?a) and 10.24~16g/(m2?a), respectively, to guarantee that 80% of the ecosystem was protected, and the short-term and long-term S deposition control targets should be 5.12~8g/(m2?a) and 7.68~12g/(m2?a) respectively, to guarantee that 95% of the ecosystem was protected. If S and BC deposition were controlled simultaneously, when BC deposition was 6.4~12.8g/(m2?a), the short-term and long-term S deposition control targets should be 2.56~4g/(m2?a) and 5.12~8g/(m2?a), respectively, when BC deposition was 4.8~9.6g/(m2?a), the S deposition control target should be 2.56~4g/(m2?a), to guarantee that 80% of the ecosystem was protected; when BC deposition was 6.4~12.8g/(m2?a), the short-term and long-term S deposition control targets should be 0.64~1g/(m2?a) and 5.12~8g/(m2?a), when BC deposition was 4.8~9.6g/(m2?a), the short-term and long-term S deposition control targets should be 0.64~1g/(m2?a) and 2.56~4g/(m2?a), to guarantee that 95% of the ecosystem was protected. When BC deposition was reduced to 2~4g/(m2?a), S deposition should be controlled to 0.64~1g/(m2?a) to make sure that 80% and 95% of ecosystem was protected. Restoration measures should be taken at the same time.

acid deposition;regional control target;VSD;cumulative frequency distribution

X517

A

1000-6923(2012)03-0411-08

2011-04-27

國家“863”項目(2006AA06A306)

* 責任作者, 教授, sdxie@pku.edu.cn

趙和春(1985-),男,遼寧沈陽人,北京大學環境科學與工程學院碩士研究生,主要從事酸沉降控制研究等方面的研究.發表論文1篇.

猜你喜歡
保護率控制目標點位
科興新冠疫苗對印尼醫護保護率94%
基于結構光視覺的鉆孔點位法矢檢測技術研究
提高林場造林成活率和保存率的措施
牛脾轉移因子聯合核酸疫苗對虹鱒魚IHN病的保護率影響研究
血糖控制目標應“因人而異”
兒童接種水痘疫苗效果及影響因素探析
大盤仍在強烈下跌趨勢中
公路路基工程施工成本控制和管理
基于空間網格的機器人工作點位姿標定方法
中國首份行政事業單位內部控制實施情況白皮書(三)
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合