一种用于克服动态问题影响的高动态范围图像融合方法

文档序号:10572957阅读:336来源:国知局
一种用于克服动态问题影响的高动态范围图像融合方法
【专利摘要】一种用于克服动态问题影响的高动态范围图像融合方法,用于解决拍摄场景动态以及成像系统晃动影响的高动态范围的图像融合问题。其中,采用基于亮度映射函数标定相机响应函数的方法克服动态问题给相机响应函数标定带来的影响,然后结合总体最小二乘法引入校正矩阵,对误差进行联合补偿,从而无需高精度的图像配准即能够获得成像系统的相机响应函数;进而通过分析多曝光图像序列特点选择权重函数,结合标定的相机响应函数构建高动态范围图像权重融合函数,利用输入的多曝光图像序列获得高动态范围图像融合结果。本发明方法能够克服动态问题给相机响应函数标定带来的困难,进而提高了多曝光高动态范围融合技术的适用性。
【专利说明】
一种用于克服动态问题影响的高动态范围图像融合方法
技术领域
[0001] 本发明属于图像处理技术领域,涉及一种适用于解决拍摄场景动态以及成像系统 晃动影响的高动态范围图像融合方法。
【背景技术】
[0002] 现实环境及战场环境中往往会出现动态范围较大的场景,即同一视场中同时出现 较亮和较暗的信息。如果拍摄场景的动态范围超过了成像器件的动态范围,就会造成采集 的图像出现过饱和或欠饱和的情况,从而影响目标信息的采集及识别。如发射的导弹由于 尾焰过亮会造成采集图像饱和,从而无法看清弹体等现象,影响目标的分析。因此,为了便 于目标观测,对观测设备、发射时间窗口选取以及天气环境等因素均提出了较高的要求。
[0003] 高动态范围成像技术能够避免因为拍摄方向(如逆光情况)和曝光量不足而使图 像存在壳度?目息缺失和偏色等问题,不影响真实场景的?目息米集,有利于复杂环境下获得 更高的成像质量,因此被广泛应用于模式识别、智能交通系统、深空探测、视频监控跟踪、机 器人视觉导航、工业检测、遥感遥测、军事监视侦察以及武器导航制导等众多领域,具有重 要的研究价值。
[0004] 成像系统的相机响应函数能够建立拍摄图像强度信息与场景辐亮度之间的严格 映射关系,是高动态范围图像融合的关键技术,拍摄场景的动态以及成像系统的晃动也会 影响成像系统相机响应函数的标定,而多曝光图像序列配准很难达到较高的精度,从而给 动态目标高精度的高动态范围图像融合带来较大困难。
[0005] 目前已有的高动态范围图像融合方法主要是利用不同曝光量的已配准的图像序 列作为输入,在已知图像之间的曝光比例的情况下,建立了图像各通道图像强度信息与场 景辐照度之间的映射关系,从而获得反映拍摄场景的高动态范围图像,但通常需要拍摄场 景保持静态并假设图像序列照度一致。例如2015年计算机工程与科学上的《基于相机响应 曲线的高动态范围图像融合》一文,文中提出利用已配准的多曝光图像序列每个像素点的 信息标定相机响应函数,进而利用相机响应函数建立图像亮度值与拍摄场景辐照度之间的 映射关系,从而获得能够反映拍摄场景高动态范围的辐照度图像。该算法需要精确配准的 静态场景多曝光量图像序列作为输入进行相机响应函数的标定,当拍摄场景存在动态问题 影响或相机晃动的情况下会给图像配准带来较大困难,从而无法进行相机响应函数标定, 也即该算法无法克服动态问题带来的影响;而且,利用该算法获得的高动态范围图像融合 存在偏色,需要结合相应的色彩校正算法对输出结果进行白平衡。

【发明内容】

[0006] 本发明解决的技术问题是:克服现有技术的不足,提供了一种克服动态问题影响 的高动态范围图像融合方法,能够利用彩色图像各颜色通道的亮度统计特征获得亮度映射 函数,从而克服动态问题影响完成相机响应函数的标定,进而结合输入图像序列特征构建 权重融合函数进行高动态范围的图像融合。
[0007] 本发明的技术解决方案是:一种用于克服动态问题影响的高动态范围图像融合方 法,包括如下步骤:
[0008] (1)利用相机获得目标不同曝光量的曝光图像序列;
[0009] (2)根据目标的曝光图像序列,获取图像的亮度映射函数T(m),具体为:
[0010] (21)将曝光图像序列按曝光量从低到高进行排列,然后从第一帧图像开始每两帧 图像组成一个输入图像对;
[0011] (22)针对每一个输入图像对,获得该输入图像对中两帧图像各颜色通道的亮度信 息的累积直方图;
[0012] (23)利用输入图像对的两帧图像各颜色通道亮度信息的累积直方图,计算两帧图 像各颜色通道之间的亮度映射函数T(m),
[0013]
[0014] 其中m为亮度值,HA(m)为输入图像对中的第一帧图像的某颜色通道亮度值为m的 累积直方图函数,/Y^(W)为输入图像对中的第二帧图像的对应颜色通道亮度值为m的累积 直方图逆函数;
[0015] (3)利用亮度映射函数T(m)对相机响应函数进行标定,获得相机响应逆函数的近 似表达式g;
[0016] (4)结合相机响应逆函数的近似表达式g及曝光图像序列曝光时间信息△ t构建图 像权重融合函数,获得彩色图像各颜色通道的辐照度信息,从而获得反映拍摄场景信息 的高动态范围图像。
[0017] 所述的图像权重融合函数^为:
[0018]
[0019] 式中i为图像中像素点的空间索引,j为曝光图像序列中的图像索引,Mij为输入曝 光图像序列中第j帧图像上像素点i的亮度信息,w为三角帽权值函数,为曝光图像序列 中第j帧图像的曝光时间,P为曝光图像序列中的图像总帧数。
[0020] 所述步骤(3)中获得相机响应逆函数的近似表达式g的方法为:
[0021] (31)获取两帧不同曝光量图像对应关系如下:
[0022] g(T(MA)) =kg(MA)
[0023] 其中Μα图像A上被估计的像素点的图像亮度信息,T是图像B和图像A之间的亮度传 递函数,k是两帧不同曝光量图像中图像A和图像B之间的曝光比;
[0024] (32)加入权重函数,将步骤(31)的结果修改为:
[0025]
[0026] 其中Ca是第一帧图像A中亮度值为Μα的像素个数,Cb是第二帧图像B中亮度值为Mb 的像素个数;
[0027] (33)利用多项式对相机响应逆函数进行建模,将不同输入图像对获得的亮度传递 函数以及相机响应逆函数多项式代入步骤(32)的表达式中获得一组超定方程;
[0028] (34)利用总体最小二乘法求解所述的超定方程,获得相机响应逆函数的模型系 数,由此得到相机响应逆函数的近似表达式g。
[0029] 本发明与现有技术相比的优点在于:
[0030] (1)本发明根据图像亮度信息的统计特征,提出了基于累积直方图建模的亮度映 射函数计算方法,进而利用亮度映射函数而不是图像各像素点信息进行相机响应函数标 定,从而克服了拍摄场景动态问题以及成像系统的晃动给相机响应函数标定带来的困难, 同时降低了相机响应函数标定算法方计算量,提高了计算效率;
[0031] (2)本发明利用总体最小二乘法求解相机响应逆函数标定模型系数,由于总体最 小二乘法引入校正矩阵对存在的误差进行联合补偿,抑制了误差及噪声对矩阵方程求解的 影响,能够获得更高精度的相机响应函数标定结果;
[0032] (3)本发明利用三角帽函数作为权重,结合标定的相机响应函数共同构建高动态 范围融合权重函数,可有效的提取输入不同曝光量图像序列中的有效信息,从而有效避免 高动态范围图像融合产生偏色的问题;
[0033] (4)本发明利用亮度映射函数标定的相机响应函数进行高动态图像融合,由于利 用的是图像亮度值的统计特征而不是各像素点的信息,因此能够消除高精度相机响应函数 标定对输入图像序列绝对空间一致的要求,从而克服场景动态问题以及成像系统自身晃动 带来的影响,提升了高动态范围图像融合技术的应用范围。
【附图说明】
[0034]图1为本发明方法的流程框图;
[0035]图2为本发明拍摄的存在相机晃动的多曝光图像序列;
[0036]图3为本发明利用不同输入图像对获得的亮度映射函数;
[0037] 图4为本发明标定的相机响应逆函数结果;
[0038] 图5为本发明高动态范围图像融合结果。
【具体实施方式】
[0039] 本发明方法的基本思想是:由于拍摄目标光照条件不变的条件下,针对拍摄场景 中某一点,比其亮与比其暗的像素点个数比例是恒定的,基本不受成像系统晃动的影响,即 场景的动态问题影响不会明显的改变图像亮度信息的统计特征,场景辐亮度的总体直方图 基本保持恒定。所以可利用不同曝光量的图像通过累积直方图规则化对输入图像对进行直 方图建模以获得图像之间的亮度映射函数,然后结合总体最小二乘法引入校正矩阵对误差 进行联合补偿,从而能够利用亮度映射函数而非采集序列图像每个像素点的信息实现成像 系统的高精度相机响应函数标定,进而分析输入的多曝光图像序列特征,选择合适的权重 函数结合标定的相机响应逆函数构建高动态范围图像融合权重函数,最终获得拍摄目标的 高动态范围图像。这种处理方式能够消除高精度相机响应函数标定对输入图像序列绝对空 间一致的要求,克服拍摄目标的动态问题以及成像系统的晃动带来的影响,同时降低计算 量,扩展高动态范围图像融合的应用范围。
[0040] 如图1所示,为本发明方法的流程框图,主要包括如下几个步骤:
[0041 ] (1)根据拍摄目标特点,获得目标多曝光图像序列;
[0042] 首先利用相机根据拍摄场景目标需求,进行相应的不同曝光量设置,采集目标的 多曝光图像序列,曝光量由低到高为...,EP。
[0043] 假设f表示相机响应函数,利用相机响应函数f表示场景辐照度值I对应采集图像 上亮度值Μ的关系,如下所示:
[0044] M = f(I)
[0045] 根据相机响应函数的单调性,存在反函数8 = = =
[0046] 针对同一场景拍摄的两帧不同曝光量图像对应辐照度之间的映射关系如下所示:
[0047] g(mB)=kg(mA)
[0048] 式中mA是图像A上被估计的像素点的图像亮度信息,mB是图像B上对应点的亮度信 息,k是图像A和图像B之间的曝光比。
[0049] (2)根据目标多曝光序列图像,获取图像的亮度映射函数。
[0050] 该步骤又包括以下几个步骤:
[0051] (21)根据曝光量将图像序列分为每两帧图像组成的多组输入图像对。
[0052] 将输入图像序列按曝光量从低到高进行排列,即E^Es,...,EP,由第一帧图像开始 依次每两帧图像组成输入图像对(Ei,E 2),. . .,(En,EP)。
[0053] (22)针对每一组输入图像对,获得该两帧图像各颜色通道亮度统计特征信息; [0054]适当的目标以及场景运动的情况下,场景辐亮度的总体直方图基本保持恒定,因 此可利用不同曝光的图像通过累积直方图规则化对输入图像进行直方图建模,获得图像对 之间的亮度映射函数。
[0055] 假设输入图像的各颜色通道亮度信息对应的直方图没有空缺部分,彩色图像的各 颜色通道(分别为红绿蓝)的累积直方图Hr(M)、Hg(M)及Hb(M)如下所示:
[0056]
[0057]式中Μ是图像各颜色通道累积直方图的最高亮度值,m是累积直方图对应亮度范围 中的亮度信息,hR(m)、hc(m)以及hB(m)分别为红绿蓝三个颜色通道的直方图。
[0058] (23)利用输入图像对的各颜色通道亮度信息的累积直方图规则化计算两帧图像 各颜色通道之间的亮度映射函数。
[0059] 由于拍摄目标光照条件不变的条件下,针对拍摄场景中某一点,比其亮与比其暗 的像素点个数比例是恒定的,基本不受成像系统晃动的影响。根据该性质假设输入的两帧 图像是在没有大幅度相机晃动或场景变化的情况下获得的,以图像红色通道为例,利用红 色通道的亮度映射函数T建立输入图像对中的第二帧图像红色通道的亮度值Mb与第一帧图 像中红色通道的亮度值Μα的对应关系,即Mb = T(Ma)。由于它们对应于拍摄场景中相同的点, 第一帧图像中亮度值小于Μα的像素点数目与第二帧图像中亮度值小于Mb的像素点数目应相 等,可利用累积直方图表示,即Hra(Ma) =Hrb(Mb) =Hrb(T(Ma))。令MA=m,则:
[0060]
[0061 ]式中Tr为两帧图像之间红色通道的亮度映射函数,Hra为第一帧图像红色通道的累 积直方图,Hrb为第二帧图像红色通道的累积直方图。
[0062]彩色图像各颜色通道的亮度映射函数T可通过上式结合输入图像的累积直方图获 得。
[0063] (3)利用亮度映射函数对成像系统的相机响应函数进行标定,建立成像系统的相 机响应逆函数标定模型;利用总体最小二乘法对相机响应逆函数标定模型系数进行求解, 从而获得成像系统的相机响应逆函数。具体又分为如下步骤:
[0064] (31)将获得的表述两帧图像对应像素点亮度值之间关系的亮度映射函数Mb = T (Ma) : =g4(kg(MA))带入g(mB)=kg(mA),得到两帧不同曝光量图像对应关系如下:
[0065] g(T(MA)) =kg(MA)
[0066] (32)加入权重函数。
[0067]用于计算图像之间亮度传递函数的两帧图像的累积直方图Ha和Hb并不是一致的, 这意味着任意一个颜色通道的亮度值对应的亮度传递函数依赖于图像中该亮度值对应的 像素数目,因此需要对约束方程g(T(M A)) =kg(MA)加入相应权重以提高相机响应函数标定 精度并克服相机响应函数的分形自相似性。
[0068] 对于映射对(Ma,T(Ma))而言,hA(MA)=CA是第一帧图像A中亮度值为Μα的像素个数, hB(Mb) = hB(T(Ma) ) = CB是第二帧图像Β中亮度值为Mb的像素个数。为了求解后续的最小二乘 问题,将不同曝光量图像结合相机响应函数及相应的约束权重得到:
[0069]
[0070] 由于直接在方程两端乘以相应的直方图权重,过大的权重差异会造成方程结果震 荡,给结果带来较大的误差,所以对直方图进行开根号处理降低两端权重差异,从而获得稳 定的求解结果。
[0071] (33)建立相机响应逆函数超定方程,并利用总体最小二乘法求解相机响应逆函 数。
[0072] 可以利用六阶多项式对相机响应函数进行建模,如下所示:
[0073] g=m〇+mix+m2X2+m3X3+m4X4+m5X 5+m6X6
[0074] 将不同输入图像对获得的亮度传递函数结合相机响应函数六阶多项式模型带入 另
庆得一组超定方程如下所示:
[0075]
[0076] 由于总体最小二乘法相较于最小二乘法能够实现对误差的联合补偿,从而能够抑 制误差及噪声对矩阵方程求解的影响,从而精确地求解系数{πιο,πη,...,m 6},获得相机响应 逆函数的近似表达式g。
[0077] (331)相机响应逆函数超定方程模型线性化。
[0078] 利用变量代换法将相机响应逆函数标定的超定方程模型进行简化处理,如下所 示:
[0079]
[0080]
[0081]
[0082] (332)获得零均值的数据矩阵。
[0083]
[0084] 式中
A为构建的 零均值数据矩阵。
[0085] (333)求解相机响应逆函数模型系数。
[0086] 对ATA进行特征值分解如下所示:
[0087]
[0088] 式中〇〇>〇1>.">〇6为特征值,Ul,U2,…,U6为特征值对应的特征向量。
[0089] 则相机响应逆函数模型系数[mo mi ... m6]T向量取特征值〇6对应的特征向量,即:
[0090] [m0 mi · · · m6]T=U6T
[0091 ]将此参数进行代入刘姐相机响应逆函数模型,从而实现相机响应逆函数的标定。
[0092] (4)利用标定的成像系统的相机响应逆函数建立彩色图像各颜色通道亮度信息与 拍摄场景辐照度信息之间的映射关系,针对多曝光图像序列特点构建高动态范围图像权重 融合函数进行高动态范围图像融合。
[0093] 相机响应逆函数是阐述亮度信息与能够反映真实拍摄信息的辐照度之间的映射 关系的。通过将彩色图像分解为不同的颜色通道,分别对不同的颜色通道中的亮度信息通 过相应的相机响应逆函数进行映射。对于动态范围较大的拍摄场景,一帧低动态范围图像 无法完全反映拍摄场景包含的信息,但通过调节成像系统的曝光量,可获得多曝光图像序 列,每一帧不同曝光量的图像中包含不同的动态范围信息。所以可通过利用相机响应逆函 数构建高动态范围融合权重函数结合输入多曝光图像序列实现高动态范围图像融合。 [0094] (41)选择权重函数。
[0095] 高动态范围图像融合中权重函数的选取应弱化欠曝光以及过饱和部分包含信息 量少的区域,强化曝光量合适的部分,即能够有效表达图像信息的部分,为此采用三角帽函 数作为权重函数:
[0096]
[0097]式中w(M)为亮度值Μ的三角帽权重函数值,Mmin为输入图像特定颜色通道最小亮度 值,Mmax为输入图像同一特定颜色通道最大亮度值。
[0098] (42)构建高动态范围图像融合权重函数。
[0099]针对彩色图像序列,将序列中的每一帧图像分解为R、G、B三个颜色通道,分别对各 颜色通道的亮度信息进行处理。由于相机响应逆函数能够获得亮度信息与辐照度信息之间 的映射关系,但每一帧图像的亮度信息由于其曝光量不同,所包含的有效信息不同,需要利 用选择的权重函数对每一帧图像中各像素点的亮度信息进行加权融合,同时考虑各帧图像 不同曝光量的影响,构建高动态范围图像融合权重函数如下所示:
[0100]
[0101]式中i为图像中像素点的空间索引,j为图像序列中的图像索引,g为标定的相机响 应逆函数,Mi j为输入多曝光图像序列中第j帧图像上像素点i的亮度信息,w为三角帽权值函 数,Atj为图像序列中第j帧图像曝光时间。
[0102]上式即为构建的高动态范围图像权重融合函数,通过多曝光图像序列中对应空间 位置像素点的亮度信息获得的相应权重以及相机响应逆函数映射值,结合相应输入帧图像 的曝光时间代入上式,即能够获得拍摄场景对应点的辐照度信息。
[0103 ] (43)获得拍摄场景的高动态范围图像。
[0104] 获得彩色图像R、G、B各颜色通道的亮度信息,结合标定的相机响应逆函数以及构 建的高动态范围图像权重融合函数,针对每一个颜色通道亮度信息进行逆映射,获得相应 颜色通道的辐照度信息,即能够获得通过单帧图像反映拍摄场景丰富信息的高动态范围图 像。
[0105] (431)获得彩色图像各颜色通道的亮度信息。
[0106] (432)通过多曝光输入图像序列中对应空间位置像素点的亮度信息获得的相应权 重以及相机响应逆函数映射值,结合相应输入帧图像的曝光时间代入(42)中构建的高动态 范围图像权重融合函数,即能够获得拍摄场景对应点的辐照度信息。
[0107] (433)分别获得各颜色通道的辐照度信息,即获得需要的反映拍摄场景高动态范 围的RGB颜色模型高动态范围图像。
[0108] 实施例
[0109] 下面以一个具体实例对本发明方法的工作原理和过程作进一步说明。
[0110] (1)获得目标多曝光图像序列。利用佳能D90相机获得的一组含有成像系统晃动的 目标多曝光序列图像如图2所示,成像系统晃动给图像之间带来的影响不超过五个像素。图 2中,图2(a)的曝光时间为Ι/lOs,图2(b)的曝光时间为l/5s,图2(c)的曝光时间为l/2s,图2 (d)的曝光时间为0.77s,图2(e)的曝光时间为2s,图2(f)的曝光时间为5s,图2(g)的曝光时 间为8s,图2(h)的曝光时间为13s。
[0111] (2)分解彩色图像。将获得的目标多曝光彩色图像按照R、G、B颜色通道进行分解, 获得各颜色通道的亮度信息。
[0112] (3)从序列图像中选取图像对获得亮度映射函数,选取图2(a)和图2(b)、图2(c)和 图2(d)、图2(e)和图2(f)、图2(g)和图2(h)组成四组图像对,分别获得四组图像对各颜色通 道相应的亮度传递函数曲线如图3所示。图3(a)、图3(b)、图3(c)分别为红色通道、绿色通 道、蓝色通道不同输入图像对的亮度映射函数。
[0113] (4)利用亮度传递函数进行相机响应函数标定,相机响应逆函数标定结果如图4所 不。
[0114] 利用均方根误差对存在动态问题影响条件下的相机响应函数标定结果进行定量 的评价,同时采用目前成熟的Debevec和Malik提出的相机响应标定算法对相同的输入图像 进行相机响应函数标定,对比结果如下表所示:
[0115]
[0116] 通过均方根误差可以观察到,针对存在动态问题影响到的图像序列进行相机响应 函数标定,基于亮度映射函数的相机响应函数标定方法仍然能够获得较为理想的标定结 果,而Debevec和Malik提出算法无法对存在动态问题影响的输入图像进行相机响应函数标 定。
[0117] (5)利用图2中多曝光输入图像序列结合构建的高动态范围图像权重融合函数进 行高动态范围图像融合,结果如图5所示。
[0118] 从实验结果可以看出,在同一帧结果图像中易于过曝光的照明灯处以及易于欠曝 光无法观察的办公桌下阴暗处均能在同一帧结果图像中得以显示,且克服了相机晃动带来 的影响,从而实现了高动态范围的图像融合。
[0119] 本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。
【主权项】
1. 一种用于克服动态问题影响的高动态范围图像融合方法,其特征在于包括如下步 骤: (1) 利用相机获得目标不同曝光量的曝光图像序列; (2) 根据目标的曝光图像序列,获取图像的亮度映射函数T(m),具体为: (21) 将曝光图像序列按曝光量从低到高进行排列,然后从第一帧图像开始每两帧图像 组成一个输入图像对; (22) 针对每一个输入图像对,获得该输入图像对中两帧图像各颜色通道的亮度信息的 累积直方图; (23) 利用输入图像对的两帧图像各颜色通道亮度信息的累积直方图,计算两帧图像各 颜色通道之间的亮度映射函数T(m),其中m为亮度值,HA(m)为输入图像对中的第一帧图像的某颜色通道亮度值为m的累积直 方图函数,为输入图像对中的第二帧图像的对应颜色通道亮度值为m的累积直方图 逆函数; (3) 利用亮度映射函数T(m)对相机响应函数进行标定,获得相机响应逆函数的近似表 达式g; (4) 结合相机响应逆函数的近似表达式g及曝光图像序列曝光时间信息△ t构建图像权 重融合函数I1,获得彩色图像各颜色通道的辐照度信息,从而获得反映拍摄场景信息的高 动态范围图像。2. 根据权利要求1所述的一种用于克服动态问题影响的高动态范围图像融合方法,其 特征在于:所述的图像权重融合函数I1S:式中i为图像中像素点的空间索引,j为曝光图像序列中的图像索引,M1伪输入曝光图 像序列中第j帧图像上像素点i的亮度信息,w为三角帽权值函数,△ tj为曝光图像序列中第 j帧图像的曝光时间,P为曝光图像序列中的图像总帧数。3. 根据权利要求1或2所述的一种用于克服动态问题影响的高动态范围图像融合方法, 其特征在于:所述步骤(3)中获得相机响应逆函数的近似表达式g的方法为: (31) 获取两帧不同曝光量图像对应关系如下: g(T (Ma))= kg (Ma) 其中Ma图像A上被估计的像素点的图像亮度信息,T是图像B和图像A之间的亮度传递函 数,k是两帧不同曝光量图像中图像A和图像B之间的曝光比; (32) 加入权重函数,将步骤(31)的结果修改为:其中Ca是第一帧图像A中亮度值为Ma的像素个数,Cb是第二帧图像B中亮度值为Mb的像素 个数; (33) 利用多项式对相机响应逆函数进行建模,将不同输入图像对获得的亮度传递函数 以及相机响应逆函数多项式代入步骤(32)的表达式中获得一组超定方程; (34) 利用总体最小二乘法求解所述的超定方程,获得相机响应逆函数的模型系数,由 此得到相机响应逆函数的近似表达式g。
【文档编号】H04N5/235GK105933617SQ201610335204
【公开日】2016年9月7日
【申请日】2016年5月19日
【发明人】都琳, 孙华燕, 张廷华, 王帅, 徐韬祜, 李迎春, 梁丹华, 高宇轩
【申请人】中国人民解放军装备学院
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1