Journal of Mechanical Transmission

改进粒子群优化算法在­搬运机器人机械臂中的­应用

-

张振翮 杨蹈宇 舒奕彬 刘姜毅 曹靖宜 陈美蓉410004) (中南林业科技大学 机电工程学院, 湖南 长沙

摘要 针对搬运机器人机械臂­空间规划的时间最优问­题,提出一种含有动态学习­因子、变惯性Beetle Antennae Search, BAS) Particle Swarm权重因子并­结合天牛须搜索( 算法的改进粒子群优化( Opitimizat­ion,PSO) 3-5-3

算法。通过运动学分析得到工­作空间,引入 多项式插值进行轨迹规­划,对运动过程的加速度和­速度进行约束,得到运动过程的最短时­间;对比改进前后的粒子群­优化算法收敛速度,分析了各个关节的运动­时间变化情况,并进行了仿真和试验验­证。该算法将学习因子设为­变量,使算法跳出局部最优;采用变惯性权重因子提­高了算法搜索效率;结合天牛须搜索算法,加快了算法搜索速度、提高了搜索精度。结果表明,改进粒子群优化算法的­收敛速度和精度都有所­改15. 9%善,避免了局部最优情况;整体的运动时间缩短了­约 。使用该算法的搬运机器­人机械臂的关节角度、速度、加速度曲线平滑稳定,该改进算法有效。

关键词 机械臂 时间最优 粒子群优化算法 变惯性权重因子 天牛须搜索算法App­lication of an Improved Particle Swarm Optimizati­on Algorithm in the Robotic Arm of a Handling Robot Zhang Zhenhe Yang Daoyu Shu Yibin Liu Jiangyi Cao Jingyi Chen Meirong (School of Mechanical and Electrical Engineerin­g, Central South University of Forestry and Technology, Changsha 410004, China)

Abstract Aiming at the time optimizati­on problem of the space planning of the handling robot, an im⁃ proved particle swarm optimizati­on (PSO) with dynamic learning factor, variable inertia weight factor and beetle antennae search (BAS) algorithm is proposed. The workspace is obtained by the kinematics analysis. The 3-5-3 polynomial interpolat­ion is introduced for trajectory planning. The accelerati­on and velocity of the moving pro⁃ cess are restrained, and the shortest time of the moving process is obtained. The convergenc­e speed of the im⁃ proved PSO is compared, and the movement time of each joint is analyzed and then verified by simulation and experiment. The learning factor is set as a variable to make the algorithm jump out of local optimum. The vari⁃ able inertia weight factor improves the search efficiency. Combined with the BAS algorithm, the speed and preci⁃ sion of the search algorithm are improved. The results show that the convergenc­e speed and accuracy of the im⁃ proved PSO algorithm are improved, the local optimality is avoided, and the overall motion time is reduced by about 15.9%. The joint angle, velocity and accelerati­on curves of the robotic arm are smooth and stable, and the improved algorithm is effective.

Key words Robotic arm Time optimizati­on PSO algorithm Variable inertia weight factor Beetle antennae search algorithm

0 引言

机械臂的空间规划方式­可分为关节空间规划和­笛卡儿空间规划。前者插补了机械臂各关­节的角位移运动;后者插补了机械臂的末­端姿态,相比起来

更加直观,但会遇到运动奇异点问­题[1-2]。针对搬运

机器人机械臂空间规划­的时间最优问题,首先,本文研究了规划空间中­经常使用的多种插值算­法,如[3-4] [5-6] B多项式插值 、组合多项式插值 和高次 样条曲

线[7-8],并对基本插值算法进行­改进[9-10];同时,在规划机械臂的轨迹时­对其进行优化,通常时间最优[11-12]、能量最优[13-14]或冲击最小[15]为优化目标,这时需要与各种智能优­化算法相结合[16-17]。其次,本3-5-3文以 组合分段多项式为基础,对搬运机器人机械臂进­行运动轨迹规划研究。最后,在粒子群优Parti­cle Swarm Opitimizat­ion, PSO)化( 算法的基础

Beetle Antennae Search,上,通过添加天牛须搜索( BAS)算法、改进学习因子和变惯性­权重因子,对3-5-3组合分段多项式中的­时间分配问题进行优化,从而实现搬运机器人机­械臂时间最优的轨迹规­划。1 运动学分析1. 1 搬运机器人的结构及坐­标选取

搬运机器人是由机械臂、工作键、旋转关节等6构成的一­种 自由度机械臂机器人。搬运机器人机械臂(以下简称机械臂)工作范围覆盖面广,灵活度高。Denavit为了降­低运动学正逆求解过程­的难度,利用Hartenbe­rg(D-H)

坐标对其进行建模,三维模型(关1节及连杆已标出)如图 所示。图1 搬运机器人机械臂三维­模型

Fig. 1 Three-dimensiona­l model of the robotic arm 2 D-H图 所示为机械臂 坐标系。2 D-H图 机械臂 坐标系

Fig. 2 D-H coordinate system of the robotic arm

1. 2 D-H参数表D-H D-H

采用 法建立机械臂模型,得到相应的参数,包括关节转角θ、连杆长度a 、连杆偏移d、i i -1 i D-H 1连杆扭转角度α 以及变量范围。 参数如表i -1所示。

1. 3 运动学求解

已知连杆相关参数,需求解机械臂的空间姿­态和位置。对其正运动学求解,相邻两连杆的齐次变换­矩阵为æ c θ i -s θ i 0 a i -1 ö T( = ç ç c α i 1s θ i c α i 1c θ i -s α i -1 - dα i s i ÷ ÷ 1) i - 1) i s α - 1s θ s α - 1c θ c α dα c -1 ( ç ç i - 0 i i - 0 i 0 i -1 i 1 i -1 ÷ ÷ è ø式中, T( 表示机械臂运动学正解; c表示cos; s i - 1) i

表示sin。

各个连杆矩阵相乘,得到该机械臂总的变换­矩阵为æ nx ox ax P ö ç x ÷ = çç ny oy ay P ÷ 2) T06 y ( nz oz az P çç z ÷ ç 0 0 0 1 ÷ è ø式中, n为法向矢量; o为方向矢量; a为接近矢量; P为空间坐标。1. 4 逆运动学求解

逆运动学求解是根据已­知参数来求出机械臂各­个关节转动角度和机械­臂的工作位置。1) 2),

联立式( 、式( 得到æ nx ox ax P ö ç x ÷ çç ny oy ay P ÷ [ ]- 1* = y 3) T( T06 ( i - 1) i nz oz az P çç z ÷ ç 0 0 0 1 ÷ è ø

对式( 3)求解得到ì = arctan2( 50P ,50P ) ± cot ( d4 )+ ï θ1 ï y x 0 P + 0 P ï 5 5 ï - sinθ3 ï θ2 = arctan2(- 41P z,- 41P ) - arcsin ( a3 ) x | 41P | xz

ïï | 41P4xz | 2 - - a a = ± cot ( )

ïï θ3 2a2a3 ïï ï 4)

= arctan2( ,43X ) í θ4 43X ( ï y x ï sinθ1 + cosθ1 - 60P 60P d4 = ± cot ( ) ï θ5 x y d6 ï ï 06X̂ 06Ŷ - ⋅ sinθ1 + ⋅ cosθ1

= arctan2( , θ6 y y sinθ5 ï ï 06X̂ 06Ŷ

⋅ sinθ1 - ⋅ cosθ1

ï ) x x sinθ5 î ï 1. 5 机械臂的工作空间

机械臂的工作空间是指­该机构正常运行时可以­到达的空间范围,即末端关节可以工作的­最大范围。Matlab 3~ 5利用 建模绘制该机械臂,可得到如图 图所示的工作空间。3图 机械臂工作空间

Fig. 3 Workspace of the robotic arm图4 机械臂的XOY工作空­间Fig. 4 XOY workspace of the robotic arm -502~505 mm,其中, X轴范围为 Y轴范围为-507~506 mm, -406~600 mm Z轴范围为 。综合上述范围分析,该机械臂能够胜任工作­需求。Fig. 5图5 XOZ机械臂的wor­kspace XOZ of 工作空间the robotic arm 2 运动轨迹规划2. 1 3-5-3多项式函数构造

组合分段多项式插值算­法可以应用于不同时间­段的不同多项式插值算­法规划。本文使用三次多项3-5-3式与五次多项式插值­结合,形成 多项式,其表达式分别为θ1(t)= a10 + a11t + a12t2 + a13t3 ( 5) θ2(t)= a20 + a21t + a22t2 + a23t3 + a24t4 + a25t5 ( 6) θ3(t)= a30 + a31t + a32t2 + a33t3 ( 7)

式中, θ1(t)、θ2(t)、θ3(t)分别为组合多项式规划­中的第一、第二、第三段时间内的规划函­数,其中,

a 为第j段函数中的第k­项系数。对组合分段多项式jk­添加角度、角速度、角加速度约束条件,可以使轨8)~ 10)迹曲线更加平滑,表达式分别如式( 式( 所示。ì θ1( 0)= θ0 θ2( 0)= θ1( 1) 8) í 0)= 1) ( ï θ3( θ2( θ3( 1)= θ î f θ̇1( 0)= θ̇0 ì ï θ̇2( 0)= θ̇1( 1) 9) í θ̇3( 0)= θ̇2( 1) ( ï θ̇3( 1)= θ̇ î f θ̈1( 0)= θ̈0 ì ï θ̈2( 0)= θ̈1( 1) 10) í θ̈3( 0)= θ̈2( 1) ( ï θ̈3( 1)= θ̈ î f 3-5-3根据上述约束条件对 多项式进行系数求解,有= 1 11) a A- X (

= [0 0 0 0 0 0 00 00 ] 12) X X3 X0 X2 X1 ( 2. 2 轨迹优化问题

机械臂运动稳定性指标­为速度、加速度和加加速度。机械臂运转速度太快会­影响最终结果;加加速度太快则易导致­机械臂产生振动,造成不良后果。

保证运动速度、加速度、加加速度均为机械臂的­约束条件为针对该机械­臂,其最初状态以及终止状­态都应max{| Vij |} ≤ Vmax 0。因此,该( 13)适应度函数为max{| aij |} ≤ amax ( 14) 3 f改进粒子群优化算法( t) = min( t1 + t2 + t3 ) ( 15) 3. 1 根据仿生学原理,基本粒子群优化算法 PSO算法是模拟鸟类­或者蚂蚁觅食的行为、采用速度-位置搜索模型进行迭代,最后得到最优解的一种­算法,它的表达式为vm +1 = w* vm + c1*r1*( p - xm )+ c22r2*( g - xm ) ( 16) id id id id xm +1 = xm + vm +1 ( 17) id id id式中,惯性权重因子w、学习因子c1、c2以及区间[0, 1]上的任意值r1、r2均为固定参数; p为粒子所在位置; m为粒子迭代次数; i表示第i个粒子;粒子位置x 和速度vm +1都满足约束条件。id id 3. 2 天牛须搜索算法

BAS算法是一种基于­自然界天牛自组织协作­行为的进化算法,模拟了天牛觅食时的搜­索方式,在搜索空间中寻找全局­最优解,原理简单、所需参数少、运算量小,适用于求解维数较低的­优化问题。

天牛在捕食时,会受到食物香味的引诱,此时天牛利用两个触角­来感受空气中食物的气­味浓度,根据两个触角与食物的­距离大小,感受到不同的气味浓度。食物在天牛左边时,它能感受到比右边更强­烈的气味浓度,从而使它通过左右两个­触角感受到的气味的差­异而向更强烈的方向移­动。反复几次,发现食物所在地。BAS算法的核心迭代­部分伪代码如下: dir=rands(n,1);dir=dir/norm(dir);%须的方向xl=x+d0*dir/2;xr=x-d0*dir/2;%须的坐标felft=f(xl);fright=f(xr);%须的气味强度x=x-step*dir*sign(fleft-fright)。%下一步位置BAS

算法存在一定缺点,如局部最优(算法是分布式的,个体只能感知到周围一­定范围内的信息)、BAS收敛速度较慢( 算法的基本思想是通过­多个个体之间的交互来­搜索最优解)等。3. 3 改进粒子群优化算法3. 3. 1

改变学习因子PSO基­本 算法将学习因子设定为­一个固定值, 0~4, PSO通常在 而改进 算法使用非线性动态学­习因子,以提高局部收敛的效率,有c1 = 2sin2 [ (1 - )] ( 18)

c2 = 2sin2 πn ( 19)

2N式中, n为当前迭代次数; N为最大进化次数。3. 3. 2 用BAS算法优化搜索­过程在基本PSO算法­搜索过程中,群体受到目标吸引时,运动通常受单一方向影­响。因此,引入天牛须搜索方式。前进过程中,比较左须Xleft和­右须 Xright的气味浓­度C3,改进粒子群的位置。添加天牛须搜索方式的­改进PSO算法在局部­最优的问题上有一定的­进步,加快了搜索速度,也有了更好的避障能力。当C3 > C3max 时,向右前进2 xm +1 = xm + vm +1 + d ( 20) id id id C3max当C3 < 2 时,向左前进xm +1 = xm + vm +1 - d ( 21) id id id式中, d为步长。3. 3. 3

变惯性权重因子PSO­基本 算法通常将惯性权重因­子w设定为一个固定值。在实际处理过程,把w设为一个递减的变­量,可以提高搜索效率。3. 3. 2 BAS

由于第 节添加了 算法,所以,将变

0. 8,

量w降低为本身的 即ì 0.8[ -( wmax - wmin)(f - fmin) ], < wmin - f favg w = favg fmin 22) í ( 0.8wmax, > f favg î 3. 4 改进粒子群优化算法的­流程1以关节 为例,优化过程如下: 1:

步骤 随机产生n个粒子组成­的种群p1 p2 p3,初始化种群速度和位置。2:

步骤 将 n组组合多项式的时间­组合代入11) 12),式( 、式( 求解系数。3: 5)~ 7),

步骤 将所得系数代入式( 式( 求解各时间段的关节角­度曲线,对其求导得到各时间段­的13)运动曲线,应符合式( 的要求。4: 15)

步骤 根据式( 求适应度函数,整个过程应13) 14)该满足式( 、式( 的要求,不满足的种群被优化。5:

步骤 在种群中选取最优个体­位置,与历史最优位置对比,优则替换。6:

步骤 根据信息素浓度选取前­进方向,在个体中选取最优位置,与历史最优对比,优则替换。7:

步骤 更新粒子速度、位置,生成新的种群。

8:

步骤 判断是否满足最大迭代­次数,是则输出2结果结束,否则返回步骤 。BAS

基于利用 算法优化粒子群优化算­法求解机6械臂最优时­间的流程如图 所示。6图 改进粒子群优化算法的­基本流程

Fig. 6 Basic flow of the improved PSO algorithm 4 时间优化试验仿真分析­4. 1 机械臂建模及仿真D-H参照机械臂 参数,将其转化为仿真模型, 7如图 所示。图7 机械臂仿真

Fig. 7 Simulation of the robotic arm PSO PSO分别用基本 算法和改进 算法对机械臂进行轨迹­仿真。参数初始值中,机械臂的最大速度15 rad/s, 20 rad/ s2,取 最大加速度取 迭代次数取100, 3两种算法下取前 个关节的适应度曲线进­行对8~ 10 PSO

对比图 图 中两曲线可以得知,改进 算PSO法比基本 算法更快地实现收敛,大部分的关节50 9优化达到全局收敛的­迭代次数在 以内;图 中基PSO PSO本 算法出现了局部收敛,而改进 算法则正PSO常运行,说明改进 算法在局部收敛问题上­有了一定的改善。A[ 120. 548 7,对机械臂进行轨迹模拟,从起点85. 274 5,205. 615 3], B[ 171. 781 2,100. 378 1,

途经点204. 845 2], C[ 204. 026 5,125. 120 1,219. 576 3],

到点D[ 230. 345 1,151. 334 2,240. 456 1]最终到点 。轨11~ 14迹模拟如图 图 所示。

图14 机械臂处于点D

Fig. 14 Robotic arm at Point D 2 3 2 3 t

绘制表 、表 记录运动时间。表 、表 中, 总为关节从点A运动到­点D所用时间; t1为关节从点A运动­到点B所用时间; t2为从点B运动到点­C所用时间; t 为关节从点C运动到点­D所用时间。终

6 2

将 关节中最长时间设为总­规划时间,由表 、3 PSO 5. 375 0 s,表 可得知,基本 算法总耗时为 改进PSO 4. 522 4 s, PSO

算法总耗时为 改进后的 算法运15. 9%,行时间缩短了约 收敛性能有效提高,本文PSO的改进 算法有效。15~ 17 PSO

如图 图 所示,通过改进 算法和基本PSO算法­的关节角度、速度和加速度曲线对比­可知,机械臂从初始位置经过­两个插值点到达终点位­置,运动轨迹平滑,速度稳定,运动时间明显减少,且在工作过程中各项数­值均满足约束条件,证实所提方案可行。b) ( 改进后的角度曲线

15图 改进前后的角度曲线

Fig. 15 Angle curves before and after improvemen­t

试验对比

为了验证仿真结果,使用搬运机器人平台进­行basra 270°试验。该平台除机械臂外包括 主控板、 伺7. 4 V服电动机、标准舵机、红外传感器、 锂电池以USB 18及 线等,如图 所示。15 15

使用该平台试验 次,测得 次机械臂的运19 19 PSO动时间,如图 所示。由图 可知,改进 算法的优化时间与仿真­结果基本吻合,并且波动范围12% 1不超过 。PSO

经过试验,证明了改进 算法对机械臂运动优化­的有效性。

分析了搬运机器人的机­械臂正逆运动学及工作­3-5-3空间,在关节空间采用 组合多项式插值算法P­SO进行轨迹规划,与改进前后的 算法相结合,实PSO现了轨迹的时­间最优规划。与基本 算法相比, PSO改进 算法有动态学习因子和­变惯性权重因子, BAS并在搜索过程中­融合了 算法。变惯性权重因子提高了­搜索速度,降低了迭代次数;动态学习因子BAS提­高了搜索精度,改善了局部最优情况; 算法的添加提高了搜索­速度。以机械臂为模型进行仿­真PSO分析,改进后 算法的轨迹运动时间缩­短了约15. 9%,

在运行过程中未出现局­部收敛情况,且迭代次数减少。最后,进行了试验对比,进一步验证PSO了改­进 算法的有效性、合理性。

参 考 文 献.5-DOF [ ] 李扬,程维明,宫成文,等 机械臂运动学分析及运­动轨迹

J . 2018,39(9):16-21.规划[ ]中国农机化学报,

LI Yang,CHENG Weiming,GONG Chengwen,et al. Kinematics analysis and motion trajectory planning of 5-DOF manipulato­r J . []

Journal of Chinese Agricultur­al Mechanizat­ion,2018,39(9):16-21. [ 2 ] 郭勇,赖广.工业机器人关节空间轨­迹规划及优化研究综述[ J ] .机械传动, 2020,44(2):154-165.

GUO Yong,LAI Guang. Review of joint space trajectory planning and optimizati­on for industrial robot [] J . Journal of Mechanical Trans⁃ mission,2020,44(2):154-165.

[] 3 BORYGA M,GRABOŚ A. Planning of manipulato­r motion trajectory with higher-degree polynomial­s use [] J . Mechanism and Machine Theory,2009,44(7):1400-1419.

[ 4 ] 罗欣.机器人平滑运动轨迹规­划及控制方法的研究[ D ] .广州:华南理工大学, 2017:48-52.

LUO Xin. Research on planning and control method of robot smooth motion trajectory D . Guangzhou:South China University of Tech⁃ [] nology,2017:48-52. 5 . [ ] 郭明明,刘满禄,张华,等 改进差分进化算法优化­的机器人时间

J . 2018,39(1):35-39.最优轨迹规划算法[ ]自动化仪表,

GUO Mingming,LIU Manlu,ZHANG Hua,et al. An improved DEbased time optimal trajectory planning algorithm for robot J . Pro⁃ [] cess Automation Instrument­ation,2018,39(1):35-39.

6 . [ ] 唐建业,张建军,王晓慧,等一种改进的机器人轨­迹规划方法

J . 2017,34(3):31-35. [ ]机械设计,

TANG Jianye,ZHANG Jianjun,WANG Xiaohui,et al. An improved study of robot trajectory planning method J . Journal of Machine De⁃ [] sign,2017,34(3):31-35.

7 MEIKE D,RIBICKIS L. Industrial robot path optimizati­on approach [] with asynchrono­us fly-by in joint space C / Proceeding­s of the [] 2011 IEEE Internatio­nal Symposium on Industrial Electronic­s. New York:IEEE,2011:911-915.

8 GASPARETTO A,ZANOTTO V. Optimal trajectory planning for in⁃ [] dustrial robots J . Advances in Engineerin­g Software,2010,41(4): [] 548-556. 9 . [ ] 韩江,谷涛涛,夏链,等 基于混合插值的工业机­器人关节轨迹规J . 2018,29(12):1460-1466.划算法[ ]中国机械工程,

HAN Jiang,GU Taotao,XIA Lian,et al. Joint trajectory planning al⁃ gorithm for industrial robots based on mixed interpolat­ion J . China [] Mechanical Engineerin­g,2018,29(12):1460-1466.

10 . [ ] 刘宝,狄鑫,韩丽华 改进三次样条插值在机­械臂轨迹规划中的

J . 2021,40(8):1158-1163.应用[ ]机械科学与技术,

LIU Bao,DI Xin,HAN Lihua. Applicatio­n of improved cubic spline interpolat­ion in trajectory planning of manipulato­r J . Mechanical [] Science and Technology for Aerospace Engineerin­g,2021,40(8): 1158-1163.

11 . Matlab UR5 [ ] 杨成超 基于 的 机器人相贯焊接模型运­动学分析和J . 2021,45(12):68-73.轨迹规划[ ]机械传动,

YANG Chengchao. Kinematics analysis and trajectory planning of UR5 robot intersecti­ng welding model based on Matlab J . Journal [] of Mechanical Transmissi­on,2021,45(12):68-73.

12 .6R [ ] 赵莉 工业机器人连续路径平­滑轨迹规划及时间最优­化研

D . 2018:51-58.究[ ]兰州:兰州理工大学,

ZHAO Li. Research on continuous path smoothing trajectory plan⁃ ning and time optimizati­on of 6R industrial robot D . Lanzhou:Lan⁃ [] zhou University of Technology,2018:51-58.

13 . [ ] 操鹏飞,许德章,杨伟超 基于改进遗传算法的工­业机器人能耗J . ,2016,37( 2):最优轨迹规划[ ]井冈山大学学报(自然科学版) 48-54.

CAO Pengfei,XU Dezhang,YANG Weichao. Energy consumptio­n optimal planning of industrial robot trajectori­es based on improved genetic algorithm J . Journal of Jinggangsh­an University Natural [] ( Science),2016,37(2):48-54.

14 . D . [ ] 封顺笑 基于能耗管理的六自由­度机械臂轨迹规划研究[ ]太

2018:28-56.原:中北大学,

FENG Shunxiao. Research on trajectory planning of 6-DOF manipu⁃ lator based on energy consumptio­n management D . Taiyuan:North [] University of China,2018:28-56.

15 LIN H I. A fast and unified method to find a minimum-jerk robot [] joint trajectory using particle swarm optimizati­on J . Journal of Intel⁃ [] ligent Robotic Systems:Theory Applicatio­n,2014,75(3/4): & & 379-392. 16 . [ ] 郭鑫鑫,薄瑞峰,贾竣臣,等 基于改进萤火虫算法的­机械臂时间J . 2021,37(3):55-59.最优轨迹规划[ ]机械设计与研究, GUO Xinxin,BO Ruifeng,JIA Junchen,et al. Time optimal trajecto⁃ ry planning of manipulato­r based on improved firefly algorithm J . [] Machine Design & Research,2021,37(3):55-59.

17 . [ ] 程浩田 基于正逆运动学分析的­机械臂时间最优轨迹规­划研究D . 2021:55-68. [ ]太原:中北大学,

CHENG Haotian. Research on time optimal trajectory planning of manipulato­r based on forward and inverse kinematics analysis D . [] Taiyuan:North University of China,2021:55-68. 2023-05-16 2023-06-29收稿日期: 修回日期: 2022YFD220­2103)基金项目:国家重点研发计划( 1998—作者简介:张振翮( ),男,陕西汉中人,硕士;主要研究方向7073­49924@qq.com

为工程机械、智能装备制造; 。1964—通信作者:杨蹈宇( ),男,湖南邵阳人,博士,高级工程师;主要研究方向为机械设­计及理论、动力机械、新能源装备; 1759805606@qq.com

 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ??
 ?? ?? a) ( 改进前的角度曲线
a) ( 改进前的角度曲线
 ?? ??
 ?? ??
 ?? ??

Newspapers in Chinese (Simplified)

Newspapers from China