?

考慮碳排放成本的長江鐵礦石運輸研究

2021-10-08 04:34趙曉胡鴻韜
上海海事大學學報 2021年3期
關鍵詞:碳排放算法

趙曉 胡鴻韜

摘要:為降低鋼鐵企業采購和運輸鐵礦石的綜合成本,考慮船型、航速和碳排放等因素,建立長江鐵礦石運輸的混合整數非線性規劃模型。將模型中的非線性項轉化為線性項,運用CPLEX求解器求解。針對CPLEX求解中大型規模算例的局限性,引入標準粒子群優化(particle swarm optimization, PSO)算法提高求解精度。針對標準PSO算法容易陷入局部最優的問題,提出一種基于自適應策略的改進PSO算法,動態調整慣性權重,提高算法的收斂性和全局尋優能力。通過數值實驗發現,改進后的算法在全局尋優能力和收斂能力上有一定的提高。

關鍵詞:? 多層級運輸網絡; 碳排放; 混合整數非線性規劃; 粒子群優化(PSO)算法

中圖分類號:? U695.2

文獻標志碼:? A

收稿日期: 2021-03-29

修回日期: 2021-04-26

基金項目: 國家自然科學基金(71771143)

作者簡介:

趙曉(1996—),男,四川巴中人,碩士研究生,研究方向為內河航運,(E-mail)1515992779@qq.com;

胡鴻韜(1981—),男,四川綿陽人,教授,博導,博士,研究方向為港口和航運管理與優化、供應鏈網絡設計與優化,(E-mail)hthu@shmtu.edu.cn

Meeting of the Waterborne Transport Division, World Transport Convention 2021 (WTC 2021)

Research on iron ore transportation in the Yangtze River considering carbon emission cost

ZHAO Xiaoa, HU Hongtaob

(a. Institute of Logistics Science & Engineering; b. Logistics Engineering College, Shanghai Maritime University, Shanghai 201306, China)

Abstract: In order to reduce the overall cost of steel companies purchasing and transporting iron ore, taking into account factors such as ship type, navigation speed and carbon emission, a mixed-integer nonlinear programming model of iron ore transportation in the Yangtze River is established. In the model, the nonlinear terms are converted into the linear terms, and the CPLEX solver is used to solve the model. Aiming at the limitation of CPLEX in solving medium and large scale cases, the standard particle swarm optimization (PSO) algorithm is used to improve the accuracy of solutions. For the problem that the standard PSO algorithm is easy to fall into the local optimum, an improved PSO algorithm based on adaptive strategy is proposed, which can dynamically adjust the inertia weight to improve the algorithms convergence and global optimization ability. Through numerical experiments, it is found that the improved algorithm has a certain improvement in the global optimization ability and convergence ability.

Key words: multi-level transportation network; carbon emission; mixed-integer nonlinear programming; particle swarm optimization (PSO) algorithm

0 引 言

近年來,以長江為紐帶的城市在快速發展的同時也帶動了長江航道的進一步發展。這些長江沿岸城市建有許多鋼鐵廠,需要通過長江航道來運輸冶煉鋼鐵所需的鐵礦石。從國外采購鐵礦石再運輸至國內,涉及采購價格、船型、航速等問題,因此本文研究在鐵礦石的采購和運輸過程中,如何合理選擇鐵礦石供應商、安排船型、制定航速以及考慮綠色內河航運因素,將總成本降到最低。

為實現更高質量的綠色內河航運,在考慮碳稅政策和碳排放的前提下,林貴華等[1]針對航速優化、航線配船、貨物分配問題,建立了船公司收益最大化的混合整數非線性模型。隋志超等[2]認為征收碳稅和實施碳排放交易機制是實現碳減排目標的有效手段。李晗[3]把政府征收的碳稅作為班輪公司碳減排行為的主要影響因素,并建立了考慮碳排放的集裝箱班輪航線配船優化模型。呂靖等[4]認定征收碳稅會增加集裝箱班輪公司的運營成本。李旭等[5]認為在引入碳稅政策后,需要綜合考慮單路徑碳排放約束、整個海運網絡碳排放約束和不同地區碳稅差異。苗紅云等[6]以航道和港口設計通過能力為條件,將碳排放轉化為低碳成本,并考慮減載運輸來建立內河散貨運輸網絡優化模型,并通過算例分析證明了該模型可以有效降低運輸中的碳排放。

在研究內河運輸時,學者們發現內河運輸的總成本不僅包含運營成本、碳排放成本,還包括裝卸成本、碼頭??砍杀镜?。LI等[7]選取長江沿岸的鋼鐵企業為對象,研究如何使企業運營成本最低,而YANG等[8]通過優化模型讓總運輸成本與未滿足所有需求的罰款的總和達到最小。ZHENG等[9]研究了集裝箱運輸軸輻式網絡設計,并提出一個最小混合整數線性規劃模型。YANG等[10]建立了以成本最低為目標的長江集裝箱班輪運輸網絡優化模型。KONNINGS等[11]探討了鹿特丹港集裝箱駁船裝卸作業的改進方法以達到成本最低的目標。MENG等[12]所提出的軸輻式網絡以最小化承運商的總運輸成本和樞紐運營商的運營成本為目標。MOON等[13]提出一個在軸輻式網絡中船隊部署的最小不定期船航線模型。而在內河航運的服務研究中,TAN等[14]針對船舶排程設計和航行速度最佳化問題,建立了成本最小化雙目標規劃模型。

本文在以往研究的基礎上,考慮船舶運輸鐵礦石的周期、運輸鐵礦石所使用船舶的類型以及船舶在通過不同航段時的航速、碳排放成本等因素,以鋼鐵企業的鐵礦石采購和運輸總成本最低為目標建立混合整數非線性規劃模型。將模型中的非線性項轉化為線性項,引入標準粒子群優化(particle swarm optimization, PSO)算法提高求解精度。提出一種基于自適應策略的改進PSO算法,使模型更容易實現、求解精度更高、收斂速度更快。

1 問題描述和模型

一家位于長江邊的鋼鐵企業在計劃周期t內將從國外供應商處購買的鐵礦石通過船舶運輸到國內指定港口進行卸船和轉運(見圖1)。鐵礦石不能直接從供應商處運輸至轉運港口之外的其他港口,因為受其他港口的水深限制,鐵礦石需先在規定的轉運港進行卸船并裝載到載重較小的船舶上,然后運輸至長江沿岸上游其他港口??紤]不同類型船舶的載重,決策各航段運輸所需使用船舶的類型及數量。以鋼鐵企業總運營成本為目標函數,并同時考慮與船舶航速、船舶載重等因素相關的碳排放成本。

本文所考慮的長江內河散貨運輸問題,不僅要考慮相應參數的動態性、不同類型船舶的載重量、船舶航速、碳排放成本等相關因素,而且要滿足每個接收鐵礦石的港口對鐵礦石的需求,決策每個周期內需從供應商處采購的鐵礦石數量、各接收鐵礦石的港口之間需運輸鐵礦石的數量和運輸鐵礦石所需的船舶數量。

1.1 條件假設

已知所選定鋼鐵企業的船舶資源及相關航線信息,針對案例提出以下6點假設:(1)因為本文涉及小規模算例,其周期不會太長(不超過一個月),且通常散貨船的租金一般以日為單位進行收取,所以本文周期以日為單位[15];(2)港口根據其實際地理位置排序;(3)鐵礦石只能從下游港口運輸至上游各港口,不能從上游港口運輸至下游各港口;(4)每個港口的需求是預先設定的,而且必須得到滿足;(5)船舶的靠港時間和在各港口的裝卸作業時間以及由此所產生的成本歸入船舶運輸鐵礦石的時間和成本中;(6)運輸鐵礦石的船舶都為滿載狀態。

1.2 模型建立

1.2.1 參數與變量

參數:t為周期,t=1,2,…,T;v為船型,v=1,2,…,V;s為航速,s=1,2,…,S;h為供應商,h=1,2,…,H;j為中轉港,j=1,2,…,N-1;k為目的港,k=j+1;Qv為船型為v的船舶容量大小;ths為船舶以航速s從供應商h到港口1的運輸時間;pht為在周期為t時從供應商h處采購鐵礦石的單價;chvst為在周期為t時船型為v的船舶以航速s從供應商h到港口1的滿載運輸成本;rhvst為在周期為t時船型為v的船舶以航速s從供應商h到港口1的滿載碳排放成本;Ik0為在周期開始時港口k的庫存;Ikt為在周期t末港口k的剩余庫存;dkt為在周期為t時港口k對鐵礦石的需求量;lkt為在周期為t時港口k對鐵礦石的單位存儲成本;tjks為船舶以航速s從港口j到港口k的運輸時間;cjkvst為在周期為t時船型為v的船舶以航速s從港口j到港口k的滿載運輸成本;rjkvst為在周期為t時船型為v的船舶以航速s從港口j到港口k的滿載碳排放成本;M表示足夠大的正數。

決策變量:xhvt,在周期為t時船型為v的船舶從供應商h到港口1的鐵礦石運輸數量;yhvst,在周期為t時以航速s從供應商h到港口1運輸鐵礦石所需的船型為v的船舶的數量;xjkvt,在周期為t時船型為v的船舶從港口j到港口k的鐵礦石運輸數量;yjkvst,在周期為t時以航速s從港口j到港口k運輸鐵礦石所需的船型為v的船舶的數量;zhvst是0-1變量,若在周期為t時船型為v的船舶以航速s從供應商h運輸鐵礦石到港口1,則其值為1,否則為0;zjkvst是0-1變量,若在周期為t時船型為v的船舶以航速s從港口j運輸鐵礦石到港口k,則其值為1,否則為0;Yhvst與Zjkvst是輔助決策變量,其值均大于0。

1.2.2 總成本最小化模型

本文所建立的總成本最小化模型為混合整數非線性規劃(mixed-integer nonlinear programming, MINLP)模型:

minHh=1Vv=1Tt=1phtxhvt+Hh=1Vv=1Ss=1Tt=1chvstyhvstzhvst+

N-1j=1Nk=j+1Vv=1Ss=1Tt=1cjkvstyjkvstzjkvst+

Hh=1Vv=1Ss=1Tt=1rhvstyhvstzhvst+

N-1j=1Nk=j+1Vv=1Ss=1Tt=1rjkvstyjkvstzjkvst+Nk=1Tt=1Iktlkt

(1)

s. t.

Hh=1Vv=1t-thsτ=1xhvτ+I0=

tτ=1d1τ+It+Nk=2Vv=1tτ=1x1kvτ

(2)

k-1j=1Vv=1t-tjksτ=1xjkvτ+Ik0=

tτ=1dkτ+Ikt+Nj=k+1Vv=1tτ=1xkjvτ

(3)

yhvstQv≥xhvt

(4)

yjkvstQv≥xjkvt

(5)

Vv=1Ss=1zhvst=1

(6)

Vv=1Ss=1zjkvst=1

(7)

xhvt,yhvst,xjkvt,yjkvst,Ikt≥0

(8)

zhvst,zjkvst∈{0,1}

(9)

式(1)為目標函數,其第一項為從供應商處采購鐵礦石的成本,第二項為從供應商處運輸鐵礦石到第一個港口的運輸成本,第三項為在各港口之間運輸鐵礦石的運輸成本,第四項為從供應商處運輸鐵礦石到第一個港口的碳排放成本,第五項為在各港口之間運輸鐵礦石的碳排放成本,最后一項為鐵礦石在港口的倉儲成本;式(2)和(3)分別為鐵礦石從供應商處向港口1運輸和鐵礦石在各港口之間運輸的流平衡;式(4)和(5)分別為鐵礦石從供應商處向港口1運輸和鐵礦石在各港口之間運輸時的船舶數量限制;式(6)和(7)分別為鐵礦石從供應商處向港口1運輸和鐵礦石在各港口之間運輸時船型和航速的決策;式(8)為非負變量的數值定義;式(9)為0-1變量的數值定義。

1.3 模型線性化

本文模型中包含了決策船型和航速,因此目標函數中存在非線性項。為便于求解,引入兩個輔助決策變量將目標函數中的非線性項轉化為線性項,改進后的模型如下:

minHh=1Vv=1Tt=1phtxhvt+Hh=1Vv=1Ss=1Tt=1chvstYhvst+

N-1j=1Nk=j+1Vv=1Ss=1Tt=1cjkvstZjkvst+Hh=1Vv=1Ss=1Tt=1rhvstYhvst+

N-1j=1Nk=j+1Vv=1Ss=1Tt=1rjkvstZjkvst+Nk=1Tt=1Iktlkt

(10)

s. t.

式(2)~(9)

Yhvst≥(zhvst-1)M+yhvst(11)

Zjkvst≥(zjkvst-1)M+yjkvst(12)

改進后的模型比原模型多了兩個約束(式(11)和(12))。這兩個約束建立起兩個相乘的決策變量之間的線性關系,從而將模型轉化為CPLEX能夠求解的一般線性規劃模型。

1.4 參數設置

本文研究的問題是從鋼鐵企業的角度需要解決的問題。實驗所使用的不同規模的數據是通過查閱相關航運資料、寶鋼鋼鐵企業相關資料并結合實際設置的,不同規模下的集合數據見表1。模型中的其他參數是通過利用程序代碼在合理且符合實際的范圍下隨機產生的,數據如下:(1)運輸鐵礦石的各類船舶的容量在[2,20]萬t內離散均勻產生;(2)從供應商處到港口1的運輸時間為其距離與航速的比值;(3)各供應商的鐵礦石價格在[500,900]元/t內隨機產生;(4)根據相關資料,鐵礦石運輸成本c與船舶容量Qv、船舶航速s、船舶運價p1、運價轉化系數p2存在關系c=p1Qv+p2(s+10),其中船舶在從供應商處到港口1和從港口1到其他港口的運價p1分別在[50,100]元/t和[20,50]元/t內隨機產生,船舶在從供應商處到港口1和從港口1到其他港口的運價轉換系數p2分別在[10,30]和[10,15]內隨機產生;(5)根據相關資料,碳排放成本r與船舶容量Qv、船舶航速s、單位碳排放成本k1、碳排放成本轉化系數存在關系r=k1Qv+k2(s+10),其中k1在[200,400]元/t內隨機產生,k2在[30,50]內隨機產生;(6)港口鐵礦石周期初、末的庫存均在[1 000,4 000]t內均勻產生;(7)每個港口的鐵礦石需求在[0.5,2]萬t內隨機產生;(8)每個港口的鐵礦石倉儲成本在[1,7]元/t內隨機產生;(9)從港口1到其他港口的運輸時間為其距離與航速的比值。

2 求解算法

2.1 PSO算法

PSO算法首先對粒子進行初始化,每個粒子代表一個解。在一個B維搜索空間有n個粒子,粒子i(i=1,2,…,n)的位置和速度分別為xi和vi。在每次迭代中,粒子通過跟蹤兩個極值pi和gi來更新自己,這里,pi和gi分別為個體極值和全局極值。在找到這兩個最優值后,粒子更新自己的速度和位置的方法如下:

vi=wvi+c1ri(pi-xi)+c2ri(gi-xi)(13)

xi=xi+vi(14)

式中:w為慣性權重,w∈[0.4,0.9],w值越大,說明粒子全局尋優能力越強,局部尋優能力越弱;ri為(0,1)內的隨機數;c1和c2為學習因子,通常c1=c2=2。粒子速度的最大值用vmax表示,如果粒子速度超過vmax,則取vi=vmax。

2.2 改進的PSO算法

在PSO算法中,權重w對全局尋優和局部尋優有著很大的影響,因此選擇合適的w是提高算法尋優能力的關鍵。為提高PSO算法的收斂性以及避免PSO算法陷入局部最優,引入自適應權重[15](權重w隨著每個粒子的位置、速度的更新進行自適應調整)。本文w取值與當前迭代次數d和迭代總次數D有關。

通過自適應調整權重w的PSO算法更新公式如下:

wid=wmin+(wmax-wmin)f(xid)-fdminfdavg-fdmindD,

f(xid)≤fdavg

wmax, f(xid)>fdavg

(15)

式中:wmin和wmax為預設的最小速度權重和最大速度權重,一般取wmin=0.4,wmax=0.9;f(xid)為粒子i第d次迭代時的適應度值;fdavg為第d次迭代時所有粒子的平均適應度值;fdmin為所有粒子第d次迭代時的最小適應度值,fdmin=min{f(x1d),f(x2d),…,f(xnd)}。適應度值越小,說明距離最優解越近,此時更需要局部搜索;適應度值越大,說明距離最優解越遠,此時更需要全局搜索。

改進的PSO算法流程見圖2。

2.3 算法求解過程

粒子群中每個粒子代表一個解,因為粒子是隨機生成的,所以可以保證粒子的多樣性。選取0-1變量zhvst和zjkvst進行編碼,共生成kt個zhvst和k(k-1)t/2個zjkvst作為生成的初始粒子群,經過算法迭代得到種群中最優的港口k和周期t所代表的位置,最后獲得粒子的最優位置和目標函數值。

算法中每個粒子的解P包含整數部分和小數部分,P的取值范圍為[0,vs-1]。由P/S得到的整數代表zhvst的v;P/S的余數代表zhvst中的s,其小數部分如在[0,0.5]范圍內則向下取整,如大于0.5則向上取整。在j、k、t、h分別確定的情況下,唯一的v、s也就隨之確定。

3 算 例

3.1 3種規模算例實驗

本次算例實驗平臺的CPU為Intel Core i5-9400 2.9 GHz,運行內存為8 GB,Windows 10 64位處理系統。采用Python 3.7進行編程,線性規劃求解部分調用CPLEX 12.6.2求解器進行求解。標準PSO算法參數如下:粒子群種群規模為100,迭代次數為50,權重w取wmin與wmax的均值0.6;改進PSO算法中權重w通過自適應進行調整,其余參數的設置與標準PSO算法的一致。本文設置的小規模、中規模、大規模算例分別包括3個港口、8個港口、20個港口,每個規模生成了12組數據。將每組數據輸入模型中進行求解,求解時間及結果見表2~4。由表2可知,CPLEX的求解時間較短,說明在進行小規模算例求解時,采用CPLEX進行求解是最為合適的。隨著算例規模的增大,CPLEX的求解時間也增加,當算例規模增大到一定程度時CPLEX無法在合適的時間內找出最優解,因此,中大規模算例已不適合用CPLEX進行求解。對于中大規模算例,標準PSO算法和改進的PSO算法均能在合適的時間范圍內求解出對應的最優解,而且求解出的解的質量較好。相比之下,中規模算例采用改進的PSO算法的求解效果比大規模算例的好。

3.2 參數靈敏度分析

模型中涉及的一系列參數,例如鐵礦石的采購價格、鐵礦石的倉儲成本等,分別分布在成本函數的各項成本中,其余參數均與之相關。下面對鐵礦石采購價格進行靈敏度分析,假設鐵礦石的供應商為3個,其分別位于澳大利亞、巴西、南非,且有載重量分別為60 000 t和80 000 t的兩種靈便型散貨船可供選擇,航速可選擇12 kn和15 kn。在其他參數都一樣的情況下,改變鐵礦石的采購價格,得到鋼鐵企業總成本的變化情況,見圖3。

由圖3可以看出,隨著鐵礦石采購價格的不斷提高,鋼鐵企業的總成本也在相應增加。由此可見各國的鐵礦石市場價格很大程度上也會影響鋼鐵企業的運營。模型中其他參數如各港口對鐵礦石的需求、鐵礦石的倉儲成本等發生變化產生的影響均與鐵礦石市場價格的影響一致,鋼鐵企業的總運營成本會隨著這一系列參數的變化而發生改變。

3.3 算法迭代過程對比

為驗證本文所用算法的有效性,選取中規模、大規模算例進行實驗,中(大)規模算例的相關參數設置:港口數為8(20)個、規劃周期為30 d(80 d)、船型數為6種(15種)、航速為6 kn(15 kn)以及供應商數為6個(15個)。標準PSO算法和改進的PSO算法的迭代過程對比見圖4和5。兩種算法在求解過程中逐步尋優,并在經過不同的迭代次數后收斂于不同的穩定值。由圖4和5可以知,兩種算法求解時,雖然標準PSO算法的收斂速度快于改進的PSO算法,但是改進的PSO算法所得到的結果優于標準的PSO算法。

4 結 論

針對長江內河鐵礦石運輸問題,在考慮碳排放成本的情況下,以一家位于長江邊的鋼鐵企業為研究背景,以鋼鐵企業的鐵礦石采購和運輸總成本最低為目標建立混合整數非線性規劃模型。采用標準粒子群(PSO)算法和改進的PSO算法進行求解并評價兩種算法的適用性。本文所提出的優化模型,能夠讓鋼鐵企業鐵礦石的采購和運輸總成本有明顯降低,為鋼鐵企業運營決策提供一定的參考,如選擇相對大的船型可以增加每艘船的鐵礦石運輸量,從而減少運輸次數、降低運營成本。

參考文獻:

[1]林貴華, 高潔, 李雨薇, 等. 考慮ECA和碳稅政策的班輪航線及貨物分配優化[J/OL]. 工業工程與管理: 1-12[2021-03-28]. http://kns.cnki.net/kcms/detail/31.1738.T.20200525.1834.012.html.

[2]隋志超, 呂靖. 碳稅和碳排放交易機制對班輪船隊經營的影響[J]. 中國水運, 2019, 19(2): 64-66.

[3]李晗. 考慮碳排放的集裝箱班輪航線配船研究[D]. 大連: 大連海事大學, 2019.

[4]呂靖, 毛鶴達. 硫排放控制區和碳排放限制下的班輪航線配船模型[J]. 大連海事大學學報, 2017, 43(1): 101-105. DOI: 10.16411/j.cnki.issn1006-7736.2017.01.016.

[5]李旭, 孟燕萍. 碳稅對集裝箱海運網絡運輸成本的影響[J]. 上海海事大學學報, 2019, 40(2): 12-17. DOI: 10.13340/j.jsmu.2019.02.003.

[6]苗紅云, 楊家其. 考慮低碳成本的內河散貨運輸航線網絡優化模型[J]. 武漢理工大學學報(交通科學與工程版), 2017, 41(5): 839-843. DOI: 10.3963/j.issn.2095-3844.2017.05.025.

[7]LI Feng, YANG Dong, WANG Shuaian, et al. Ship routing and scheduling problem for steel plants cluster alongside the Yangtze River[J]. Transportation Research Part E, 2019, 122: 198-210. DOI: 10.1016/j.tre.2018.12.001.

[8]YANG Dong, WANG Shuaian. Analysis of the development potential of bulk shipping network on the Yangtze River[J]. Maritime Policy & Management, 2017, 44(4): 512-523. DOI: 10.1080/03088839.2016.1275863.

[9]ZHENG Jianfeng, YANG Dong. Hub-and-spoke network design for container shipping along the Yangtze River[J]. Journal of Transport Geography, 2016, 55: 51-57. DOI: 10.1016/j.jtrangeo.2016.07.001.

[10]YANG Zhongzhen, SHI Haiping, CHEN Kang, et al. Optimization of container liner network on the Yangtze River[J]. Maritime Policy & Management, 2014, 41(1): 79-96. DOI: 10.1080/03088839.2013.780217.

[11]KONNINGS R, KREUTZBERGER E, MARA V. Major considerations in developing a hub-and-spoke network to improve the cost performance of container barge transport in the hinterland: the case of the port of Rotterdam[J]. Journal of Transport Geography, 2013, 29: 63-73. DOI: 10.1016/j.jtrangeo.2012.12.015.

[12]MENG Qiang, WANG Xinchang. Intermodal hub-and-spoke network design: incorporating multiple stakeholders and multi-type containers[J]. Transportation Research Part B, 2011, 45: 724-742. DOI: 10.1016/j.trb.2010.11.002.

[13]MOON I K, QIU Z B, WANG J H. A combined tramp ship routing, fleet deployment, and network design problem[J]. Maritime Policy & Management, 2015, 42(1): 68-91. DOI: 10.1080/03088839.2013.865847.

[14]TAN Zhijia, WANG Yadong, MENG Qiang, et al. Joint ship schedule design and sailing speed optimization for a single inland shipping service with uncertain dam transit time[J]. Transportation Science, 2018, 52(6): 1297-1588. DOI: 10.1287/trsc.2017.0808.

[15]于桂芹, 李劉東, 袁永峰. 一種結合自適應慣性權重的混合粒子群算法[J]. 哈爾濱理工大學學報, 2016, 21(3): 49-53. DOI: 10.15938/j.jhust.2016.03.010.

(編輯 賈裙平)

猜你喜歡
碳排放算法
國際主流軋差算法介紹:以CHIPS的BRA算法為例
Travellng thg World Full—time for Rree
學習算法的“三種境界”
算法框圖的補全
算法初步知識盤點
濟南市公共交通低碳發展路徑探索
新疆碳排放與經濟增長實證研究
新疆碳排放與經濟增長實證研究
寧夏碳排放與經濟增長的脫鉤關系研究
重慶市碳排放現狀及低碳發展路徑分析
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合