CN111583298B - 一种基于光流法的短时云图追踪方法 - Google Patents

一种基于光流法的短时云图追踪方法 Download PDF

Info

Publication number
CN111583298B
CN111583298B CN202010332583.3A CN202010332583A CN111583298B CN 111583298 B CN111583298 B CN 111583298B CN 202010332583 A CN202010332583 A CN 202010332583A CN 111583298 B CN111583298 B CN 111583298B
Authority
CN
China
Prior art keywords
cloud
image
time
optical flow
sky
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010332583.3A
Other languages
English (en)
Other versions
CN111583298A (zh
Inventor
齐冬莲
闫云凤
张淞杰
董哲康
杨心怡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Yuanjian Information Technology Co ltd
Zhejiang University ZJU
Original Assignee
Hangzhou Yuanjian Information Technology Co ltd
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hangzhou Yuanjian Information Technology Co ltd, Zhejiang University ZJU filed Critical Hangzhou Yuanjian Information Technology Co ltd
Priority to CN202010332583.3A priority Critical patent/CN111583298B/zh
Publication of CN111583298A publication Critical patent/CN111583298A/zh
Application granted granted Critical
Publication of CN111583298B publication Critical patent/CN111583298B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/207Analysis of motion for motion estimation over a hierarchy of resolutions
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/13Edge detection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10004Still image; Photographic image
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10024Color image
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30181Earth observation
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04SSYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00Systems supporting electrical power generation, transmission or distribution
    • Y04S10/50Systems or methods supporting the power network operation or management, involving a certain degree of interaction with the load-side end user applications

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于光流法的短时云图追踪方法。首先通过全天空成像仪采集实时天空图片,并将图像中的建筑等无关因素除去;其次利用图像的RBG色彩信息将云和天空的像素区分开;最后通过对两帧的同一云像素进行对比得出云的边缘位置和速度矢量,预测下一时刻云的边缘位置并与太阳位置做对比实现预测。本发明能够实现对云图位置的识别预测,具有较高的准确率,并且具有速度快,通用性高等优点,具有良好的鲁棒性,能够应用于辐照强度预测或太阳能光伏功率输出预测系统。

Description

一种基于光流法的短时云图追踪方法
技术领域
本发明涉及了一种短时云图追踪方法,尤其是涉及一种基于光流法的短时云图追踪方法。
背景技术
云层的漂移可能会遮挡太阳,从而导致太阳能发电系统出现明显的功率输出波动。预测云层的位置可以提供足够的时间来平滑功率输出。目前常用的光伏功率预测方法都是基于历史数据进行数据分析实现的功率预测,时间尺度基本都是以天为单位,因此直接忽略云层移动对功率造成的影响。而传统的天气预报等气象数据只能预测一个时间段内的温度、降雨概率、风速等数据,时间尺度都是以小时为单位,要实现预测因云层的遮挡带来的短期功率波动,需要实时预测短时内的云层位置变化,而短时云层的位置预测只能通过图像实现。因此,如果能够实现基于图像的短时云图追踪,对于实现电力系统的智能化具有重要意义。而现有的云图追踪方法通常是将云视为形状不变的物体,通过设定质心进行跟踪预测,不符合现实中的云层运动规律。此外,现有的云图追踪方法通常都是直接选取某一种云像素区分算法对所有种类云图进行区分,并没有针对不同云种类进行算法的选择和优化,因此区分效果常常不太准确,而本发明可以实现更好的区分效果,从而提高追踪的准确性。
发明内容
为了解决背景技术中的问题,本发明提出了一种基于光流法的短时云图追踪方法,能快速预测云图位置并具有良好准确性。
本发明采用以下技术方案:
一种基于光流法的短时云图追踪方法,包括以下步骤:
1)通过全天空成像仪设备采集实时天空图片;
2)对采集到的所有天空图片进行预处理,消除无关的图像信息(如建筑物,屋顶,避雷针等),只保留图片中天空的区域;
3)将步骤2)获得的天空区域通过颜色信息进行云像素与晴空像素的区分。颜色信息主要分为三个通道,即R(红色)通道,G(绿色)通道和B(蓝色)通道,云像素和晴空像素在这三个通道的特性是不同的。通过这三个通道,可以定义一些指标来区分云和晴空像素;所使用的这些指标在面对不同种类的云时拥有不同的识别效果,因此要先获得每个指标在面对不同种类的云时对应的最优阈值,并选取最优阈值作为区分云和晴空像素的阈值;
4)对步骤3)得到的RGB信息区分云像素和晴空像素,以1代表云像素,0代表晴空像素,将图像变成二值图像,从而得到每张图片的二进制图像,通过对相邻两张图像的边缘云像素进行检测,可以得到云的边缘位置和速度矢量;
5)通过步骤4)中得到的云的边缘位置和速度矢量,可以预测下一时刻的云的边缘位置,之后将下一时刻云的边缘位置和太阳做对比,判断太阳能否被云遮挡。
上述技术方案中,进一步地,所述的全天空成像仪的光学系统前部安装有全画幅鱼眼镜头。
进一步地,所述步骤2)中,对采集到的所有天空图片进行预处理,具体过程如下:
2.1)首先使用彩色图像边缘提取算法生成彩色边缘图像,并通过加权平均方法获得彩色图像相应的灰度图像;
2.2)采用填充算法对步骤2.1)得到的图像进行处理,生成二元掩模图像,并通过点积运算去除无关的图像信息。其中,由于全天空相机的位置是固定的并且周围建筑物的边界形状相对不变,因此掩模在同一场景中是通用的,可以通过只对一张图片进行处理得到二元掩模图像,并用于该全天空成像仪捕获的所有图像,
2.3)由于边缘失真,因此将天顶角设定为70°(为了在基于全天空图像的应用中进行最佳评估,建议使用70°的天顶角)。在本算法中仅包括FOV(视场角)≤70°的像素区域,将大于70°天顶角的图像范围也通过掩模隐去。
进一步地,所述步骤3)中进行云像素与晴空像素的区分采用的方法为:目前最常用的快速区别云像素和晴空像素的方法有:RBR,BRD,BRBG,BRBGB以及MCE方法。上述5种经典算法均为基于阈值的方法,通过设定合适的阈值,这些方法可以作为区分云像素和晴空像素的指标。每种算法均存在最优阈值区间,该最优阈值区间受到诸多因素的影响(例如拍摄设备、云的类别等等)。
更进一步地,为了获得最优阈值区间,需要参考二进制图像。通过Photoshop,可以人为合成晴空的参考二进制图像。同时,每个图像的RGB层也可用于获得云图的参考二进制图像。
更进一步地,所述的获得最优阈值区间,具体为:针对不同的云种类,对实验数据遍历阈值,找到属于每种云种类的最优阈值区间,并通过比较不同算法分别对不同种类云的云/晴空像素识别准确率,找到5种算法中最适合每一种云的图像识别方法。
进一步地,所述步骤4)具体为:
对于在时间t=n和t=n+1时拍摄的两个连续图像,基于两个连续图像之间的云像素变化信息,可以得到t=n+1时的云图的速度矢量,该速度矢量可以预测云遮挡太阳的时间。为了实现这一点,选择光流法算法来追踪云像素的变化。
更进一步地,所述的光流法为Lucas-Kanade光流法。
进一步地,所述步骤5)中,预测下一时刻的云的边缘位置具体为:通过当前时刻检测到的云的边缘以及速度矢量,结合太阳所在位置,可以用卡尔曼滤波或概率密度等方法估计每个边缘像素点经过多长时间将会覆盖太阳,并以最短时间作为估计的云层遮挡太阳的时间。
本发明中,所述的光流法具体为:
为了计算像素点在时间(t+1)时的位置,一般通过视觉特征进行像素点的跟踪。目前常用的视觉跟踪算法是光流法。光流法一般通过一系列图像计算每个图像中每个像素的移动速度和移动方向。例如,当第t帧中的某一个点的位置是(x,y)并且第(t+1)帧中的位置是(x+u,y+v)时,可以确定位移是(u,v)。同时,也可以通过像素点的位移速度来估计它与观察者之间的距离,因为在视线中远距离的物体比近距离的物体移动得慢很多。
因为Lucas-Kanade算法在相邻图像之间的同一像素位移速度较慢时效果最好,因此本发明采用光流法中的Lucas-Kanade(LK)算法。LK光流方法是基于梯度的局部参数化光流估计方法。它假设光流向量在空间极小的邻域中是恒定的,并通过加权最小二乘法估计光流。LK光流法计算每个像素在时间t到t+Δt的两个帧的位移。由于LK光流法基于图像信号的泰勒级数,因此该方法属于微分方法的光流估计。
Lucas-Kanade光流法基于以下三个假设:
假设1:亮度恒定。为了追踪图像中的对象,对象的灰度值应该在短时间内保持不变。假设I(x,y,t)和I(x+dx,y+dy,t+dt)是两张连续图像中的同一像素点,图像约束方程定义如下:
I(x,y,t)=I(x+dx,y+dy,t+dt)
其中x和y是图像中像素的位置,t表示图像的时间序列。
假设2:时间连续或位移很小。当位移足够小时,可以计算灰度值相对于位置和时间的偏导数,该假设适用于计算瞬时速度而不是平均速度。如果没有这个假设,就无法跟踪像素的位置,也意味着Lucas-Kanade光流算法无法实现。同时通过使用泰勒级数作为图像约束方程,可以得到:
Figure BDA0002465495240000041
其中H.O.T代表二阶或更高阶。由于位移很小,H.O.T的值可以省略,因此约束方程也可以如下表示:
Figure BDA0002465495240000042
缩写该方程,可以得到如下表达式:
Figure BDA0002465495240000043
Figure BDA0002465495240000044
Vx=dx/dt,并且Vy=dy/dt。上述表达式可以用如下矩阵表示:
Figure BDA0002465495240000045
假设3:邻域的光流是一致的。上式矩阵中有两个未知数Vx和Vy而只有一个方程,因此无法求解Vx和Vy的值。为了解决这个问题,需要引入第三个假设:假设光流(Vx,Vy)在大小为m*m(m>1)的小窗口中为常数。此时可以从像素点1,2,...,n获得以下集合(其中n=m2):
Ix1Vx+Iy1Vy=-It1
Ix2Vx+Iy2Vy=-It2
Figure BDA0002465495240000046
IxnVx+IynVy=-Itn
对于上述集合,存在两个未知数,而存在两个以上的方程,这意味着该方程组是超定的,即方程中存在冗余,该方程有解。
本发明的有益的效果是:
本发明基于已有的云图识别算法,设计了针对不同的云种类选取不同的像素区分算法,并设定对应的最优阈值的云图识别方法。然而现有的云图追踪算法基本都将云层视为形状不变的物体进行追踪预测,当云层发生形变后预测准确率便大幅降低。但是本发明方法基本不受云图形状的影响,且在云图形状变化的情况下也能够准确地预测云的边缘形状,具有较强的适应能力。
本发明方法考虑了云与太阳的遮蔽关系对光伏功率输出的影响,能够有效地捕获太阳被云层遮蔽时光伏功率的下降突变以及云层离开太阳时光伏功率的上升突变,对电力系统安全、稳定的运行具有重大的意义。
附图说明
图1为实施例数据采集设备图片。
图2为实施例本方法步骤图。
图3为实施例本方法在不同指标下的云图检测与对比。
图4为实施例本方法前后两帧云图对比图。
图5为实施例本方法实现的云图运动追踪。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细说明。
本发明的实施例如下:
本发明中使用的图片由香港理工大学实验室屋顶上的商业全天相机(型号:SRF-02)获得(纬度:22.30,经度:114.18),如图1所示。从图1中可以看出,整个系统被封装在防水防风箱中,玻璃穹顶指向天顶,使其能够在白天(早上6点至晚上18点)在任何天气条件下工作。值得注意的是,安装在摄像机入射光学系统前部的特殊全画幅鱼眼镜头负责扩展视场角(FOV)。摄像机被编程为每5秒捕获一对图像。一个是正常曝光图像,另一个是具有相同场景的相关曝光不足图像。所有图像都以24位(每个RGB通道为8位)JPEG格式存储,最大分辨率为1600x1200。我们设定了六个数据集来验证光流法,由于摄像机被编程为每5秒捕获一次图像,因此实际时间的精度以5秒为单位。
将上述数据集中不同类型云使用不同的指标进行云像素识别处理,识别而得到的结果如图3所示。可以观察到RB Ratio指标在识别层云时具有比其他指标更好的效果,而BRDiff在识别积云时具有更好的效果,BRBGB在识别卷云方面表现更好,MCE在识别晴空时表现更好。同时,上述几种云中效果差异最明显的云是卷云。
采用上述数据集测试图片,得到的结果如表1和表2所示:
表1不同指标的最优阈值
Figure BDA0002465495240000051
Figure BDA0002465495240000061
表2光流法预测结果
数据集 预测时间 实际时间 误差率(%)
1 30 35 14.3
2 30 40 25.0
3 30 30 0.0
4 30 25 20.0
5 30 35 14.3
6 30 40 25.0
表1显示了在本测试数据下不同方法的最佳阈值。很明显,每个标准的最佳阈值变化很大。该表还显示云的类型影响阈值。由此可得出以下结论:当到达不同的云时,最佳阈值会发生变化。因此,建议在云像素识别之前执行云分类任务。
表2显示了预测时间与云遮挡的实际时间以及预测误差的结果。对于给定的六个数据集(即,数据集1至数据集6),预测误差均小于25%,六个数据集的平均误差为16.4%,误差在可接受范围内。
并且,本发明实施例对多种不同云图情况进行实施,实施例结果如图3所示,可以看到针对不同的云的种类,可以实施不同的评价指标来实现云像素和晴空像素的检测。
图5为常规情况下本发明方法的云图运动追踪结果示意。
由此可见,本发明能够实现云图的自动识别,具有较高的准确率,并且具有稳定性好,抗干扰能力强,通用性高等优点,同时也能对云图的运动轨迹轨迹进行实时追踪,能够应用于光伏功率预测系统或者是智能微电网系统。
上述具体实施方式用来解释说明本发明,而不是对本发明进行限制,在本发明的精神和权利要求的保护范围内,对本发明作出的任何修改和改变,都落入本发明的保护范围。

Claims (8)

1.一种基于光流法的短时云图追踪方法,其特征在于,包括以下步骤:
1)通过全天空成像仪采集实时天空图片;
2)对采集到的所有天空图片进行预处理,消除无关的图像信息,只保留图片中天空的区域;
3)将步骤2)获得的天空区域通过图像的RBG颜色信息进行云像素与晴空像素的区分;
4)通过步骤3)得到每张图片的二值图像,通过对相邻两张图像的边缘云像素进行追踪,得到云的边缘位置和速度矢量;
5)根据步骤4)中得到的云的边缘位置和速度矢量,预测下一时刻的云的边缘位置并与太阳位置做对比,判断太阳能否被云遮挡;
所述步骤2)中对采集到的所有天空图片进行预处理,具体采用以下方法:
2.1)首先使用彩色图像边缘提取算法生成彩色边缘图像,并通过加权平均方法获得彩色图像边缘相应的灰度图像;
2.2)采用填充算法处理步骤2.1)得到的图像,生成二元掩模图像,并通过点积运算去除无关的图像信息;对一张图片进行处理得到的二元掩模图像可用于该全天空成像仪捕获的所有图像;
2.3)由于边缘失真,因此将天顶角设定为70°,大于70°天顶角的图像范围通过掩模隐去;
为了获得最优阈值区间,需要参考二进制图像,所述的参考二进制图像可通过Photoshop人为合成或根据每个图像的RGB层获得。
2.根据权利要求1所述的一种基于光流法的短时云图追踪方法,其特征在于:所述的全天空成像仪的光学系统前部安装有全画幅鱼眼镜头。
3.根据权利要求1所述的一种基于光流法的短时云图追踪方法,其特征在于:所述步骤3)中,所述的区分云像素和晴空像素的方法为:RBR,BRD,BRBG,BRBGB或MCE方法,每种算法均存在最优阈值区间。
4.根据权利要求3所述的一种基于光流法的短时云图追踪方法,其特征在于:所述的获得最优阈值区间,具体为:针对不同的云种类,对实验数据遍历阈值,找到属于每种云种类的最优阈值区间,并通过比较不同算法分别对不同种类云的云/晴空像素识别准确率,找到最适合每一种云的图像识别方法。
5.根据权利要求1所述的一种基于光流法的短时云图追踪方法,其特征在于:所述步骤4)具体为:
对于在时间t=n和t=n+1时拍摄的两个连续图像,基于两个连续图像之间的云像素位置变化信息,得到t=n+1时的云图的速度矢量,该速度矢量可用于预测云遮挡太阳的时间;
所述的云像素的位置变化采用光流法进行追踪。
6.根据权利要求5所述的一种基于光流法的短时云图追踪方法,其特征在于:所述的光流法为Lucas-Kanade光流法。
7.根据权利要求1所述的一种基于光流法的短时云图追踪方法,其特征在于:所述步骤5)中,预测下一时刻的云的边缘位置,具体为:通过当前时刻检测到的云的边缘以及速度矢量,结合太阳所在位置,用卡尔曼滤波或概率密度方法估计每个边缘像素点经过多长时间将会覆盖太阳,并以最短时间作为估计的云层遮挡太阳的时间。
8.根据权利要求1-7任一项所述的一种基于光流法的短时云图追踪方法,其特征在于:可应用于辐照强度预测或太阳能光伏功率输出预测系统。
CN202010332583.3A 2020-04-24 2020-04-24 一种基于光流法的短时云图追踪方法 Active CN111583298B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010332583.3A CN111583298B (zh) 2020-04-24 2020-04-24 一种基于光流法的短时云图追踪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010332583.3A CN111583298B (zh) 2020-04-24 2020-04-24 一种基于光流法的短时云图追踪方法

Publications (2)

Publication Number Publication Date
CN111583298A CN111583298A (zh) 2020-08-25
CN111583298B true CN111583298B (zh) 2022-09-06

Family

ID=72112534

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010332583.3A Active CN111583298B (zh) 2020-04-24 2020-04-24 一种基于光流法的短时云图追踪方法

Country Status (1)

Country Link
CN (1) CN111583298B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115170619B (zh) * 2022-06-10 2023-08-15 山东电力建设第三工程有限公司 一种基于密集光流法的云遮挡预测方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10444406B2 (en) * 2014-04-17 2019-10-15 Siemens Aktiengesellschaft Short term cloud coverage prediction using ground-based all sky imaging
CN110514298B (zh) * 2019-08-30 2022-04-05 河海大学常州校区 一种基于地基云图的太阳辐照强度计算方法

Also Published As

Publication number Publication date
CN111583298A (zh) 2020-08-25

Similar Documents

Publication Publication Date Title
CA3100569C (en) Ship identity recognition method base on fusion of ais data and video data
CN109416413B (zh) 太阳能预报
CN108596129B (zh) 一种基于智能视频分析技术的车辆越线检测方法
CN102252623B (zh) 基于视频差异分析的输电线路导/地线覆冰厚度测量方法
CN101586956B (zh) 基于单目摄像机的河流水位监测方法
WO2014190651A1 (zh) 一种基于地基云图的光伏功率预测方法
CN112288736B (zh) 一种基于图像的能见度估测方法
Najiya et al. UAV video processing for traffic surveillence with enhanced vehicle detection
CN113159466B (zh) 一种短时光伏发电功率预测系统及方法
Alonso et al. Short and medium-term cloudiness forecasting using remote sensing techniques and sky camera imagery
CN108804992B (zh) 一种基于深度学习的人群统计方法
CN105611244A (zh) 一种基于球机监控视频的机场外来异物检测方法
CN109711256B (zh) 一种低空复杂背景无人机目标检测方法
CN111879292B (zh) 一种海岸线动态监控方法、监控设备及存储介质
CN102663385B (zh) 一种星上点目标检测方法
CN115170619B (zh) 一种基于密集光流法的云遮挡预测方法
Chang et al. Cloud tracking for solar irradiance prediction
CN111583298B (zh) 一种基于光流法的短时云图追踪方法
Dissawa et al. Cross-correlation based cloud motion estimation for short-term solar irradiation predictions
CN113936031A (zh) 一种基于机器视觉的云影轨迹预测方法
Zhang et al. Intrahour cloud tracking based on optical flow
El Jaouhari et al. Cloud tracking from whole-sky ground-based images
CN115984768A (zh) 一种基于固定型单目摄像头的多目标行人实时检测定位方法
CN107527003B (zh) 球型摄像机镜头粘附灰斑的视频质量诊断方法
CN103473787B (zh) 一种基于空间几何关系的桥梁上运动目标检测方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CB03 Change of inventor or designer information

Inventor after: Yan Yunfeng

Inventor after: Qi Donglian

Inventor after: Zhang Songjie

Inventor after: Dong Zhekang

Inventor after: Yang Xinyi

Inventor before: Qi Donglian

Inventor before: Yan Yunfeng

Inventor before: Zhang Songjie

Inventor before: Dong Zhekang

Inventor before: Yang Xinyi