季丰民，王仙同，万　星，沈建眯

(1. 宁波工程学院 材料学院，宁波 315016； 2. 宁波工程学院 理学院，宁波 315211)



比辐射率反演的一个新方法

季丰民1，王仙同2，万星2，沈建眯2

(1. 宁波工程学院 材料学院，宁波 315016； 2. 宁波工程学院 理学院，宁波 315211)

物体的比辐射率是反映物体热辐射性质的一个重要参数，比辐射率反演不仅提供了一种获得物体比辐射率的便捷途径，并且涉及许多有趣的理论问题.目前该反问题只有一种基于Hermite函数展开的普适函数方法.考虑到这类反问题固有的不适定性问题，发展另外一种独立的方法非常有必要.本文基于Laguerre函数展开为比辐射率反问题发展了一种新的方法.对一些理论谱的数值计算表明，新方法不仅拥有老方法的所有优点，并且跟老方法相比拥有更高的计算精度，能够计算结构更加复杂的理论谱.本文最后通过比较和分析揭示了新方法优越性的根本原因.

物理学中的反问题； 比辐射率； 渐近行为控制条件； 普适函数方法

反问题研究已经成为现代物理学的一个新的生长点，例如在统计物理学中出现的一类反问题开辟了理论物理研究的广阔领域，它们包括比热-声子谱反问题[1-6]、黑体辐射反问题[7-11]、比辐射率反问题[12]、透射率反问题[13]等.它们都可以表示成Planck型积分方程，即有关化学势为零的玻色系统(如声子、光子等)的物理性质的积分方程.本文着重探讨比辐射率反演问题，这个反问题为人们提供了一种获得物体比辐射率的便捷途径.

就物体辐射电磁波的能力来讲，一般的物体总是比黑体辐射要弱一些，我们称之为“灰体辐射”.灰体的辐射能力和黑体相比之后得到的比值定义为灰体的比辐射率g(ν)，其中ν为辐射电磁波的频率.由此灰体的比辐射率总有0≤g(ν)≤1.根据Planck黑体辐射理论，物体单位表面积的辐射功率谱为

(1)

其中：h，c，kB分别为Planck常数，光速和Boltzmann常数，T是热力学温度.假如我们知道了灰体的比辐射率，那么物体单位表面积的辐射总功率J(T)就是辐射功率谱在整个频率空间的积分:

(2)

基于反遥测领域的研究需要，人们提出了该辐射问题的反问题[12]： 能否通过测量灰体单位表面积的辐射总功率，求解上述积分方程而得到该灰体的比辐射率？这个反演问题被称为比辐射率反问题.大家知道，反遥测技术是一项有重要实用价值的工程技术，而比辐射率的研究是反遥测技术中的重要内容.在反遥测技术中，如何在特定的背景里隐藏自己、躲避比如红外探测等的侦察至关重要，一个很重要的途径就是降低己方物体的比辐射率.考虑到一些物体的比辐射率还可能跟温度有关，因此比辐射率应该写成g(ν，T).根据测得的灰体辐射总功率求解温度有关的比辐射率称为广义比辐射率反问题，由明灯明等[14]首次提出并做了一些近似计算.本文仅限于讨论狭义比辐射率的反演问题.

比辐射率可以在实验上利用谱分析的方法获得，但是这样的测量需要相当复杂的仪器，对于不同的频率测量区间，适合的探测器就有所不同，要制作适应于不同频率的探测器也是比较困难的.因此通过反演的方法获取比辐射率数据就成为一个比较实际可行的途径.

人们曾经基于Hermite函数展开，为比辐射率反问题发展了一种普适函数方法(Universal Function Set method based on Hermite functions， UFS-H)[12]，并利用这种方法计算了一个单峰的理论谱，计算结果跟理论符合得很好.但是，UFS-H方法在双峰理论谱上的表现不尽如人意，即使计算至很高的阶数，计算结果依然跟理论值存在可观的偏差.考虑到比辐射率反问题涉及第一类Fredholm型积分方程，在求解过程中不可避免的会遇到所谓的不适定性问题[15-18]： 微小的输入差别可能导致巨大的输出不同.如果要计算实际体系的比辐射率反演，UFS-H方法恐怕会更加难以胜任.其次，对于比辐射率反问题而言，两种独立方法之间的互相检验也是非常有必要的.本文主要目的是基于Laguerre函数集在比辐射率反演问题当中建立一种新的方法.

1　严格解公式简要回顾

在提出比辐射率反问题之后，人们参照相关反问题[2，8]的解决思路，基于傅里叶变换推导得到了一个封闭(积分)形式的严格解公式[12]：

(3)

(4)

另外，参数s>1的引入是为了消除傅里叶变换的发散问题，并且避开数学上著名的黎曼猜测[19].值得说明的是，在严格解公式中输入函数是辐射总功率J(T)，输出的是物体的比辐射率g(ν).

可以证明[12]，如果辐射总功率满足如下条件:

(5)

那么比辐射率反演问题的解存在并且唯一.上述条件中的记号o[…]表示高阶无穷小，这个条件被称为渐近行为控制(Asymptotic Behavior Control， ABC)条件.

2　一种新的反演方法

在发展一种新方法的过程中，ABC条件起着关键的引导作用.它是辐射总功率必须满足的条件，但是由于任何实验测量数据必然存在误差，而且测量的温度区间不可能无穷大，导致ABC条件不可能被实验测量数据所直接满足.所以新方法必须具有自动筛选功能，把符合ABC条件的有效成分予以保留，把不符合ABC条件的测量噪音剔除在外.这样就对展开函数有了极其严格的要求，不仅要求它们的傅里叶变换每一项都严格可解，并且它们的渐近行为都要求满足ABC条件.

我们曾经参照以前的UFS-H方法，试图直接用扩展了定义域之后的Laguerre函数展开Q0(x)，以此发展出一种新的方法.但是我们很快发现，与Hermite函数的情况不同，Laguerre函数的傅里叶变换并不满足ABC条件，我们必须另辟蹊径！

(6)

这里的ξ是一个实参数，展开系数Cm是复数，Laguerre函数Lm(ξk)定义为

(7)

(8)

并且Laguerre函数的渐近行为有如下表达式：

(9)

(10)

(11)

则比辐射率反问题存在唯一解.对比前面Laguerre函数的渐近行为，我们发现只要实参数ξ取值大于π，Laguerre函数的每一项都能满足该反问题的ABC条件.

(12)

其中Cm是展开系数，并且

(13)

不难发现Gm(ν)是一个纯数学的、与具体体系无关的函数集，这一点和以前的UFS-H方法相似，我们也称它为新方法中的普适函数集，它可以事先用计算机计算完成.这种新方法我们称之为基于Laguerre函数的普适函数方法(Universal Function Set method based on Laguerre functions， UFS-L).

根据Laguerre函数的正交性以及傅里叶变换的定义，我们有

(14)

因为Laguerre函数Lm(ξk)是实函数，上式中的函数Dm(x)可以表示为

(15)

接下来我们必须解析完成Dm(x)函数里面的积分，找到一个尽可能简单的表达式.很幸运我们得到这个新函数的简单形式

(16)

在推导过程的最后一步，我们应用了二项式定理完成了求和过程.

到现在为止，一种基于Laguerre函数展开的新的普适函数方法(UFS-L)已经建立起来了.根据测量得到的灰体辐射总功率实验数据，通过以下步骤即可反演计算灰体的比辐射率： 首先，将辐射总功率离散的实验数据拟合成一个连续光滑的函数，根据(4)式转换成Q0(x)；其次，根据(14)式计算展开系数；再次，根据(13)式计算普适函数集；最后，将展开系数和普适函数代入(12)式通过简单的求和得到比辐射率.

3　新方法的检验

大家知道，在统计物理中严格解是比较稀少的.在比辐射率反问题中，我们很幸运地找到如下一组严格解，它们能够很好地检验反演方法的可行性以及它们的优劣:

(17)

其中:A和μ均为实参数，η为频率的标度参数，ζ(z，q)为广义黎曼函数.

上述严格解(理论谱)的一个单模曾经用UFS-H方法做过计算，计算结果和理论值符合很好[12].作为检验，我们用新的UFS-L方法计算J1(T)，相关参数和常数的取值为A=3.0×1013，μ=2，η=2.86×10-15Hz-1，s=4.0，T0=2.5×104K，ξ=9.0.根据新的UFS-L方法，要得到比辐射率，只需要计算普适函数和展开系数.其中普适函数与具体的体系无关，因此可以事先运用计算机进行计算.计算好展开系数之后，把普适函数和展开系数代入公式(12)即得该单模对应的比辐射率g1(ν).图1显示了UFS-L方法的计算结果以及与理论值的比较，从中可见新方法的确是可行的，并且相当成功，展开级数从第10阶左右开始稳定收敛.

计算单峰的理论谱无论对于之前的UFS-H方法还是现在的UFS-L方法都是相对比较容易的，计算结构更加复杂的理论谱却对这些方法构成相当大的挑战.

假设有集合an={0.66，-12，98，-427，1120，-1820，1848，-1142，394，-58}，则以下比辐射率由两个峰构成:

(18)

(19)

以上的gn(ν)和Jn(T)定义参照(17)式.

这个双峰谱的计算比前面的单峰谱要困难一些，展开级数的收敛要比单峰谱慢一些，不过用UFS-L方法计算至36阶已经和理论符合得很好.为了和以往的UFS-H方法做比较，我们也用UFS-H方法计算了这个双峰谱，结果发现UFS-H方法在这个问题上收敛更慢，一直算到66阶(这几乎是PC电脑机计算普适函数的极限阶数)，UFS-H方法计算的结果和理论值依然存在微小的偏差.两种方法的计算结果以及理论值如图2所示.相关参数和常数的取值为：B=1.0×1013，μ=5，η=1.0×10-12Hz-1，s=3.5，T0=300K，ξ=9.0.通过图2，可以看出新方法在计算理论谱的时候，具有比以往方法更高的计算精度，因而它能够精确计算具有更加复杂结构的理论谱.通过两种方法的比对，我们发现关键在于UFS-L方法中的函数集Dm(x)，它比以往方法中的Hermite函数集更加简单.这个新的函数集只包含一些初等函数，甚至它们具有一个共同的包络线(如图3所示，参数取值ξ=10.0).

(20)

4　讨论与总结

比辐射率反演是统计物理学中出现的具有工程应用背景的反问题.本文基于数学上著名的Laguerre函数集，为比辐射率反问题发展了一种新的反演方法： UFS-L方法.新方法是继现有的UFS-H方法之后的又一种普适函数方法.其中的普适函数都是与实际体系无关的纯数学表达式，设定好几个参数之后即可以用计算机完成计算.与实际体系有关的只是展开系数，获得展开系数之后通过简单的求和过程即可获得比辐射率.由于比辐射率反问题固有的不适定性问题的存在，在处理实际问题时，反演结果往往是很难预估的，两种独立方法之间的相互验证是很有必要的.

对单峰理论谱和双峰理论谱的数值计算表明，新方法不仅拥有老方法的所有优点，并且跟老方法相比拥有更高的计算精度.究其原因，是因为新方法中导出的Dm(x)函数属于初等函数，比老方法中的特殊函数Hermite函数更加简单，计算更加容易.因此，这个新函数也可运用在统计物理学中出现的其他反问题(比如比热-声子谱反问题、黑体辐射反问题、透射率反问题等)中.

由于这种新方法中每一个展开项都满足ABC条件，所以它能自动筛选出辐射总功率测量数据中的有效信息成分，同时剔除由于测量误差以及有限测量区间引起的信息噪音，于是为真实体系的比辐射率反演铺平了道路.今后如果有更加精密和完善的辐射总功率实验测量数据，运用新方法有望实现首个实际体系的比辐射率反演.

[1]CHEN N X. Modified Möbius inverse formula and its application in physics [J].PhysRevLett， 1990，64(11)： 1193-1195.

[2]DAI X X， XU X W， DAI J Q. On a specific heat-phonon spectrum inversion problem. exact solution， unique existence theorem and riemann hypothesis [J].PhysLettA， 1990，147(8/9)： 445-449.

[3]DAI X X， WEN T， MA G C，etal. A concrete realization of specific heat-phonon spectrum inversion for YBCO [J].PhyLettA， 1999，264(1)： 68-73.

[4]WEN T， MA G C， DAI X X，etal. Evenson. Phonon spectrum of YBCO obtained by specific heat inversion method for real data [J].JPhys：CondensMatter， 2003，15(2)： 225-238.

[5]MING D M， WEN T， DAI J X，etal. A unified solution of the specific heat-phonon spectrum inversion problem [J].EurophysLett， 2003，61(6)： 723.

[6]JI F M， DAI X X， STEVENS R，etal. Thermodynamic functions of ZrW2O8from its heat capacity [J].SciChina-PhysMechAstron， 2012，55(4)： 563-567.

[7]BOJARSKI N N. Inverse black body radiation [J].IEEETransAntennasPropagation， 1982，30(4)： 778-780.

[8]DAI X X， DAI J Q. On unique existence theorem and exact solution formula of the inverse black-body radiation problem [J].IEEETransAntennasPropagation， 1992，40(3)： 257-260.

[9]YE J P， JI F M， DAI X X，etal. The black-body radiation inversion problem， its instability and a new universal function set method [J].PhysLettA， 2006，348(3-6)： 141-146.

[10]CHEN N X. A new method for inverse black body radiation problem [J].ChinesePhysLett， 1987，4(8)： 337-340.

[11]JI F M， DAI X X. A new solution method for black-body radiation inversion and the solar area-temperature distribution [J].SciChina-PhysMechAstron， 2011，54(11)： 2097-2102.

[12]WEN T， MING D M， DAI X X，etal. Type of inversion problem in physics： An inverse emissivity problem [J].PhysRevE：RapidCommunication， 2001，63(4)： 045601.

[13]JI F M， YE J P， SUN L，etal. An inverse transmissivity problem， its möbius inversion solution and new practical solution method [J].PhysLettA， 2006，352(4/5)： 467-472.

[14]MING D M， WEN T， DAI X X，etal. Generalized emissivity inverse problem [J].PhysRevE， 2002，65(4)： 045601(R).

[15]PHILLIPS D L. A technique for the numerical solution of certain integral equations of the first kind [J].JACM， 1962，9(9)： 84-97.

[16]TIKHONOV A N. Solution of incorrectly formulated problems and the regularization method [J].SovMathDolk， 1963，4： 1035-1038.

[17]BAKER C T H， FOX L， MAYERS D F，etal. Numerical solution of Fredholm integral equation of the first kind [J].ComputJ， 1964，7： 141-148.

[18]LAKHTAKIA A. The inverse black body radiation problem [J].ModernPhysicsLetterB， 2011，5(5)： 491-497.

[19]RIEMANN B. On the number of prime numbers less than a given quantity [J].MonatsberBerlinerAkad， 1859，145： 671.

A New Method for Emissivity Inversion

JI Fengmin1， WANG Xiantong2， WAN Xing2， SHEN Jianmi2

(1.SchoolofMaterialsScienceandEngineering，NingboUniversityofTechnology，Ningbo315016，China；2.SchoolofScience，NingboUniversityofTechnology，Ningbo315211，China)

How to obtain the emissivity of objects is an important research content of the anti telemetry technology. The inverse emissivity problem provides us a convenient way to obtain the emissivity. The purpose of this inverse problem is to get the emissivity according to measured total radiant power， by solving an integral equation. Although an exact solution formula of this inverse problem was derived， it can never be used directly to obtain emissivity in practical inversion. The integral equation requires that the total radiant power must satisfy a so-called asymptotic behavior control condition. The experimental measured total radiant power is impossible to satisfy this ABC condition because of the ill-posed problem. In order to realize the inversion of real system， some concrete methods are necessary. At present， there is only one existing method， that is the universal function set method based on the Hermite functions. In this paper， starting from the famous Laguerre functions， an independent new method is developed for inverse emissivity problem. The new method has all the advantages of the old method. Moreover， numerical calculation shows that the new method has higher accuracy than the old one. It can calculate the theoretical spectrum with more complex structure. Each expansion term in the new method satisfies the asymptotic behavior control condition， so it can automatically filter out the effective information in the experimental total radiant power data. At the same time， the noise of information can be eliminated. This new method is suitable for the inversion of real system. If there are more sophisticated and perfect radiation power measurement data， this new method can be applied to the emissivity inversion for the first real system.

inverse problems in physics； emissivity； asymptotic behavior control condition； universal function set method

0427-7104(2016)04-0403-07

2015-11-13

国家自然科学基金(11447211，10675031)，浙江省新苗人才计划(2013R422034)

季丰民(1979-)，男，讲师，博士，E-mail： JiFengm@163.com.

O 411.1

A