倪育才 不确定度之1如何建立数学模型.docx

上传人:b****6 文档编号:4615944 上传时间:2022-12-07 格式:DOCX 页数:21 大小:230.41KB
下载 相关 举报
倪育才 不确定度之1如何建立数学模型.docx_第1页
第1页 / 共21页
倪育才 不确定度之1如何建立数学模型.docx_第2页
第2页 / 共21页
倪育才 不确定度之1如何建立数学模型.docx_第3页
第3页 / 共21页
倪育才 不确定度之1如何建立数学模型.docx_第4页
第4页 / 共21页
倪育才 不确定度之1如何建立数学模型.docx_第5页
第5页 / 共21页
点击查看更多>>
下载资源
资源描述

倪育才 不确定度之1如何建立数学模型.docx

《倪育才 不确定度之1如何建立数学模型.docx》由会员分享,可在线阅读,更多相关《倪育才 不确定度之1如何建立数学模型.docx(21页珍藏版)》请在冰豆网上搜索。

倪育才 不确定度之1如何建立数学模型.docx

倪育才不确定度之1如何建立数学模型

如何建立数学模型

讲授人:

中国计量科学研究院研究员倪育才

     应《中国计量》编辑部的邀请,根据近几年来在各种“测量不确定度评定”学习班上讲课的经验和对初学者经常感到困惑的问题或不容易理解的问题,现写若干短篇并在《中国计量》上连载。

为使这些短篇相互独立,每一篇讲述一个问题,合起来看又有一定的系统性。

经与编辑部商讨,初步确定下述几个问题:

   1  如何建立数学模型

      2  极差法和贝塞尔法之间的比较

      3  被测量Y可能值分布的判定

      4  包含因子k的选择

      5  测量不确定度评定在不同应用中的差别

      6  测量误差的基本概念

      7  测量不确定度的基本概念

      8  测量误差和测量不确定度的差别

    在测量不确定度评定中,建立数学模型也称为测量模型化,目的是要建立满足测量不确定度评定所要求的数学模型,即建立被测量Y和所有各影响量X间的函数关系,其一般形式可写为:

    Y=f(X1,X2,…,Xn)

    可以说,建立数学模型是进行测量不确定度评定最关键的第一步,也是许多初学者在进行测量不确定度评定时遇到的第一个困难。

    《测量不确定度表示指南》(GUM)在摘要介绍测量不确定度评定步骤时,首先就提到要建立数学模型,并说:

“Thefunctionfshouldcontaineveryquantity,includingallcorrectionsandcorrectionfactors,thatcancontributeasignificantcomponentofuncertaintytotheresultofmeasurement.”。

其意是数学模型f中应包含所有对测量结果的不确定度有影响的修正值和修正因子。

也就是说,数学模型中应包含所有应该考虑的影响量,而每一个影响量将对测量结果贡献一个值得考虑的不确定度分量。

因此一个好的数学模型,其中所包含的影响量和此后不确定度评定中所考虑的每一个不确定度分量应该是一一对应的。

这样建立起来的数学模型,既能用来计算测量结果,又能用来全面地评定测量结果的不确定度。

    要找出每一个影响量与被测量之间的函数关系,往往是很困难的,有时简直不可能得到两者关系的解析表达式。

于是许多初学者往往将测量中用来获得被测量的计算公式作为数学模型而列出。

例如在各种测量中,最经常采用的方法之一是比较测量。

将被测量值y和参考标准所提供的标准量值s相比较,通过测量两者之差Δ可以计算出被测量y。

于是在已经发表的各种测量不确定度评定的文章中,经常见到将y=x+Δ作为数学模型的情况。

但在进行不确定度评定时,则又往往脱离数学模型而重新考虑各个不确定度分量。

这样的数学模型对测量不确定度评定实际上毫无帮助。

    在某些特殊情况下(例如某些检测项目)将计算公式作为数学模型可能是允许的,但一般说来不要把数学模型简单地理解为就是计算测量结果的公式,也不要理解为就是测量的基本原理公式。

两者之间经常是有区别的。

    从原则上说,似乎所有对测量结果有影响的输入量都应该在计算公式中出现,但实际情况却不然。

有些输入量虽然对测量结果有影响,但由于信息量的缺乏,在具体测量时无法定量地计算它们对测量结果的影响。

也有些输入量由于对测量结果的影响很小而被忽略,故在测量结果的计算公式中也不出现,但它们对测量结果的不确定度的影响却可能是必须考虑的。

因此如果仅从计算公式出发来进行不确定度评定,则上述这些不确定度分量就可能被遗漏。

当然,在某些特殊情况下如果所有其他不确定度贡献因素的影响都可以忽略不计时,数学模型也可能与计算公式相同。

    对于不同的被测量和不同的测量方法,数学模型的具体形式可能差别很大,但实际上都可以用一种比较系统的方式来给出数学模型,或者说可以给出数学模型的通式。

    根据测量误差的定义:

误差=测量结果-真值。

同时误差又可以分为随机误差和系统误差两类,且三者之间的关系为:

误差=系统误差+随机误差。

于是可以得到:

    真值=测量结果-误差

          =测量结果-系统误差-随机误差

    由于修正值等于负的误差,于是上面的关系式就成为:

    真值=测量结果-系统误差-随机误差

          =测量结果+系统误差的修正值+随机误差的修正值

    实际上,真值就是想得到的被测量的测量结果,于是上式可写成

    被测量=测量结果+系统误差的修正值+随机误差的修正值

   例1:

对于常见的量块比较测量,若ls为标准量块的长度,Δl为测得的两量块的长度差,于是被测量块长度lx的计算公式为:

    lx=ls+Δl

    由于测量时量块的温度通常会偏离标准参考温度20℃,考虑到温度和线膨胀系数对测量结果的影响,计算公式成为:

    lx=ls+Δl+lsδαθx+lsαsδθ

    式中α和θ分别表示线膨胀系数和对标准参考温度20℃的偏差;脚标“s”、“x”分别表示标准量块和被测量块;以及δθ=θs-θx和δα=αs-αx。

    考虑到量块测量点可能偏离量块测量面中心点对测量结果的影响,数学模型成为:

    lx=ls+Δl+lsδαθx+lsαsδθ+δl

    将此数学模型和上面给出的通式相比较就可以发现,等式右边的第一、二项ls+Δl即是由测量得到的未修正的测量结果。

等式右边的第三、四项lsδαθx+lsαsδθ是对由温度偏差所引入的系统误差的修正值,在本例中这两项的数值十分小而可以忽略,但它们对测量结果不确定度的影响是必须考虑的。

等式右边的最后一项δl,是表示由于测量点可能偏离量块中心对测量结果的影响。

测量点的偏离对测量结果引入随机误差,因此最后一项实际上是对该随机误差的修正值。

由下图可见两者之间的对应关系。

   例2:

砝码校准,将被测砝码的质量与具有相同标称值的标准砝码相比较。

若被校准砝码和标准砝码的折算质量分别为mx和ms,测得两者的质量差为Δm,于是被校准砝码折算质量mx的计算公式为:

    mx=ms+Δm

    考虑到标准砝码的质量自最近一次校准以来可能产生的漂移Δmd,质量比较仪的偏心度和磁效应的影响Δmc,以及空气浮力对测量结果的影响δB后,其数学模型成为:

    mx=ms+Δm+δmd+δmc+δB

    模型中等式右边的第一、二项为未修正的测量结果。

该测量不存在值得考虑的系统误差,也就是说,在数学模型中不存在对系统误差的修正值。

等式右边的第三、四、五项为对三项随机误差分量的修正量。

与数学模型通式之间的对应关系为:

   在建立数学模型时,未修正的测量结果和系统误差的修正值通常都能比较容易地得到解析形式的数学表达式。

惟有随机误差的修正值无法得到其解析形式的表达式。

因此只能在数学模型中简单地加上一项,表示对随机误差的修正值。

根据随机误差的定义,无限多次测量结果的随机误差的平均值等于零,因此这些项的数学期望为零。

也就是说,增加这些修正值后不会对被测量的数值有影响。

需要知道的是这些修正值的可能取值范围,通常可以由测量者的经验或辅助的实验测量得到。

再由假定的概率分布,可以通过B类评定估算出它们的标准不确定度。

    有些测量,其计算公式中可能仅包含各影响量的积和商,即被测量可以用下述函数形式表示:

   

   式中的系数c并非灵敏系数,而是比例常数,且指数pi可以为正数或负数。

在这种情况下,需要增加的不是修正值,而是相乘的修正因子。

此时,数学模型的通式可以表示为:

被测量等于未修正测量结果的计算公式乘以由于系统误差引入的修正因子(它们的数学期望值不等于1),再乘以由于随机误差引入的修正因子(它们的数学期望值等于1)。

    有些领域,例如化学分析领域,经常出现这种类型的数学模型。

   例3:

在用原子吸收光谱法测定陶瓷容器中镉的溶出量的实例中,被测量为被醋酸溶液浸泡的容器单位表面积镉的溶出量r,它可以表示为:

    

    式中:

ρ0——醋酸浸取液中镉的质量浓度;VL——醋酸浸取液体积;aV——被醋酸溶液浸泡的容器表面积。

    考虑到还有三项随机误差在上述公式中未反映出来,它们分别是浸泡温度、浸泡时间和醋酸的体积分数对测量结果的影响,于是最后采用的数学模型成为:

    

    在该数学模型中,

是未修正的测量结果,ftemp、ftime和facid分别是相对于三项随机误差的修正因子,它们的数学期望均等于1。

在本例中不存在值得考虑的系统误差。

    由此可见,写出符合要求的数学模型并不难,关键还是要找到所有能影响测量结果的误差来源。

一般先根据测量的最基本原理导出被测量的基本计算公式,然后考察该计算公式是否已经对所有的系统误差进行了修正,否则就补充加入其余未考虑的系统误差分量的修正值(或乘以修正因子),最后再加上对所有随机误差分量的修正值(或乘以修正因子)。

只要对测量工作有一定程度的了解,写出计算公式和系统误差修正值的函数形式对大部分测量人员并不困难,因此要做的仅是简单地将所有需要考虑的随机误差的修正值(或修正因子)补充进入数学模型。

    必须注意,即使对于相同的被测量和相同的测量方法,数学模型也不是一成不变的。

随着所选择的影响量的不同,对测量不确定度评定所要求的严密程度的不同,其数学模型也可能会有所不同。

    此外,对于测量仪器和量具的常规检定或校准来说,还必须注意两者在数学模型上可能存在的微小差别。

当被测对象是测量仪器时,由于仪器本身一般不提供标准量值,其量值需要用其他测量标准进行标定。

故在进行测量不确定度评定时,被测量应该是测量仪器的示值误差Ex,因此其数学模型需写成示值误差的形式,即“Ex=……”。

当被测对象是实物量具时,由于实物量具本身能提供一个标准量值,故在进行测量不确定度评定时,被测量既可以是其相对于标称值的偏差(相当于示值误差),也可以是它所提供的量值。

也就是说,其数学模型既可以写成“y=……”的形式,也可以写成“Ex=……”。

由于两者之间仅相差一个标称值,而标称值是一个规定值而不存在不确定度,因此两种数学模型在不确定度评定时毫无差别。

极差法和贝塞尔法之间的比较

讲授人:

中国计量科学研究院研究员倪育才

  标准不确定度的A类评定定义为:

“用对观测列进行统计分析的方法,来评定标准不确定度”。

国家计量技术规范JJF1059-1999《测量不确定度评定与表示》中介绍了两种A类评定的方法,贝塞尔法和极差法。

    1.贝塞尔法

    当在重复性或复现性条件下,对被测量X进行n次独立观测。

若得到的测量结果分别为x1,x2,……,xn,n次测量的平均值为

于是用贝塞尔公式可以求出单次测量结果xi的实验方差s2(xi)和实验标准差s(xi)。

    

    2.极差法

      当在重复性或复现性条件下,对被测量X进行n次独立观测。

若n个测量结果中最大值和最小值之差为R(称为极差),在可以估计X接近正态分布的条件下,单次测量结果的实验标准差s(xiv)可近似地表示为:

    s(xi)=R/C=u(xi)

      式中系数C为极差系数。

极差系数之值与测量次数n的大小有关。

表1给出极差法的极差系数和自由度与测量次数的关系。

    

  

     既然随机变量X的标准偏差可以用两种方法得到,就不可避免地会提出两种方法孰优孰劣的问题。

无疑,极差法具有计算简单的优点。

但在计算机应用已经十分普及的今天,用贝塞尔公式计算也已变得相当容易。

因此关键问题还在于用何种方法估算得到的不确定度更为准确。

    表面上看来,用贝塞尔公式进行计算时使用了全部n个测量结果,而极差法只用了一个极大值和一个极小值,其余数据均弃之不用,因此用贝塞尔法得到的实验标准差应该比极差法更为可靠。

比较两种方法的自由度也可以看出,极差法的自由度比贝塞尔法小(贝塞尔法的自由度为n-1,而极差法的自由度

于是可以得到同样的结论,贝塞尔法比极差法更为可靠。

    但实际上问题并没有这么简单。

根据定义,用标准偏差表示的不确定度称为标准不确定度。

因此从理论上说,应该计算的是标准偏差σ,而不是实验标准差s。

但标准偏差是一个总体参数,也就是说,要进行无限多次测量才能得到。

在实际工作中只能用样本参数来代替总体参数,即用实验标准差s来作为标准偏差σ的估计量。

理论上可以证明,实验标准差s并不是标准偏差σ的无偏估计量。

这就是说,当用实验标准差s来代替标准偏差σ时,除了实验标准差s本身是一个随机变量外,它的数学期望值(即无限多次测量结果的平均值)相对于标准偏差σ还有一个与测量次数有关的系统性偏差。

测量次数越少,其系统性偏差就越大。

因此可以对贝塞尔公式作一无偏差的修正。

经过无偏差修正后的贝塞尔公式为:

      

    上式中修正因子Mn的数值见表2。

由表2可知,当测量次数n≤6时,随着测量次数减少,偏离系数Mn将明显加速偏离1。

    

 

   也可以分别计算出用贝塞尔公式和极差法得到的实验标准差的相对标准不确定度,其计算结果见表3。

由表3可以看出,当测量次数n=10时,两种方法得到的实验标准差准确程度几乎相同。

当n>10时,贝塞尔法优于极差法;当n<10时,极差法优于贝塞尔法。

至于修正的贝塞尔公式,相比而言虽然最为准确,但因比较麻烦实际上很少使用。

这就是为什么国家计量技术规范JJF1059-1999中在给出极差系数及自由度表后指出“一般在测量次数较小时采用该法”,以及国家计量技术法规统一宣贯教材《测量不确定度评定与表示指南》中同时还指出“测量次数以4~9次为宜”。

      

上面的分析,仅是针对实验标准差而言的。

在大部分的测量不确定度评定中,测量不确定度A类评定仅是其中的一个或几个分量。

他们还将与其他B类评定的分量合成,才能得到合成标准不确定度。

合成的方法是方差相加。

虽然实验标准差s并不是标准偏差的无偏估计量,但却可以证明实验方差s2是总体方差σ2的无偏估计量。

因此,若A类评定需要和其他B类分量合成,且A类评定分量不占优势时,则无论测量次数的多少,贝塞尔法将优于极差法。

    因此笔者认为结论应该是:

    

(1)当A类评定不确定度分量不是合成标准不确定度中惟一占优势的分量时,则无论测量次数多少,贝塞尔法优于极差法。

      

(2)当A类评定不确定度分量是合成标准不确定度中惟一占优势的分量时,则两种方法的优劣与测量次数有关。

当测量次数n<10时,极差法优于贝塞尔法;当测量次数n≥10时,贝塞尔法优于极差法。

测量结果与不确定度表示

讲授人:

李慎安

  JJF1059第8.13节指出输入量和输出量的估计值,应修约到与它们的不确定度的位数一致。

这里所谓的位数实指其末位所到达的位数。

例如,当测量结果及其不确定度以相同的计量单位给出时,其末位应对齐。

也就是说不能达不到,也不能多出。

其中更需注意的是所报告的测量结果(输出量的最佳估计值),应与所报告的扩展不确定度U或Up的末位对齐。

多数情况下是:

确定了扩展不确定度取几位(一或两位)之后,按这一修约间隔来修约所报告的测量结果。

但有时也会碰到,特别是通过数字显示式仪器的一次测量结果作为被测量的最终结果时,评定出的扩展不确定度的末位已小于所显示的末位。

这时,对测量结果是否能采用补零的方式使其末位对齐?

专家们对不同意见进行了讨论,例如:

通过数字式电压表一次测量的结果为220.043V,其扩展不确定度U=2.5mV(k=2),U修约成两位,末位达到0.1mV,但测量结果只到1mV,专家们认为这时的测量结果应报告成:

220.0430V。

写成V=(220.0430±0.0025)V,其末位是对齐的。

应该认为,表明测量结果可靠程度的不是所给出的结果本身而是其不确定度。

那种认为物理实验结果只能保留一位不可靠的值(只有末位不可靠而不能有两位是不可靠的)的观点和做法,与当今不确定度的表述并不一致。

现在认为不确定度可以有两位有效数,从而测量结果的末两位均为可疑值了。

    关于所报告的扩展不确定度(U,Up和Urel,Uprel)应采取何种规则进行修约,在JJF1059第8.13节给出两种方法均可以用,其一为“只进不舍”,其二为通用的修约规则,即大于半个修约间隔则进,小于半个修约间隔则舍,正好等于半个修约间隔则看前面一位是奇数还是偶数而定。

根据第一种方法,如果对U=0.1112修约成为一位有效数,按只进不舍,就成为U=0.2,比修约前增大了几乎一倍,虽不违反规则,但显然并不可取。

如果U=0.3112,也只取一位有效数而给成为U=0.4,比修约前也大了1/4左右,似亦不可取。

专家们推荐采用:

当第一个有效数为1和2时,取两位有效数为好,至于3以上,既可取一位也可取两位,对于一般测量,可均只取一位。

至于是按上述两种修约方法中的哪一种,评定人员可自行选用。

上述的这种建议,在JJF1059以及GUM中都未提及,只是在某些国家的标准中提到,例如DIN,不无道理,未必不可以参照使用。

    现在在一些检定证书或是校准证书上,给出了测量结果(校准结果、某些检定点或校准点的示值误差或修正值)。

对于校准(自愿行为),给出校准值及其不确定度,是符合JJF1059中8.2节要求“证书上的校准结果或修正值应给出测量不确定度。

”但是在检定中,例如:

对压力表、千分尺、台案秤等类衡器,按检定规程,其证书上是不给出测量结果的,现在也要求给出检定结果,有时甚至也给出其不确定度。

从测量仪器的使用上来说,这些内容不起任何作用,因不能按测量结果修正使用。

惟一的作用是让使用者知道这些仪器距离不合格还有多远。

专家们认为,究竟在证书上如何给出和给出什么,应按有关规程处理,至于自愿的校准要求,则可按用户需要。

    关于测量仪器特性评定问题,目前仍按JJG1027-1991技术规范中的有关规定处理。

计量司官员在会上表示,用于代替该技术规范这部分的内容的新的技术规范现已审定通过,处于报批之中,预计今年内可发布。

其中规定了测量仪器特性评定的基本原则、通用方法、准确度等、级、响应特性、灵敏度、鉴别力、稳定性、漂移、响应时间等性能的评定以及有关不确定度问题。

关于测量仪器重复性的评定,该规范给出了基本方法,即按重复性条件下通过重复观测,采用贝塞尔公式计算出单次结果的实验标准差s。

s的相对标准不确定度:

    

  

    式中:

n——重复观测次数。

    对于只有一个被测量来说,上式也就是:

    

  

    式中:

ν——标准差s的自由度。

    该标准还给出了最大残差法用于测量仪器重复性的评定。

上述内容可作为JJF1059的补充。

    按该新规范,可以采用MPEV作为测量仪器的最大允许误差绝对值的符号。

包含因子K的选择

讲授人:

中国计量科学研究院研究员倪育才

  当得到合成标准不确定度后,获得扩展不确定度Up的前提是确定包含因子k的数值。

包含因子的确定方法取决于被测量的分布,因此当被测量Y的分布不同时,应采用不同的方法来确定包含因子。

一、当无法判断被测量Y的分布时

   当无法判断被测量Y的分布时,不可能根据分布来确定包含因子k。

由于大部分测量均规定要给出扩展不确定度U而不是Up,因此只能假定取k=2或3,绝大部分情况下均取k=2(JJF1059-1999规定,当包含因子取其他值时,应说明其来源)。

于是扩展不确定度成为:

    U=2uc

    由于不知道被测量的分布,故无法建立置信概率p和包含因子k之间的关系。

此时的k值是假设的,而不是由置信概率p导出的,也就是说,无法知道此时所对应的置信概率。

二、当被测量Y接近于某种非正态分布

   当被测量接近于某种已知的非正态分布时,例如矩形分布、三角分布、梯形分布等,则绝不应该按上面的方法直接取k=2或3,也不能按正态分布的方法。

根据计算得到的有效自由度νeff,并由t分布表得到kp。

此时应根据已经确定的被测量Y的分布,由其概率密度函数具体计算出包含因子k。

    

(1)当可以判定被测量Y接近于矩形分布时,由其概率密度函数可以计算得到包含因子k与置信概率p之间的关系为:

    

    当p=0.95时,k95=1.65;

    当p=0.99时,k99=1.71。

    

(2)当可以判定被测量Y接近于三角分布时,通过类似的计算可以得到包含因子k与置信概率p之间的关系为:

    

    当p=0.95时,k95=1.90;

    当p=0.99时,k99=2.20。

    (3)当可以判定被测量接近于梯形分布时,通过计算可以得到包含因子k与置信概率p之间的关系为:

    

    

    式中β为梯形的角参数,即梯形的上底和下底之比。

    表1给出当所要求的置信概率分别为95%和99%时,由上式计算得到的不同β值梯形分布的包含因子k。

三、被测量Y接近于正态分布时

   在被测量接近正态分布的情况下,不能直接取正态分布所对应的k值。

由于标准不确定度定义为以标准偏差表示的不确定度。

而标准偏差是一总体参数,只有通过无限多次测量才能够得到,正态分布也是对应于无穷多次测量的总体分布。

也就是说,只有当用总体标准偏差σ来作为标准不确定度时,才能采用正态分布的k值。

但由于在实际测量中不可能进行无限多次测量,只能用有限次测量的实验标准差s作为σ的估计值,并且这一估计必然会引入误差。

由于该误差的存在,如果仍采用正态分布的k值,将达不到所要求的置信概率。

反过来说,为了得到对应于所规定置信概率的扩展不确定度,必须适当增大k值。

并且随着测量次数的减少,用实验标准差代替标准偏差可能引入的误差将越来越大,包含因子k的值也必将随之增加。

因此,这时的包含因子k将是一个与测量次数有关的变量。

    在数学上,这相当于总体分布满足正态分布时,其样本分布满足t分布。

t分布是表征正态分布总体中所取子样的分布。

不同的子样大小,对应于不同的t分布,其包含因子k也将不同。

因此当被测量Y接近于正态分布时,仅仅根据所要求的置信概率还不足以得到包含因子k,还必须再知道一个与所取样本大小有关的参数,这个参数就称为“自由度”,一般用希腊字母ν表示。

对于不同的自由度,包含因子kp=tp(ν)的数值可以由所规定的置信概率p和估计得到的有效自由度νeff通过查表得到。

    JJF1059-1999规定,当可以判断被测量Y接近于正态分布时,可以采用以下方法得到扩展不确定度。

    通过计算被测量Y的有效自由度νeff,并根据有效自由度和所要求的置信概率p由t分布临界值表得到包含因子kp=tp(νeff),于是扩展不确定度Up等于合成标准不确定度uc和包含因子kp=tp(νeff)的乘积,即:

    Up=kp·uc

    此时,由于包含因子k的数值由置信概率通过t分布表得到,因此得到的扩展不确定度必须以Up表示,表明所给的是对应于置信概率为p的扩展不确定度。

    有时为简单起见,也可以不必考虑比较麻烦的Y的分布情况。

此时可采用k=2或3,即:

    U=kuc  (k=2或3)

    此时,由于包含因子k的数值是假定的,因此得到的扩展不确定度只能以U表示,表明无法知道与所给扩展不确定度对应的置信概率。

    这种情况就是无法判断被测量Y的分布差所采取的做法,在已知正态分布的情况下,只要有效自由度不太小,当k分别取2或3时,他们大体上对应于95%或99%的置信概率。

    但如可估计为正态分布,而有效自由度较小时,在k取2(或3)时,则所给扩展不确定度对应的置信概率可能会与95%(或99%)相去甚远。

因此笔者建议仅在确保有效自由度不太小的情况下(例如不小于15)采用该法,除非该领域统一规定直接

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 工程科技 > 冶金矿山地质

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1