?

兩方零和馬爾科夫博弈策略梯度算法及收斂性分析

2024-03-12 12:48王卓李永強馮宇馮遠靜
浙江大學學報(工學版) 2024年3期
關鍵詞:馬爾科夫納什梯度

王卓,李永強,馮宇,馮遠靜

(浙江工業大學 信息工程學院,浙江 杭州 310000)

兩方零和博弈(two-player zero-sum game)是博弈論中最基本的,同時也是最常見的博弈模型之一.其特點是一個玩家獲得的收益必然等于另一個玩家所失去的收益,即參與博弈的2個玩家的收益之和為零.此類博弈問題近年來引起了經濟、交通、軍事等各個領域的高度重視,其中一些應用效果取得了里程碑式的突破,如針對完全信息圍棋博弈的AlphaGo[1],非完全信息回合制撲克博弈的Pluribus[2].

納什均衡 (Nash equilibrium)[3]是博弈論中的一個重要概念.在均衡狀態下,任何玩家都不能僅通過更改自己的策略來獲得更多的收益[4].現有研究中,通常根據博弈特征與信息完備性的差異將兩方零和博弈分為擴展式博弈(extensive-form game)和馬爾科夫博弈(Markov game).馬爾科夫博弈(Markov games)被視為馬爾科夫決策過程(Markov decision process,MDP)在博弈場景中的推廣,適用于描述完全信息、同時移動博弈.其特點為在博弈過程中,下一步的狀態只依賴于當前狀態,而不依賴于之前的狀態.從原理上,多智能體博弈強化學習方法可以分為2類: 基于值函數方法和基于策略方法.對于兩方零和馬爾科夫博弈,基于值函數方法旨在找到最優值函數,并由其獲取相應的平穩納什均衡策略.在狀態轉移概率未知的情況下,Littman[5]通過歷史數據近似估計Bellman算子,提出Minimax-Q學習,將單智能體強化學習的Q學習算法[6]擴展到兩方零和馬爾科夫博弈.對于大規模問題,為了解決維度災難,Fan等[7-9]將函數逼近思想引入Minimax-Q,利用函數逼近對值函數進行近似,但基于值函數的方法收斂速度較慢.神經網絡作為強大的非線性函數逼近器,具有自動學習特征的能力,使用神經網絡對值函數和策略進行近似可以從原始數據中學習到適合問題的特征表示,從而有效地解決大規模問題中維度災難的問題.

基于策略的策略梯度方法在MDP中具有最先進的性能[10-11].在策略參數化的設定下,尋找零和馬爾科夫博弈的納什均衡解通常表述為非凸非凹的極大極小優化問題[12-13].Daskalakis等[14]探討在多智能體競爭性環境中,如何使用獨立策略梯度方法(independent policy gradient,IPG)來訓練多個智能體的策略,以實現博弈最優解的學習.Nouiehed等[15]提出的嵌套梯度下降方法(nested gradient descent method)通過交替地對每個玩家的策略進行梯度下降來逐步逼近博弈的納什均衡點,與IPG(該方法中每個智能體的策略被看作是獨立的)相類似,兩者本質都是單智能體強化學習方法.在納什均衡的基礎上,Fiez等[16]提出雙時間尺度算法(two-timescale algorithm),其中追隨者使用梯度博弈更新規則,而不是精確的最佳響應策略,該算法已被證明收斂于Stackelberg均衡.然而,該均衡只對應于零和博弈的單邊均衡解,并且無法分析同時移動博弈的收斂性.Perolat等[17]針對同時移動博弈,將演員-評論家網絡(Actor-Critic network)與虛擬博弈相結合,提出基于樣本平均策略的算法并解決了虛擬博弈算法中的收斂問題.然而,在每輪迭代時都須求解一方的最佳響應策略,因此對計算資源的消耗也是巨大的.

針對上述方法的不足,本研究基于額外梯度法[18]提出雙方玩家同時進行策略梯度更新的兩方零和馬爾科夫博弈策略優化算法.1)基于參數化策略,提出馬爾科夫博弈策略梯度定理,并進行了嚴格證明.2)將回報作為狀態動作價值函數的無偏估計,基于額外梯度法提出馬爾科夫博弈的策略梯度算法,并證明該算法可以收斂到近似納什均衡.3)針對大規模的Oshi-Zumo游戲,采用神經網絡作為參數化策略,驗證了算法的有效性.

1 兩方零和馬爾科夫博弈

1.1 問題描述

兩方零和馬爾科夫博弈問題可以用一個五元組來描述[19],定義如下:

式中:O表示狀態空間;U表示玩家1的動作空間;W表示玩家2的 動作空間;p(s′,r|s,a,b):O×U×W→Δ(O×R) 表示在任意聯合動作 (a,b)∈U×W下,從任意狀態s∈O轉移到狀態s′∈O,且玩家1獲得獎勵r∈R ,玩家2獲得獎勵 -r的概率;γ為折扣因子,γ∈(0,1.0].玩家1、2的動作分別通過隨機策略 π(a|s):O→Δ(U) 和μ (b|s):O→Δ(W) 來選擇.每個時刻t,玩家1、2根據當前環境的狀態St分別同時選 擇動作At和Bt,St∈O,At~π(·|St),Bt~μ(·|St);聯合動作送入環境中執行,環境按照概率分布p進行狀態轉移并給出獎勵,即St+1,Rt~p(·,·|St,At,Bt),玩家1的獎勵為Rt,玩家2的獎勵為 -Rt.回報,即長期累積獎勵,定義為

式中:T為終止時刻.

給定聯 合策略 (π,μ) ,狀態價 值函數Vπ,μ(s):O→R 和狀態動作價值函數Qπ,μ(s,a,b):O×U×W→R都定義 為期望回報,即

在上述設定下,玩家1的目標是最大化期望回報,而玩家2的目標是最小化期望回報.在馬爾科夫博弈中,每個玩家的回報不僅僅與自己的策略有關,還取決于其他玩家的策略,因此馬爾科夫博弈最優策略的定義與MDP最優策略的定義有本質差異.馬爾科夫博弈的最優策略通常用納什均衡來描述,最佳響應策略[20]和納什均衡策略[21]的定義如下.

定義1 兩方零和博弈最佳響應策略.給定玩家2策略 μ,如果玩家1的策略沒有比πb更好的,那么稱 πb為最佳響應策略,即

相反地,給定玩家1的策略 π,玩家2的最佳響應策略 μb滿足

定義2 兩方零和博弈納什均衡策略.如果2個玩家的策略 π*和 μ*都是對方的最佳響應策略,則 (π*,μ*) 為納什均衡聯合策略,即

博弈雙方處于納什均衡中,意味著任何玩家都不能僅通過更改自己的策略來獲得更多的獎勵,因此納什均衡可以被視為“對最佳反應的最佳反應”.對于整個博弈過程來說,納什均衡聯合策略即為兩方零和馬爾科夫博弈的最優策略.在兩方零和博弈中,始終存在著納什均衡解,其等價于一個最大最小問題的優化解[22],即

因此,將納什均衡的求解過程轉化為上述最大最小問題的優化過程.

1.2 貝爾曼方程

為了描述不同時刻以及不同價值函數之間的關系,為后續證明提供理論依據,提出以下引理.

引理1對于兩方零和博弈,給定聯合策略(π,μ),當前時刻狀態價值函數與當前時刻狀態動作價值函數之間關系的貝爾曼方程為

當前時刻狀態動作價值函數與下一時刻狀態價值函數之間關系的貝爾曼方程為

證明:將式(3)所示的狀態價值函數的定義展開,結合式(4)可以得到

將回報Rt的定義(式(2))代入式(4),展開可得

對式(12)中的第2項求期望就等于下一時刻的狀態價值函數:

2 馬爾科夫博弈的策略梯度

為了實現納什均衡策略的求解,考慮參數化的策略 πθ(a|s) 和μ?(b|s) ,其中θ和?為可調參數,玩家的目標是最大化自己的期望回報,因此考慮指標函數

式中:ρ0為初始狀態的概率分布.那么,博弈問題(式(1))可以轉化為如下最大最小問題:

為了利用梯度法求解式(15),須得到策略梯度,即J(θ,?) 分別關于參數θ 和? 的梯度.將MDP的策略梯度擴展到兩方零和馬爾科夫博弈問題(式(1)),定理如下.

定理1對于兩方零和馬爾科夫博弈問題(式(1))和參數化的聯合隨機策略 (πθ,μ?),式(14)定義的指標函數J(θ,?) 分別關于參數θ和?的梯度為

證明:首先考慮玩家1的策略參數 θ,對式(9)兩邊求導得到

通過上述推導,可以計算式(14)中指標函數的梯度:??J(θ,?)的證明同理,定理1證畢.

由定理1可以得到兩方零和馬爾科夫博弈中2個玩家策略梯度,從而可以利用梯度方法求解該最大最小優化問題(式(15)).由狀態動作價值函數的定義(式(4)),可知回報Gt是的無偏估計,因此給出如下推論.

推論1對于兩方零和馬爾科夫博弈問題(式(1))和參數化的聯合隨機策略 (πθ,μ?),式(14)定義的指標函數J(θ,?) 關于參數θ和?的梯度分別為

證明:以玩家 1的梯度項為例,由式(4)、(16) 可知:

玩家2梯度項的證明同理,推論1證畢.

3 馬爾科夫博弈策略梯度算法及收斂性

3.1 求解最大最小優化問題的梯度方法

考慮凸凹函數J(θ,?)=θ? ,其中θ,?∈R,顯然?θJ(θ,?)=?,??J(θ,?)=θ ,且 (0,0) 為該優 化問題的最優解.若使用基于梯度的方法求解式(15),最簡單的方法為同時梯度上升下降法(simultaneous gradient ascent-descent),即

式中:[·;·]表示列向拼接,η為學習率.由式(28)以及遞歸思想,可以計算[θk+1;?k+1]二范數的平方:

圖1 不同更新規則的探索軌跡Fig.1 Exploration trajectories for different update rules

第2種基于梯度的方法為交替梯度上升下降(alternating gradient ascent-descent),參數 θ、? 的更新是按 順序進行的,即 θk+1=θk+η?k,?k+1=?kηθk+1,其中第2個更新式的參數?k+1在計算梯度時使用了其他參數(θk+1)的更新.該更新過程呈平穩性,不收斂于 (0,0) 并陷入極限環[18],如圖1(b)所示.

第3種改進的方法為近端點法(proximal point method),更新形式為 θk+1=θk+η?k+1,?k+1=?kηθk+1.該方法可收斂于該問題最優解,如 圖1(c)所示.下面給出簡單證明.

令θk+1=θk+η(?k-ηθk+1),?k+1=?k-η(θk+η?k+1),則有因此

使用額外梯度法解決該優化問題,相較于同時梯度上升下降、交替梯度上升下降具備良好的收斂性,相較于近端點法具有良好的收斂速度,如圖1(d)所示.額外梯度法的核心思想在于通過一個簡單的梯度更新步驟逼近即

相比梯度上升下降法,每次迭代,額外梯度法增加一步外推(extrapolation)點的計算,使用外推點的梯度完成當前點的更新.額外梯度法通過外推點的計算加大探索,使得更新迭代時更容易找到逼近最優解的正確方向,具備良好的收斂性.

3.2 馬爾科夫博弈策略梯度算法

基于額外梯度法提出適用于馬爾科夫博弈的策略梯度算法.該算法使用的策略梯度由式(26)獲得,但式中的期望在無模型設定下是無法計算的,只能獲得其采樣值.因此,在該設定下通過隨機策略梯度近似真實策 略梯度,即

基于此,提出馬爾科夫博弈策略梯度(Markov game policy gradient,MG-PG)算法,如算法1所示.

算法1 MG-PG算法

3.3 馬爾科夫博弈策略梯度算法的收斂性

當目標函數J(θ,?)滿足凸凹性時,是單調的.使用額外梯度算法求解式(15)的收斂性在長達幾十年的研究中已經被解決.但在馬爾科夫博弈問題中,目標函數J(θ,?) 是非凸非凹的,即不滿足單調性.基于此,首先提出以下假設.

假設1兩方零和馬爾科夫博弈(式(1))和參數化的聯合策略滿足以下假設:

Pethick等[18]給出了一類滿足弱Minty變分不等式(weak Minty variational inequality,MVI)的非凸非凹minimax問題的額外梯度型算法.本研究的問題設定中無約束項,因此文獻[18]的Assumption I(i)中的正則算子A=0.本研究假設1的條件1)是比一致連續更強的光滑性條件,其限制了策略梯度的局部變動幅度不能超過常量L,滿足文獻[18]的Assumption I(ii).本研究假設1的條件2) 可以確保 [θ*;?*] 是有限的,其中式(37)是MVI在兩方零和馬爾科夫博弈下的具體表述.根據推論1可知,隨機策略梯度

證明:基于假設1,根據文獻[18]定理3.5可知,當λ和αk滿足定理2中的取值范圍時,由MG-PG算法迭代生成的策略參數序列滿足

存在k∈{0,1,···,K},式(39)成立.

4 仿真實驗

4.1 實驗對象和評估指標

采用棋盤游戲Oshi-Zumo[23]來驗證算法MGPG的收斂性.該游戲中有2個玩家,是一種同時移動的零和博弈,其難點在于每次決策時的不確定性和信息不對稱性,即每個玩家無法知道對手的決策,需要更加復雜的策略來應對可能的動作.一輪博弈往往須經過多步博弈才能分出勝負,其游戲規則如下: 棋盤上有2N+1個格子一維排列,編號1,···,2N+1,在棋盤中心(第N+1個格子)有一面旗幟,旗幟的移動方向取決于每一輪博弈中玩家1、2的每步博弈結果(移動格數始終為1).如圖2所示,玩家1、2初始時各有X枚硬幣,每一步,玩家1、2同時出硬幣,記為M1和M2,然后比較M1和M2的大小,旗幟始終往出幣數少的一方移動,即若M1>M2,旗幟向右移動一個格子;若M1<M2,旗幟向左移動一個格子;若M1=M2,旗幟不移動.每一步博弈后2個玩家的總數減去出幣數M1和M2,并開始下一步博弈.比賽一直進行到雙方的硬幣都用完,或者旗幟被推出棋盤外.獲勝者根據旗幟的位置確定-距離旗幟所在的位置近的玩家輸掉比賽.如果旗幟的最后位置是棋盤中心,則比賽結果為平局.

圖2 Oshi-Zumo在不同設定下的狀態變化示意圖Fig.2 Schematic diagrams of state of Oshi-Zumo in different settings

Oshi-Zumo的狀態由3部分組成: 玩家1的剩余硬幣數、玩家2的剩余硬幣數和旗幟的位置.旗幟位置由格子編號表示,當旗幟從左移出第1個格子后,旗幟位置為0;當從右移出第2N+1個格子后,旗幟位置為2N+2.當游戲參數確定后,初始狀態也是確定的.玩家的動作就是出幣數,僅當博弈結束后才獲得獎勵.玩家1作為max玩家(期望收益為正),獲勝則獎勵為1;玩家2獲勝則獎勵為-1,平局則獎勵為0.

采用納什收斂指標[17,24-25]作為聯合策略的評價指標.給定聯合策略 (π,μ),納什收斂指標定義為

式中的最佳響應策略由單智能體強化學習算法訓練獲得.由定義1、2可知,當 NashConv(π,μ)=0 時,聯合策略 (π,μ) 達到納什均衡;當 NashConv(π,μ)<ε,?ε>0,聯合策略(π,μ) 為ε 近似納什均衡.在下面的實驗中,采用單智能體強化學習reinforce算法求解最佳響應策略,即固定一方玩家的策略,對另一方玩家的策略進行訓練,直到勝率達到90%或策略參數的更新次數達到5萬次.

4.2 表格式softmax參數化策略

在第k迭代輪次,將玩家在狀態s下的策略參數拼 接為一個參數向量其中分別表示玩家1、2在狀態s下合法動作個數.2個玩家的策略服從指數柔性最大化(softmax)分布[26],即

式中:[·]a表示第k輪次在狀態s下對所有不同的合法動作的選擇概率依次按照方括號內的規則進行計算.參數的初始值全為0,即初始策略服從均勻分布.

Oshi-Zumo規模設置如圖2(a)所示,2個玩家的幣數都為6,棋盤的格子數為5.實驗中MGPG算法的梯度步長 λ =0.3,αk=0.5,折扣因子 γ=1.0,最大更新次數kmax=106,每次更新采樣的局數(軌跡數)n=1,如表1所示.

表1 不同參數化策略設定下的算法超參數Tab.1 Algorithm hyperparameters under different parameterized policy settings

設置10組實驗,實驗的評估數據如圖3所示.圖中,k表示策略參數更新的次數,Na表示納什收斂指標的值,圓點曲線表示10組實驗納什收斂指標的均值及變化趨勢,豎線表示10組實驗納什收斂指標的離散程度,豎線的上下界分別由均值加減標準差得到.可以看出,實驗中隨著更新次數的增加,納什收斂指標雖然存在一定的方差,但其均值整體呈減緩下降趨勢,當更新次數達到約60萬次時,納什收斂指標的均值接近零,此時聯合策略達到近似納什均衡.

圖3 策略參數化下MG-PG算法的納什收斂指標Fig.3 Nash convergence of MG-PG algorithm with policy parameterized

4.3 神經網絡參數化策略

在如圖2(b)所示的游戲規模下,狀態動作空間大小急劇增長,使用表格式softmax參數化策略容易導致實驗過程中出現內存溢出的問題.因此,將神經網絡作為參數化策略,并通過MG-PG算法訓練優化神經網絡的參數,最終達到近似納什均衡.

在神經網絡作為參數化策略的實驗中,2個玩家的初始幣數都是50,棋盤中共有 2×7+1=15個格子,每個玩家最多有51個候選動作(包含出幣數為0的動作).因此首先定義維數為(51+51+15)=117的one-hot向量(狀態信息)作為神經網絡的輸入,神經網絡輸出維數為51.將神經網絡的輸出進行softmax運算,即得到輸入狀態下的策略.在該設定下由于其大規模的狀態空間和動作集,應該使用更多的隱藏層和更大的層大小.因此采用含2個隱藏層的非線性網絡,從輸入側到輸出側,隱藏層分別有256、128個節點,并使用ReLU激活函數、Adam優化器.關于算法的超參數設置,如表1所示.相較于使用表格式softmax參數化策略的算法表現,使用神經網絡作為參數化策略由于激活函數的存在,可能出現數值不穩定性的問題,如梯度爆炸或梯度消失.通過使用較小的梯度步長,可以減緩參數的更新,有助于緩解這些數值穩定性問題.同時由于游戲規模較大,其最大更新次數也較大.

策略網絡優化流程如下.采用在線學習方式,在第k輪次下,將環境狀態st作為輸入分別輸送到2個玩家的策略網絡 πθ和 μ?.依照策略網絡輸出的概 率分布,隨機得 到動作at和bt,在與環境交互后,得到新 的環境狀態st+1.重復以上步驟,得到一條博弈軌跡將該軌跡中的狀態 (st)0≤t≤T作為輸入依次輸送到策略網絡 πθ和 μ?.以策略網絡 πθ為 例,輸出概 率分布后計算梯度梯度項求和后以步長 λ 更新賦值給臨時策略網絡,臨時策略梯度的計算同理可得.最后通過臨時策略梯度以步長αk更新網絡參數[θ;?].至此通過算法迭代完成一次策略網絡更新.

在神經網絡設定下同樣設置了10組實驗,實驗的評估數據如圖4所示.與策略參數化設定相類似,盡管納什收斂指標存在一定的方差,但其均值整體呈減緩下降趨勢.當更新次數達到約500萬次時,納什收斂指標的均值接近零,此時聯合策略達到近似納什均衡.實驗結果表明,MGPG算法對于大規模博弈問題仍然適用.

圖4 神經網絡下MG-PG算法的納什收斂指標Fig.4 Nash convergence of MG-PG algorithm with neural network setting

4.4 對比實驗

在基于Oshi-Zumo的實驗中,MG-PG在2種參數設定下都展示出較好的性能.為了能夠更全面地評估和驗證本研究所提出的算法MG-PG,將現有的兩方零和馬爾科夫博弈策略梯度方法-獨立策略梯度方法[23]、嵌套梯度下降方法[24]、雙時間尺度算法[25]和演員-評論家虛擬博弈算法[8]在Oshi-Zumo中的表現與MG-PG進行對比分析.實驗的游戲規模與神經網絡下MG-PG的實驗規模相同,即2個玩家的初始幣數都為50,棋盤中共有15個格子.

獨立策略梯度方法本質上為單智能體強化學習,若2個玩家都使用常用的reinforce算法進行參數更新,則更新式為

式中:ηd為學習率.

實驗中玩家策略采用神經網絡的形式,參數ηd與MG-PG的 步長 αk保持一致,均為0.3.

嵌套梯度下降方法則將最大最小優化問題(式(15))轉化為如下最小優化問題:

因此,嵌套梯度下降方法在每一輪迭代時,先進行一次“嵌套內循環”,近似計算出 θ*后,隨后在??J(θ*,?) 的方向上更新參數?.實驗中玩家策略采用神經網絡的形式,超參數同為0.3.

將策略參數為 ? 的玩家視為領導者,策略參數為 θ 的玩家視為追隨者.根據雙時間尺度算法的思想,對領導者和追隨者之間進行時間尺度分離.在實驗中基于嵌套梯度下降方法,令領導者策略參數更新的超參數為0.1,小于追隨者的策略參數更新超參數,并使領導者策略參數更新間隔為10個迭代輪次,以達到領導者的學習速度比追隨者慢的效果.

演員-評論家虛擬博弈算法通過Actor-Critic框架對虛擬博弈過程進行隨機近似,使得算法能夠收斂到納什均衡.策略與價值函數的更新式如下.

式中:上標i表示玩家標號,ψk為算法的超參數,B?(Q)為最佳響應策略(選擇概率函數),

?(·) 為確定 性擾動,C?(Q)為B?(Q) 的平均 值,即

在實驗中步長選擇為 ηk=0.01,ψk=0.1.

每種算法同樣進行10組實驗,實驗結果如圖5所示.每種算法都迭代更新了800萬次.可以看出,獨立策略梯度方法、嵌套梯度下降方法和雙時間尺度算法3種方法在800萬的更新次數內不收斂,并且納什收斂指標也沒有明顯下降的趨勢.3種方法都使用了reinforce算法,因此在10組實驗中都呈現出較大的方差.演員-評論家虛擬博弈算法則在實驗中表現良好,納什收斂指標呈現出較為平滑的下降趨勢.但該算法在800萬次的更新中并未收斂到近似納什均衡,納什收斂指標的均值在更新后約為0.50,相較于MG-PG算法在500萬次達到約0.15還存在一定的性能差距.同時該算法在訓練過程中由于要計算最佳響應策略,對于時間的消耗也是巨大的.不過,得益于Actor-Critic框架值函數提供了對策略的評估,策略改進更加穩定,大大減小了不同實驗組策略優化過程中的方差.

圖5 4種對比算法的納什收斂指標圖Fig.5 Nash convergence index diagram of four comparison algorithms

綜上,單智能體強化學習的思想及其方法在Oshi-Zumo這種同時移動博弈游戲中的表現欠佳.演員-評論家虛擬博弈算法相較于本研究所提出的MG-PG算法收斂速度較慢,但其Actor-Critic框架帶來的小方差為本研究今后的工作提供了指導和啟示.

5 結語

提出兩方零和馬爾科夫博弈的策略定理及其相關證明,并基于隨機策略梯度和額外梯度法提出適用于馬爾科夫博弈的策略梯度算法.該算法在外推點梯度計算階段,利用已知的博弈軌跡計算外推點的梯度,改善探索方向,使得算法更容易朝正確方向更新;在策略參數更新階段,利用外推點的梯度更新策略參數.最終在假設1的條件下分析得到該算法的收斂性.

在Oshi-Zumo上的大量實驗表明,本研究算法經過一定的迭代輪次可以有效地收斂到近似納什均衡解.應用神經網絡擬合博弈策略,解決了大規模博弈問題.通過對比實驗驗證了MG-PG算法的優越性和有效性.

如何通過應用Actor-Critic Network減小MGPG算法方差以及如何解決環境部分可觀的馬爾科夫博弈問題,是未來可以繼續研究的方向.

猜你喜歡
馬爾科夫納什梯度
一個改進的WYL型三項共軛梯度法
基于疊加馬爾科夫鏈的邊坡位移預測研究
THE ROLE OF L1 IN L2 LEARNING IN CHINESE MIDDLE SCHOOLS
THE ROLE OF L1 IN L2 LEARNING IN CHINESE MIDDLE SCHOOLS
一種自適應Dai-Liao共軛梯度法
基于改進的灰色-馬爾科夫模型在風機沉降中的應用
一類扭積形式的梯度近Ricci孤立子
馬爾科夫鏈在教學評價中的應用
愛,納什博弈人生的真理
基于馬爾科夫法的土地格局變化趨勢研究
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合