基于三角形收縮法改進(jìn)的分?jǐn)?shù)傅里葉變換估計(jì)牛頓環(huán)參數(shù)

武進(jìn)敏, 姜 盛, 魯溟峰, 沈德明, 范軍芳,李亞峰, 張 峰, 陶 然

(1.北京信息科技大學(xué) 自動(dòng)化學(xué)院，北京 100192；2.北京理工大學(xué) 信息與電子學(xué)院，北京 100081)

1 引 言

由于干涉測(cè)量技術(shù)具有較高的測(cè)量精度和靈活度[1,2],在光學(xué)測(cè)量領(lǐng)域得到較為廣泛的應(yīng)用[3～8],其測(cè)量參數(shù)包括測(cè)量元件曲率半徑、光源波長(zhǎng)、介質(zhì)折射率、物體面形等關(guān)鍵物理量[9～12]。由于干涉測(cè)量技術(shù)是以干涉條紋圖像作為結(jié)果輸出,因此干涉測(cè)量技術(shù)中的關(guān)鍵環(huán)節(jié)是分析處理所采集的干涉條紋圖。牛頓環(huán)干涉條紋是一種典型的干涉條紋圖,目前,有關(guān)學(xué)者已提出針對(duì)牛頓環(huán)條紋圖分析處理的許多方法[13～19],包括:基于條紋周期的分析方法,如數(shù)環(huán)法;基于傅里葉變換的分析方法,如定心輪廓法。此外,因?yàn)榕ｎD環(huán)條紋圖其相位具有二次多項(xiàng)式分布規(guī)律,此類條紋圖在信號(hào)處理領(lǐng)域也被稱為線性調(diào)頻信號(hào)(chirp信號(hào)),所以作為處理chirp信號(hào)的最優(yōu)方法,基于分?jǐn)?shù)傅里葉變換的分析方法于2013年由Lu M F等[20]提出。對(duì)于此類方法,重要環(huán)節(jié)是搜索匹配旋轉(zhuǎn)角,為了精確檢測(cè)匹配旋轉(zhuǎn)角,需要減小步長(zhǎng)遍歷搜索范圍內(nèi)的所有旋轉(zhuǎn)角,導(dǎo)致算法時(shí)間過(guò)長(zhǎng)。為了滿足工程實(shí)際需求,目前采用非均勻搜索[21]的優(yōu)化方法在保證精度的前提下可以縮短搜索時(shí)間,但其易受局部極大值的影響。

為了進(jìn)一步解決在較短時(shí)間內(nèi)精確檢測(cè)到匹配旋轉(zhuǎn)角,本文提出基于三角收縮的優(yōu)化方法,通過(guò)分析牛頓環(huán)干涉條紋在分?jǐn)?shù)傅里葉域的幅值最大值與相應(yīng)旋轉(zhuǎn)角分布規(guī)律,基于最小二乘擬合證明這一分布規(guī)律在匹配旋轉(zhuǎn)角附近具有局部軸對(duì)稱的特性,從而通過(guò)該分布的對(duì)稱軸來(lái)確定匹配旋轉(zhuǎn)角,因此,該方法無(wú)需遍歷搜索范圍內(nèi)的旋轉(zhuǎn)角,大大減少了時(shí)間消耗。

2 理論分析及原理

2.1 牛頓環(huán)分?jǐn)?shù)傅里葉域信號(hào)處理原理

牛頓環(huán)是一種典型的具有二次相位分布的干涉條紋圖,如圖1所示。

圖1 牛頓環(huán)干涉條紋

其強(qiáng)度分布可描述為

I(x,y)=I0+Acos[φ(x,y)]

(1)

式中:I0是條紋圖背景強(qiáng)度;A是條紋圖調(diào)制強(qiáng)度;φ(x,y)表示牛頓環(huán)干涉條紋的相位,

φ(x,y)=Kπ[(x-x0)2+(y-y0)2]+π

(2)

式中:K=2/(λR);λ表示入射光波長(zhǎng);R是被測(cè)透鏡曲率半徑;(x0,y0)為牛頓環(huán)環(huán)心坐標(biāo)。

對(duì)于牛頓環(huán)條紋圖,其行(或列)也是一維chirp信號(hào),數(shù)學(xué)表達(dá)式可以表示為

I(x)=rect(x/Xm)[I0+f(x)+f*(x)]

(3)

式中:*號(hào)表達(dá)共軛;rect(x/Xm)表示持續(xù)時(shí)間;而f(x)可以表達(dá)為

f(x)=(A/2)exp[jφ(x)]

=(A/2)exp[j(πKx2+2πfcenx+φy)]

(4)

式中:fcen=-Kx0=-2x0/(λR);φy是行(或列)的固定相位。

一維chirp信號(hào)的分?jǐn)?shù)傅里葉變換定義為

(5)

式中:ux為行信號(hào)在分?jǐn)?shù)傅里葉域中的頻率;分?jǐn)?shù)傅里葉的核函數(shù),可以描述為

Kα(ux,x)=

(6)

將核函數(shù)表達(dá)式代入公式(5)可得:

(7)

當(dāng)旋轉(zhuǎn)角滿足

cotα=-K

(8)

牛頓環(huán)條紋圖每一行(或列)的Fα(ux)表示為

sinc[πXm(uxcscα-fcen)]

(9)

因此,對(duì)于一維chirp信號(hào),當(dāng)旋轉(zhuǎn)角滿足式(8)時(shí),其分?jǐn)?shù)傅里葉域幅值分布如圖2所示。由圖2可知,|Fα(ux)|在fcen=uxcscαm位置處取得最大值,αm相應(yīng)的即為匹配旋轉(zhuǎn)角。從而根據(jù)被測(cè)透鏡曲率半徑與匹配旋轉(zhuǎn)角的關(guān)系,R可表示為

(10)

圖2 一維信號(hào)分?jǐn)?shù)傅里葉域幅值分布

另外,牛頓環(huán)環(huán)心位置也可由αm下分?jǐn)?shù)傅里葉域幅值最大值的位置計(jì)算得到

(11)

2.2 三角收縮優(yōu)化方法原理

由于在每個(gè)旋轉(zhuǎn)角下,牛頓環(huán)在分?jǐn)?shù)傅里葉域幅值均有唯一最大值,不同分?jǐn)?shù)域中幅值最大值與相應(yīng)旋轉(zhuǎn)角α之間的關(guān)系如圖3所示。圖3中最大值對(duì)應(yīng)的α為αm。

圖3 分?jǐn)?shù)傅里葉域幅值最大值與相應(yīng)旋轉(zhuǎn)角關(guān)系

由圖3可知,幅值最大值分布關(guān)于αm具有局部軸對(duì)稱特性,因此通過(guò)確定對(duì)稱軸,即可得匹配旋轉(zhuǎn)角αm,根據(jù)αm可得被測(cè)透鏡曲率半徑。

為了證明幅值最大值分布具有軸對(duì)稱性,定義新函數(shù)FM(α),表示每個(gè)旋轉(zhuǎn)角下對(duì)應(yīng)的分?jǐn)?shù)域幅值最大值,數(shù)學(xué)表達(dá)式為

FM(αm)=Max(Fα)

(12)

根據(jù)公式(9),公式(11)可以寫成

FM(R)=Max(FR)

(13)

現(xiàn)將公式(13)中R改為R-Rm,分?jǐn)?shù)傅里葉域幅值最大值與曲率半徑關(guān)系如圖4所示。

圖4 分?jǐn)?shù)傅里葉域幅值最大值與曲率半徑關(guān)系

圖4中Rm表示由αm計(jì)算的曲率半徑。相應(yīng)地,公式(13)可以寫為

=a0+F1+F2

(14)

式中:n∈Z;F1是FM(R-Rm)奇次項(xiàng)之和;F2是FM(R-Rm)偶次項(xiàng)和;a0代表常數(shù)項(xiàng)。

圖5為圖4表示的原始分布與其多項(xiàng)式擬合分布對(duì)比圖。

圖5 圖4表示的原始分布與其多項(xiàng)式擬合分布對(duì)比圖

公式(14)中等式左端描述的分布規(guī)律與右端多項(xiàng)式擬合分布規(guī)律對(duì)比結(jié)果如圖5(a)所示;多項(xiàng)式擬合分布基本符合原分布規(guī)律,奇次項(xiàng)之和如圖5(b)中紅色曲線所示,偶次項(xiàng)之和如圖5(b)中藍(lán)色曲線所示。

因此,根據(jù)最小二乘擬合原理,FM(R-Rm)可近似表達(dá)為

(15)

根據(jù)公式(15)可知,FM(R)關(guān)于Rm所在軸具有局部軸對(duì)稱性。

本文提出基于三角收縮優(yōu)化方法牛頓環(huán)參數(shù)估計(jì)效率,三角形收縮法原理如圖6所示,其處理流程如下:

圖6 三角形收縮法原理

(1) 設(shè)置初值

a,Fa=FM(a)

b,Fb=FM(b)

c=(a+b)/2,Fc=FM(c)

(16)

(2) 判斷Fa,Fb,Fc大小,確定對(duì)稱軸分布區(qū)間,

① 當(dāng)FaFc>Fb時(shí),如圖6(b)所示,舍棄區(qū)間(c,b],保留區(qū)間[a,c]。

② 當(dāng)FaFb時(shí),連接3點(diǎn)構(gòu)成三角形,定義角度β為兩向量夾角可以描述為

(17)

式中:β∈[0,π];O代表(a,Fa)到(b,Fb);P代表(c,(Fa+Fb)/2)到(c,Fc)。

如圖6(c)所示,當(dāng)β=π/2時(shí),在點(diǎn)c處與縱軸平行的線即為對(duì)稱軸,c即為被測(cè)透鏡的曲率半徑;反之,當(dāng)β≠π/2時(shí),對(duì)于這種情況,需修正角度β;當(dāng)β>π/2時(shí),如圖6(d)所示,舍棄區(qū)間(c,b],保留區(qū)間[a,c];當(dāng)β<π/2時(shí),如圖6(e)所示,舍棄區(qū)間[a,c),保留區(qū)間[c,b];直到β=π/2或達(dá)到精度要求停止計(jì)算,從而通過(guò)對(duì)稱軸的確定得到被測(cè)透鏡曲率半徑(或匹配旋轉(zhuǎn)角)。

最后,可計(jì)算匹配旋轉(zhuǎn)角對(duì)應(yīng)的分?jǐn)?shù)傅里葉變換,可以進(jìn)一步獲得包含在分?jǐn)?shù)域幅值最大值位置對(duì)應(yīng)的物理參數(shù)。

3 實(shí) 驗(yàn)

為了驗(yàn)證所提方法的有效性,在配有M1芯片的計(jì)算機(jī)上使用MATLAB R2021a根據(jù)公式(1)仿真牛頓環(huán)條紋圖進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)中牛頓環(huán)背景強(qiáng)度I0=2,調(diào)制強(qiáng)度A=2,入射波波長(zhǎng)λ=589 nm,被測(cè)透鏡曲率半徑為R=0.86 m,牛頓環(huán)環(huán)心坐標(biāo)設(shè)值為條紋圖中心位置。

對(duì)于不同像素尺寸的牛頓環(huán)條紋圖,用本文方法估計(jì)曲率半徑、環(huán)心坐標(biāo)以及所需時(shí)間結(jié)果如表1所示。當(dāng)條紋圖尺寸小于640×640 pixels時(shí),所提方法處理?xiàng)l紋圖所需時(shí)間<1 s,而隨著條紋圖尺寸的增大,由表1可知,基于本文方法的曲率半徑估計(jì)值相對(duì)誤差減小,參數(shù)估計(jì)的精度也逐步提高,這主要與條紋圖中包含的條紋級(jí)數(shù)相關(guān)。

表1 不同尺寸牛頓環(huán)干涉條紋參數(shù)估計(jì)

以下列尺寸條紋圖為例:240×240,480×480,720×720 pixels;條紋圖及分?jǐn)?shù)域幅值分布如圖7所示。沿著圖中箭頭所示方向條紋圖包含的級(jí)數(shù)增加,αm下分?jǐn)?shù)域幅值分布能量聚集性更加明顯,因此,包含在αm中的曲率半徑估計(jì)更為精確。

圖7 不同尺寸牛頓環(huán)條紋圖及其分?jǐn)?shù)傅里葉域幅值分布

本文方法估計(jì)曲率半徑的相對(duì)誤差與條紋圖級(jí)數(shù)的關(guān)系如圖8所示。當(dāng)條紋圖尺寸增加,基于本文方法的參數(shù)估計(jì)時(shí)間增加,但仍可滿足工程實(shí)際需求,以處理1 080×1 080 pixels的圖像為例,估計(jì)值相對(duì)誤差為0.001%,處理時(shí)間為3.31 s。

圖8 曲率半徑估計(jì)相對(duì)誤差與條紋級(jí)數(shù)關(guān)系

為了進(jìn)一步測(cè)試所提方法的抗噪能力,實(shí)驗(yàn)中通過(guò)對(duì)尺寸為720×720 pixels的牛頓環(huán)條紋圖添加不同信噪比(SNR)的高斯噪聲來(lái)驗(yàn)證。

為了保證相同的參數(shù)估計(jì)精度,傳統(tǒng)分?jǐn)?shù)傅里葉變換采用1×104步長(zhǎng)均勻搜索方式估計(jì)高斯噪聲污損牛頓環(huán)干涉條紋物理參數(shù),相應(yīng)的實(shí)驗(yàn)結(jié)果如表2所示。

表2 高斯噪聲污損下傳統(tǒng)分?jǐn)?shù)傅里葉變換估計(jì)牛頓環(huán)干涉條紋參數(shù)

本文方法估計(jì)高斯噪聲污損牛頓環(huán)干涉條紋物理參數(shù),實(shí)驗(yàn)結(jié)果如表3所示。

表3 高斯噪聲污損下三角形收縮法估計(jì)牛頓環(huán)干涉條紋參數(shù)

實(shí)驗(yàn)結(jié)果表明:傳統(tǒng)分?jǐn)?shù)傅里葉變換和本文方法都具有很好的抗噪能力,但在保證同樣估計(jì)精度的同時(shí),傳統(tǒng)分?jǐn)?shù)傅里葉變換平均用時(shí)約為884 s,本文三角形收縮法平均用時(shí)約1.28 s。

高斯噪聲污損牛頓環(huán)條紋圖及其分?jǐn)?shù)傅里葉域幅值分布如圖9所示,對(duì)于SNR=-15的噪聲污損條紋圖,在分?jǐn)?shù)域中其幅值分布仍可以實(shí)現(xiàn)很好的能量聚集性,相應(yīng)地,采用本文方法估計(jì)曲率半徑的相對(duì)誤差為0.069%,處理時(shí)間為1.26 s。

圖9 高斯噪聲污損牛頓環(huán)條紋圖及其分?jǐn)?shù)傅里葉域幅值分布

此外,通過(guò)處理實(shí)際條紋圖,本文方法的可行性也得到了驗(yàn)證,對(duì)于實(shí)際采集的640×480 pixels的牛頓環(huán)條紋圖如圖10所示,基于本文方法估計(jì)的曲率半徑相對(duì)誤差為0.067%,處理時(shí)間為0.64 s。

圖10 實(shí)際牛頓環(huán)干涉條紋

4 結(jié) 論

針對(duì)傳統(tǒng)分?jǐn)?shù)傅里葉變換處理牛頓環(huán)參數(shù)估計(jì)時(shí)處理速度較慢這一問(wèn)題,本文提出基于三角收縮優(yōu)化的分?jǐn)?shù)域處理方法。

通過(guò)建立不同旋轉(zhuǎn)角α下牛頓環(huán)條紋圖分?jǐn)?shù)域幅值最大值分布模型,基于此模型利用初值構(gòu)造三角形,通過(guò)判斷三角形形狀,從而確定α收縮區(qū)間并對(duì)三角形形狀進(jìn)行修正,直到所構(gòu)造的三角形為等腰三角形,進(jìn)而通過(guò)確定等腰三角形對(duì)稱軸所在位置確定匹配旋轉(zhuǎn)角αm,利用αm確定包含在其中的物理參數(shù),如本文實(shí)驗(yàn)中的被測(cè)透鏡曲率半徑。該方法估計(jì)的曲率半徑相對(duì)誤差仿真結(jié)果最小可達(dá)0.001%,但是在處理實(shí)際條紋圖時(shí),其相對(duì)誤差在1%左右,實(shí)際測(cè)量時(shí)限制該方法估計(jì)精度的因素包括:①因光照不均勻?qū)е聴l紋圖對(duì)比度下降引入的誤差;②測(cè)量過(guò)程中存在像差以及被測(cè)元件變形等引入的系統(tǒng)誤差;③實(shí)際條紋圖處理時(shí)初值選取的較大偏差導(dǎo)致三角形模型判斷失誤而引入的誤差。對(duì)上述誤差源進(jìn)行校準(zhǔn),可以進(jìn)一步提升本文方法在估計(jì)牛頓環(huán)參數(shù)時(shí)的精度。

實(shí)驗(yàn)結(jié)果表明:本文方法具有可行性,對(duì)于圖像尺寸小于640×640 pixels的條紋圖,處理所需時(shí)間<1 s,隨著圖像尺寸增加,條紋圖中包含的條紋數(shù)目增加,曲率半徑估計(jì)值相對(duì)誤差降低,而處理時(shí)間仍可滿足工程實(shí)際需求,以處理1 080×1 080 pixels的圖像為例,估計(jì)值相對(duì)誤差為0.001%,處理時(shí)間為3.31 s。該方法估計(jì)720×720 pixels高斯噪聲污損的牛頓環(huán)干涉條紋圖像平均用時(shí)1.28s,約為傳統(tǒng)分?jǐn)?shù)傅里葉變換用時(shí)的1/700。