【光学工程与电子技术】

基于色彩恒常理论的多光谱图像真彩色复原技术

张星铭1,孙文邦2,岳 广3

(1.空军航空大学战勤学院研究生队, 长春 130000;2.空军航空大学战勤学院, 长春 130000)

摘要:色彩匹配函数使得多光谱数据合成彩色图像简单易行,但是直接使用色彩匹配函数加权得到的多光谱彩色图像容易偏色,且无法避免低光照度和阴影对图像质量的影响;针对这一问题,本文结合色彩恒常理论,使用双sigmoid函数作为权重优化了Frankle-McCann Retinex算法在像素亮度计算的重置策略,达到对基于颜色匹配函数的真彩色复原技术改进,实现地物的彩色复原和增强;首先,对各波段图像建立亮度比较路径,通过对路径像素计算和迭代,估计目标在各波段的亮度值,其次利用色彩三刺激值原理与多光谱相机系统参数建立彩色图像三刺激分量生成模型,最后通过色彩空间转换得到目标真实色彩;利用Cave实验室数据和航空多光谱数据进行了验证和评价,结果表明本文提出的方法能够充分利用多光谱相机各波段图像信息,能够对图像阴影部分的目标色彩增强,色彩复原效果相对传统方法有明显提升,对于航空多光谱图像目视解译工作具有重要意义。

关键词:色彩恒常;色匹配函数;真彩色复原技术;多光谱

符合人眼视觉特性[1]的彩色遥感图像可以直观表达目标信息,对于目视解译工作十分重要。彩色相机获取的图像彩色自然,但是其不能提供地物目标的光谱特性,高光谱数据真彩色复原技术虽然较为成熟,但其庞大的数据量存在大量的冗余信息,且处理成本较高,航空多光谱相机兼具光谱分辨能力和空间分辨能力,具有灵活机动、数据轻便的优势,可以弥补星载数据在时间、地域及天气方面的劣势,可快速成像为解译人员服务。因此,研究提高航空多光谱遥感图像真彩色合成质量的方法非常必要。

遥感图像真彩色合成方法主要分为基于像素和基于光谱两类。基于像素的方法以假彩色合成、模拟真彩色合成[2]、直方图均衡等基于红、绿、蓝三通道的处理方法,这类方法简单直接,但是也使得大量的波段数据被舍弃,且容易出现图像偏色的问题。在这一方面,许辉熙[3]提出了波段模拟的方法,用于解决一些蓝波段缺失的多光谱相机真彩色合成的问题,韩秀珍[4]以风云三号D卫星MERSI-Ⅱ传感器数据为例,经过大气校正和非线性拉伸等方法生成了真彩色影像,取得较好的效果,但是都具有较强的针对性,对于灵活机动的航空平台不适用。

基于光谱的方法在考古、画作复原等方面得到广泛应用,近年来随着多光谱相机光谱精度的提升,该方法也逐渐被应用于遥感图像色还原领域。黄红莲等[5]提出了基于地面靶标和大气校正模型建立图像色彩转换矩阵的方法,通过地面靶标光谱反射率校正靶标图像色彩,并在同一大气条件下对多个地区的数据进行验证,取得较好的校正效果,但该方法受大气和光照条件影响,当外界条件变化时,结果也必然会改变,此外,其色彩复原的精度依赖于对地面靶标光谱反射率的测量,因此不适用于机动性强的航空平台多光谱载荷数据。兑利涛[6]研究了基于神经网络方法的光谱反射率重建,赵楠[7]研究了基于光谱三刺激值线性分段的色彩还原输出方法,实现对感兴趣区域灰度增强或抑制,MASAHIRO等[8]提出了一种基于16波段多光谱相机的多基色宽色域显示方法,阐述了谱段数据和输出设备色域之间的关系,这些方法在色彩质量上保持了一定的精度,但都是在实验室具备较高光谱分辨率和稳定平台的硬件条件下实现的,对于航空多光谱数据则具有较大的限制。

对以上文献研究发现,多光谱图像色彩还原大多基于光谱的重建,还原精度随着波段采样精度差异变化较大,且受到光源条件的制约。本文基于航空相机模型,在色匹配方法的基础上构建了图像对图像的解决方案,利用双sigmoid函数加权优化的色彩恒常理论算法对各波段原始图像进行针对性处理,实现地面目标色彩的复原与增强。该方法可以对暗处细节具有较好的增强作用,避免了光源光功率分布对计算过程的影响,对于航空多光谱数据彩色图像的目视解译工作十分有用。

1 物体色彩的计算原理

1.1 色彩匹配

颜色视觉理论[9]和等色实验(图1)证实了通过不同比例混合红、绿、蓝3种基色光的刺激值,可以得到人眼可见的绝大多数色彩。即,任意色刺激[C]可用式(1)计算:

[C]=R[R]+G[G]+B[B]

(1)

式中“=”表示视觉上相等。

图1 等色实验原理示意图

通过遍历全部可见光波段,便可以获取可见光范围内光谱的三刺激值,称其为色匹配函数[10](Color Matching Functions,简称CMFs)。

由于最初得到的色匹配函数存在大量负值,使用起来难于理解,CIE又推出了国际色度学系统CIE1931-XYZ系统,新的系统规定光谱三刺激值XYZ均为正值,(X)、(Z)两种原色只代表色度,没有亮度,光度量只与三刺激值Y成比例。该系统的优点在于物体的色彩在处理过程中可避免环境光照影响,进而有助于实现设备无关。两者可通过式(2)的转换矩阵互相转换[11],以此实现RGB色彩的获取。

(2)

CIE 1931 RGB系统和CIE 1931 XYZ系统的颜色匹配函数曲线如图2所示。

图2 CIE 1931 RGB色匹配函数曲线(a)和
CIE 1931 XYZ色匹配函数曲线(b)

1.2 基于色匹配函数的目标连续光谱色彩计算

根据三刺激值定义,假设某地物的反射光的光谱分布函数为φ(λ),光谱三刺激函数为φ(λ)按光谱三刺激值分解,则可以得到每一波长对应三刺激值,再对其从全可见光波段范围内积分,就可得到该颜色三刺激值,见式(3):

(3)

其中k是调整因子,调节k的值三刺激值也会发生变化,通常调整k值使得白光的Y值为100,所以令

波长λ的积分范围一般从380 nm至760 nm。

XYZ值代入式(2)即可获得目标反射光谱进入人眼的R、G、B刺激值,通过计算色品坐标[12]可以获取该色光在CIE色彩空间[13]中的位置。

1.3 存在问题

基于色匹配函数计算目标三刺激值的方法结果与光照条件密切关联,即当光谱分布函数φ(λ)变化时,三刺激值XYZ必然会随之发生变化,导致图像色彩偏离本色,影响航空侦察数据的解译结果。由于航空侦察具有跨地域和临机的特性,获取目标地点的光谱分布信息并不方便,因此需要对图像自身进行分析,使目标影像的色彩质量得到改善,最大限度符合其真实色彩。

2 基于色彩恒常理论真彩色复原

人眼能够在不同环境下识别出同一种颜色,但显示设备做不到,低光照度[14]和阴影[15]条件都会使图像的色彩发生变化。色彩恒常理论认为物体的颜色是其自身对光的反射特性决定的[16],并不会随光照条件而改变,通过Retinex算法可以实现对图像自身光照估计,并消除光照不均匀[17],Frankle-McCann Retinex(FMR)算法[17]是一种基于路径的多重迭代策略亮度估计方法,如图3所示。图像中单个点的像素值取决于一条特定路径的环绕的结果,通过对螺旋式路径上的各像素点的灰度值进行多次求差、重置、平均,通过迭代使结果逼近理想值。路径上越靠近目标点采样的点数越多,保证了光照估计的结果符合实际[18]。这正是航空垂直图像彩色复原希望实现的目标。

图3 FMR算法原理示意图

2.1 FMR算法的亮度估计

将多光谱所有通道的图像数据进行对数变换,得到对数域S(x,y)。初始化常数图像矩阵r(x,y)。该矩阵作为进行迭代运算的初始值,本文选取各波段数据亮度最大值为初始值。计算路径。设(rows,cols)为图像的行列数,如图3所示,利用式(4)和式(5)可求得选取距离为D的目标点沿xy方向的最远坐标(0,D)和(D,0),为计算方便选取y轴正方向最远点为起始点顺时针依次计算。

P=fix[log2(min(rows,cols))-1]

(4)

D=2P

(5)

之后各采样点与中心点的距离缩短为上一步的一半,同时按顺时针改变路径旋转方向,即D=-D/2。则可得到沿路径排列的一组点集A1,A2,…,Am

对路径上的像素点亮度计算:中心点像素估计值等于当前像素值与路径上的像素值的差的一半与前一时刻的估计值之间的和,见式(6)、式(7)。要注意的是,采样点按照从远到近的顺序参与计算。

r(x,y)为目标点初始值,rt(x,y)为第t次估计的结果,ΔL为路径上的亮度差,max是原图像的像素最大值,则第t+1次亮度估计的结果为

(6)

(7)

FMR算法亮度估计流程见图4。

图4 FMR算法亮度估计流程框图

2.2 双sigmoid函数的最大值优化

实验表明:式(7)中当rt(x,y)的取值大于max时,取不同的重置值可以得到不同的效果,当取值为原图像素最大值max时,原图像亮部像素亮度增强幅度较大,导致色彩饱和度降低,见图5(b);当取值为均值时,图像对比度下降,图像质量较低,见图5(c);当取原像素值时,亮部分的色泽可以得到保留,但是对于暗处提升效果不明显,见图5(d)。

图5 不同重置值的效果

为达到增强图像暗区域亮度同时抑制亮区域增益的目的,本文提出了基于双sigmoid函数的加权方法对式(7)改进,改进后的函数曲线如图6所示。

图6 双sigmoid函数曲线

(8)

(9)

X轴表示图像灰度gray量化范围,函数sig1的值为亮度最大值max的权重,sig2的值为原像素值Pic(x,y)的权重,则式(7)可改为式(10)所示:

(10)

2.3 多光谱相机三刺激分量模型

对于航空多光谱相机的成像系统,设光源光谱分布为I(λ),物体的光谱反射率为ρ(λ),成像系统的光谱透过率函数τ1(λ),光学系统的光谱透过率函数为τ2(λ),CCD的光谱响应函数为S(λ),由式(3)可得到相机响应的三刺激分量,见式(11)。

(11)

其中,φ(λ)=I(λ)ρ(λ)τ1(λ)τ2(λ)S(λ)。

多光谱相机的分光系统使得接收到的自然光被分为若干透射波谱,不同的谱段在CCD器件上作用,得到不同波段的图像数据。

所以,上式积分过程可描述为在各个透射谱段内根据波长积分,在不同谱段间求和的过程,见式(12)。

(12)

假设多光谱成像系统是线性的,那么φ(λ)在物理意义上构成以λ为变量的多光谱成像系统的响应,其在不同透射带宽内积分的结果是该波段的图像灰度值Hl(i, j),见式(13)。

Hl(i, j)=φ(λ)dλ

(13)

其中,(i, j)表示像素在图像的位置,l为波段序号。

通过光谱分布计算地物三刺激值的过程就可以通过式(13)转化为对多光谱图像灰度的加权,见式(14),其中,n为波段数,权重系数即为对应波段的色匹配函数值:

(14)

2.4 复原步骤

1) 利用FMR算法对各波段图像计算,最大迭代次数设为4次,得到过渡图像FMR_Hl(i, j)。

2) 确定透射波段的光谱三刺激值通过10 nm间隔三刺激值加权表进行拟合后读取对应波段的数值,构成新的加权表。

3) 将以上数据代入式(14)中,计算图像三刺激分量,结果见式(15)。

(15)

其中,BminBmax分别表示第l波段的近端和远端。

4) 利用转换矩阵将步骤3)的结果转换成RGB三刺激值对应的图像三刺激分量。

5) 归一化。

图7所示利用本文方法对多光谱数据真彩色还原的技术流程。

图7 多光谱图像真彩色还原技术流程框图

3 实验分析

遥感影像质量评价方法分为主观评价法和客观评价法两类,主观评价方法是根据人的主观经验评价影像质量,是最直接的一种影像质量评价方法[19]。客观评价方法采用不同的指标组合来对影像进行全面的定量评价,具有可批量处理、结果可重现的特点,不会出现人为偏差。

实验使用哥伦比亚大学Cave实验室的一组暗背景、低照度的多光谱数据样本和一组实际获取的滤光片式航空多光谱相机地面数据进行分析评价,分别使用假彩色合成法、CMFs方法和本文方法对数据操作对比,其中假彩色合成选取红色中心波长700 nm、绿色选取中心波长540 nm、蓝色选取中心波长450 nm的数据。对结果采用目视观察和计算图像质量评价指标两种方法进行比较。

Cave实验室的3个样本为暗背景上的羽毛、玻璃和小球,对比数据为实验提供的彩色图像;航空数据样本为中科院的滤光片阵列式多光谱相机,可见光部分8个波段,公路、房屋、林地、土地和阴影处的雪,对比数据为非同时相的Google彩色影像。

Cave实验室样本处理结果及直方图见图8、图9。航空数据样本处理结果及直方图见图10、图11。

图8 Cave实验室数据样本彩色复原

图9 Cave实验室数据样本直方图

图10 航空遥感影像彩色复原

图11 航空遥感影像彩色复原直方图

通过对图像和直方图的主观判断,假彩色合成图像效果存在较为严重的偏色问题,CMFs方法可消除部分色彩偏差,但缺乏对暗处目标的亮度和色彩的进一步改善。本文方法较好的再现了地面目标的原有色彩,并且使得暗处的景物得到显现。

图像的客观评价指标通常基于图像的特征进行计算,包括信息熵、清晰度、信噪比、均方差等[20]。当前航空遥感领域对图像色彩质量评定缺少统一的标准,文献[21]研究了一种具有普遍适用性的等效圆的检测偏色的方法,该方法基于图像自身特征进行分析,可用于快速评价图像色彩质量。本文在选取3种常用的评价指标的基础上加入偏色度、总色度差、饱和度差三项彩色图像评价指标进行检验。

图像的信息熵是用来描述图像信息量的重要指标,计算方法见式(16):

(16)

式(16)中:pi表示像素点灰度值为i的概率,信息熵越大,表示纹理的复杂度越高,图像信息量越丰富。

图像亮度均值反映了图像的平均亮度,均值越大说明图像亮度越大,反之越小,计算方法见式(17):

(17)

标准差反映了图像像素值与均值的离散程度,标准差越大说明图像的质量越好。计算方法见式(18):

(18)

CIELab是均匀色空间,符合色彩之间的差别符合人眼的观察,因此通过计算两点的距离即可实现色彩之间差异的比较,其色差计算公式见式(19)、式(20),偏色度公式见式(21)。

总色差:

(19)

饱和度差:

(20)

偏色度:利用ab二维直方图计算图像色度中心与原点(0,0)的距离与等效圆半径的比值,用K表示,K值越大,偏色越严重。

K=D/M

(21)

D为色度ab的均值到原点的距离。M为等效圆半径。

其中为图像ab色调的平均值。

由于航空多光谱图像的对比图像来自Google Earth的非同时段拍摄数据,地物的光照度、色泽等均有一定的差异,因此色差值并不代表图像与真实地物的色彩差异,但可作为图像质量恢复效果的参考依据。图像质量评价指标计算结果见表1(对应图8、图9)、表2(对应图10、图11)。

为了使数据表达更加直观,将指标数据通过柱状图的形式显示,如图12、图13。

表1 Cave实验室样品图像评价指标

评价指标羽毛玻璃小球信息熵假彩色5.566 86.702 56.168 6CMFs5.211 46.434 84.253 6本文7.733 57.733 56.756 9标准差假彩色27.016 534.327 546.476 9CMFs23.664 128.456 627.814 6本文45.571 058.365 443.710 4均值假彩色23.158 746.237 937.703 0CMFs25.845 743.567 215.390 3本文106.036 9106.036 970.559 9偏色度假彩色0.713 3030.516 5210.781 256CMFs0.701 4060.605 7840.479 520本文0.335 6210.542 9250.432 908饱和度差假彩色-48.871 7-11.005 1-16.989 8CMFs-25.032 8-26.776 7-25.166 4本文-9.868 8-10.510 2-3.150 0总色差假彩色75.582 0 19.495 6 29.637 1CMFs38.677 233.953 530.702 0本文23.328 213.476 210.236 4

表2 航空多光谱图像评价指标

评价指标公路土地房屋林地阴影中的雪信息熵假彩色6.822 46.491 47.702 46.693 06.936 2CMFs6.807 36.259 97.631 06.623 06.919 0本文7.109 26.782 97.757 57.418 57.383 6标准差假彩色51.117 132.814 056.002 430.070 851.048 5CMFs49.687 129.722 751.435 227.887 149.774 8本文48.897 529.555 256.126 542.784 153.688 2均值假彩色66.793 185.307 6124.767 752.203 362.292 0CMFs65.993 380.902 3123.522 952.308 363.798 1本文86.063 7120.773 4139.337 2101.147 790.494 5

续表(表2)

评价指标公路土地房屋林地阴影中的雪偏色度假彩色0.211 180.950 50.217 60.121 31.039 1 CMFs0.233 610.799 40.169 80.188 31.058 7本文0.308 290.602 40.102 30.174 90.668 7饱和度差假彩色20.533 925.330 529.976 911.140 87.724 5CMFs11.947 034.163 144.904 37.022 531.749 4本文5.335 9-0.300 538.778 81.575 65.115 9总色差假彩色63.534 214.825 836.402 258.943 452.549 1CMFs19.471 912.577 760.792 912.469 165.693 4本文33.107 715.275 058.943 420.376 435.736 7

图12 Cave样本彩色复原评价指标

图13 航空遥感影像彩色复原评价指标

由以上数据分析,本文方法比假彩色的信息熵分别提高4.2%、4.5%、0.7%、10.8%、6.5%,标准差分别提升-4.3%、-9.9%、0.2%、42.3%、5.2%,灰度均值分别提升28.9%、41.6%、11.7%、93.8%、45.3%。通过柱状图观察可知,本文方法的结果自身偏色度较低,与对比图像之间的色彩饱和度差比较小,由于本文方法对于图像亮度的变化比较剧烈,所以总色差的变化比较波动,但相较于其他两种方法,其数值仍可以保持持平或降低。

4 结论

本文研究并改进了一种多光谱图像真彩色复原的方法,该方法结合色彩视觉原理、色彩恒常理论,设计了一组双sigmoid权重函数,利用其改进FMR算法的重置值,解决因重置值设置过大导致在图像复原时亮部色度损失严重的问题,然后将该方法与色匹配函数相结合,利用多光谱各波段数据求得三原色分量,实现航空多光谱遥感图像色彩复原。

通过对实验室数据和航空遥感数据处理表明,该方法较好的利用多光谱各波段段图像数据信息,色彩复原效果相对传统算法有较高提升,可以较好地还原地物真实色彩。本文方法对图像低光照度和阴影处的目标具有增强作用,提高了图像的可读性。本文方法仍需要进一步优化,例如可以进一步修复图像的饱和度要素,达到色彩更加鲜明的结果。

参考文献:

[1] 许向阳.与视觉认知过程相关的图像色貌建模的研究[D].广州:华南理工大学,2016.

[2] 王忠武,刘顺喜.资源一号02C卫星多光谱数据真彩色模拟技术研究[J].测绘与空间地理信息,2014,37(01):13-16.

[3] 许辉熙,陈云浩,薛万蓉.蓝波段缺失遥感影像真彩色模拟方法研究[J].激光与光电子学进展,2015,52(05):77-84.

[4] 韩秀珍,王峰,单天婵.风云三号D星真彩色影像合成方法研究及应用[J].海洋气象学报,2019,39(02):13-23.

[5] 黄红莲,易维宁,杜丽丽,等.基于人工靶标的多光谱遥感图像真彩色合成[J].红外与激光工程,2016,45(11):354-359.

[6] 兑利涛.多光谱图像的光谱反射率重建技术研究[D].南昌:华东交通大学,2017.

[7] 赵楠.多光谱图像色彩还原方法及在成像光谱仪中的应用[D].武汉:华中科技大学,2011.

[8] YAMAGUCHI M,TERAJI T,OHSAWA K.et al.Color image reproduction based on the multispectral and multiprimary imaging:Experimental Evaluation[C]//Color Imaging:Device-Independent Color,Color Hardcopy,and Applications Ⅶ.San Jose,CA(US):SPIE,2002:15-26.

[9] 黄敏,史春洁,李泽阳,等.色觉正常观察者辨色差异影响研究[J].光学学报,2016,36(09):341-349.

[10] 黄敏,何瑞丽,郭春丽,等.不同年龄观察者颜色匹配函数的测试及优化[J].光学学报,2018,38(03):459-466.

[11] 吕晓凯,李传荣,马灵玲,等.一种基于物理机理的高光谱遥感图像真彩色校正模型[J].遥感技术与应用,2013,28(03):467-473.

[12] WESTLAND S.CIE Chromaticity Coordinates(xyY)[J].Encyclopedia of Color Science and Technology,2015(02):1-4.

[13] 金伟其,胡威捷.辐射度、光度与色度及其测量[ M].北京:北京理工大学出版社,2006.

[14] 贾存坤,戴声奎,卫志敏.采用亮通道先验的低照度图像增强算法[J].华侨大学学报(自然科学版),2018,39(04):595-599.

[15] 杨玲,阮心玲,李畅.一种自适应Retinex的航空影像阴影消除方法[J].测绘工程,2013,22(03):1-4,15.

[16] 靳珍璐,潘泉,赵春晖,等.基于局部精确直方图匹配的无人机景象匹配导航色彩恒常算法[J].中国惯性技术学报,2015,23(05):674-680.

[17] 刘军.基于Retinex理论的彩色图像增强技术研究[D].西安:中国科学院研究生院(西安光学精密机械研究所),2015.

[18] 秦华锋,刘霞.基于稀疏自编码的手指静脉图像分割[J].重庆工商大学学报(自然科学版),2019,36(4):1-8.

[19] 严亮.基于Retinex理论的图像增强技术研究[D].武汉:华中科技大学,2016.

[20] 董胜越,孙根云,杜永明,等.高分五号全谱段光谱成像仪影像数据质量评价研究[J].遥感技术与应用,2020,35(02):381-388.

[21] 尹灵芝,朱军,蔡国林,等.遥感影像质量评价方法综述[J].测绘与空间地理信息,2014,37(12):32-35,45.

[22] 董玲.基于视觉感知特征的图像质量评价方法研究[D].武汉:武汉理工大学,2017.

True Color Restoration of Multispectral Image Based on Color Constancy

ZHANG Xingming1, SUN Wenbang2, YUE Guang1

(1.Graduate team of School of Aviation Operations and Services, Air Force Aviation University, Changchun 130000, China; 2.School of Aviation Operations and Services, Air Force Aviation University, Changchun 130000, China)

Abstract: The color matching function makes multi-spectral color restoration simple and easy, but the multi-spectral color image that is directly weighted by the color matching function is prone to color cast, and the impact of low light and shadow on the image quality cannot be avoided. In response to this problem, this paper combined the color constancy theory and used the dual sigmoid function as the weight to optimize the reset strategy of the Frankle-McCann Retinex algorithm in pixel brightness calculation, to achieve the improvement of the true color restoration technology based on the color matching function, and to achieve the ground object. Firstly, we established a brightness comparison path for each band image, and estimated the brightness value of the target in each band by calculating and iterating the path pixels. Secondly, we used the color tristimulus value principle and the multispectral camera system parameters to establish a color image tristimulus component generation model. Finally, the real color of the target was obtained through color space conversion. The experiment used Cave laboratory data and aerial multispectral data to verify and evaluate. The results show that the method proposed in this paper can make full use of the image information of each band of the multispectral camera, and can enhance the target color in the shadow of the image, and the color restoration effect is relatively traditional methods have been significantly improved, which is of great significance to the visual interpretation of aerial multispectral images.

Key words: color constant; color matching function;true color restoration technology; multispectral

doi: 10.11809/bqzbgcxb2020.11.045

收稿日期:2020-05-11; 修回日期:2020-06-20

基金项目:全军军事类研究生课题(DSSQ910252018010)

作者简介:张星铭(1987—),男,在职研究生,主要从事图像处理研究,E-mail:1055975729@qq.com。

通讯作者:孙文邦(1976—),男,博士,副教授,主要从事图像处理研究,E-mail:chwenbang@163.com。

本文引用格式: 张星铭,孙文邦,岳广.基于色彩恒常理论的多光谱图像真彩色复原技术[J].兵器装备工程学报,2020,41(11):248-256.

Citation format:ZHANG Xingming, SUN Wenbang, YUE Guang.True Color Restoration of Multispectral Image Based on Color Constancy[J].Journal of Ordnance Equipment Engineering,2020,41(11):248-256.

中图分类号:TP751

文献标识码:A

文章编号:2096-2304(2020)11-0248-09

科学编辑 刘丹凤 博士(大连民族大学讲师)

责任编辑 杨梅梅