放射源定位范例-基于極大似然估計和簡單蒙特卡洛方法的單放射源定位研究

基本介紹

本研究的目標(biāo)是展示極大似然估計和貝葉斯估計這兩種統(tǒng)計學(xué)方法在放射源定位中的應(yīng)用炉媒。在放射源定位研究中饱搏，放射源的位置僅僅是待估計的變量之一鞍帝，除此之外胎围，待定的變量還可以是放射源的數(shù)量束世、放射源的移動速度士嚎、環(huán)境材料等等背捌×ν迹考慮到本研究的目標(biāo)是展示統(tǒng)計學(xué)方法，如此多的待定變量只能增加所研究問題的復(fù)雜程度，同時針對每一個變量的統(tǒng)計過程都是一樣的显拜，變量數(shù)量的增加僅僅是提升了統(tǒng)計過程的冗余程度衡奥，對于入門選手容易迷失在多變量數(shù)學(xué)方程的構(gòu)造上，而不是統(tǒng)計學(xué)意義和估計流程的理解远荠。為此矮固，本研究將定位問題簡化，僅保留最根本的放射源坐標(biāo)譬淳，作為待估計變量档址，而將其他變量作為固定數(shù)值提供。同時邻梆，為了能更加簡明地體現(xiàn)數(shù)學(xué)思路守伸，將更具有實際意義的三維空間定位問題轉(zhuǎn)換為二維平面的定位問題，這樣處理的好處是將笛卡爾坐標(biāo)系下的放射源坐標(biāo)由三個獨立分量降低為兩個獨立分量浦妄，進一步簡化數(shù)學(xué)結(jié)構(gòu)尼摹。

本研究以放射源位置為估計目標(biāo)，設(shè)定待定位的放射源數(shù)量為1剂娄，即僅存在一個放射源在環(huán)境中蠢涝，放射源的強度已知，放射源僅存在于二維平面內(nèi)阅懦，并且在觀測期間保持在固定位置惠赫，不發(fā)生移動。

本研究選擇統(tǒng)計學(xué)中的極大似然估計方法故黑，作為定位的數(shù)學(xué)模型儿咱，以最簡單的蒙特卡羅方法作為數(shù)學(xué)模型的求解方法。

本研究的流程由三個部分組成场晶，首先是統(tǒng)計學(xué)模型構(gòu)造混埠，接著是統(tǒng)計學(xué)模型求解，最后是方案的評價诗轻。第一部分所謂統(tǒng)計學(xué)模型钳宪，是一個以觀測結(jié)果為自變量的概率函數(shù)，這個模型通過觀測結(jié)果與研究目標(biāo)的物理關(guān)系構(gòu)造得到扳炬，在本研究背景下吏颖，觀測結(jié)果就是探測器的響應(yīng)數(shù)量，研究目標(biāo)就是放射源的位置恨樟，物理關(guān)系就是 $\gamma$ 光子被探測器觀測到的概率半醉。第二部分的實際內(nèi)容是提出方案、評價方案和篩選方案劝术，以得到最優(yōu)方案為結(jié)果缩多。第三部分是對最優(yōu)方案進行評價呆奕，目標(biāo)是評價本次研究是否可以接受。

統(tǒng)計學(xué)模型構(gòu)造

探測器響應(yīng)模型

$\gamma?$ 射線從放射源發(fā)出衬吆，到被探測器接收并被響應(yīng)梁钾，這個過程可以分為三部分，并統(tǒng)一使用統(tǒng)計學(xué)知識進行描述逊抡。第一部分是 $\gamma?$ 射線在飛行過程中的衰減過程姆泻，第二部分是 $\gamma?$ 進入探測器的有效探測面，第三部分是 $\gamma?$ 被探測器響應(yīng)冒嫡。

第一部分與介于放射源和探測器之間的介質(zhì)有關(guān)麦射，使用指數(shù)衰減模型進行描述。

設(shè)定灯谣，探測器總數(shù)為 $ND$ ，探測器的編號以 $ND$ 為標(biāo)志蛔琅，坐標(biāo)使用二元向量描述胎许，第i個探測器的坐標(biāo)為 $\boldsymbol{r}_{Di}$ 。

設(shè)定罗售，放射源的總數(shù)為 $NS=1?$ 辜窑，放射源的編號以 $NS?$ 為標(biāo)志，第i個放射源的坐標(biāo)為 $\boldsymbol{r}_{Si}?$ 寨躁、放射性強度為 $I_{Si}?$ 穆碎。由于本研究僅存在一個放射源，因此這個放射源的編號為 $S0?$ 职恳，坐標(biāo)為 $\boldsymbol{r}_{S0}?$ 所禀，放射性強度為 $I_{S0}?$ 。這兩個參數(shù)中放钦， $I_{S0}?$ 為已知色徘，僅 $\boldsymbol{r}_{S0}?$ 未知，作為待估計變量操禀。

設(shè)定褂策， $\gamma?$ 射線在材料中的衰減過程為指數(shù)衰減，衰減系數(shù)為 $\mu({\boldsymbol{r})}?$ 颓屑，總宏觀界面(the total macroscopic cross sections of material)為 $\Sigma(\boldsymbol{r})?$ 斤寂， $\boldsymbol{r}?$ 表示某一指定位置，這樣設(shè)置是考慮到在放射源與探測器之間揪惦，存在各種介質(zhì)遍搞，比如空氣、混凝土墻壁器腋、樹木尾抑、土壤等等歇父，不同介質(zhì)對 $\gamma?$ 射線的衰減效應(yīng)不同。 $dis(\boldsymbol{r})?$ 表示向量 $\boldsymbol{r}?$ 的距離再愈，也就是 $\boldsymbol{r}?$ 的二范數(shù)榜苫，或者歐幾里得距離。
$I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \int_{Si}^{Di} \mu({\boldsymbol{r}}) \cdot \rho({\boldsymbol{r}}}) \cdot \rmtbbnnr7(dis(\boldsymbol{r})))$

$I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \int_{Si}^{Di} \Sigma({\boldsymbol{r}}}) \cdot \rmfvhn9ln(dis(\boldsymbol{r})))$

在本研究背景下翎冲，設(shè)定材料僅為空氣垂睬，并且不考慮 $\gamma$ 射線的能量對衰減系數(shù)的影響，也就是認(rèn)為 $\Sigma(\boldsymbol{r})?$ 已知抗悍，并且在任何位置均相等驹饺。
$I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \int_{Si}^{Di} \Sigma} \cdot \rmvjx9zz1(dis(\boldsymbol{r})))$

$I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \Sigma \cdot \int_{Si}^{Di} } \rmrn3lhhr(dis(\boldsymbol{r})))$

$I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \Sigma \cdot dis({\boldsymbol r}_{S0}-{\boldsymbol r}_{Di})})$

$I_{attenuated,Si,Di}= I_{Si} \cdot exp(- \Sigma \cdot \left\| \boldsymbol{r_{S0}}-\boldsymbol{r_{Di}} \right\|_2^2 )$

第二部分與探測器有效探測面積有關(guān)。默認(rèn)放射源發(fā)射粒子的方向是各向同性缴渊，在二維空間中赏壹，也就是 $2\pi?$ 方向發(fā)射。探測器與放射源距離為 $dis(Di,S0)?$ 衔沼，探測器的有效探測面積為 $A_{Di}?$ 蝌借，設(shè)定探測器的有效探測面的法線平行于探測器與放射源的連線 $\boldsymbol{r_{S0}}-\boldsymbol{r_{Di}}?$ ，則所有釋放出來的光子中指蚁，可以進入探測器有效探測面的數(shù)量 $I_{S0,Di}?$ 表示如下菩佑。
$I_{attenuatedAndGeo,S0,Di} = I_{attenuated,S0,Di} \cdot \frac{A_{Di}}{4\pi \left\| \boldsymbol{r_{S0}}-\boldsymbol{r_{Di}} \right\|_2^2}$
第三部分與探測器的固有探測效率有關(guān)。并不是所有進入探測器的光子都能被探測器探測凝化，對于一個能量的光子稍坯，探測器能夠探測到的光子數(shù)量比例是一個固定值，稱為固有探測效率搓劫，以 $\epsilon_{Di}?$ 表示瞧哟。被第i個探測器測量到的光子數(shù)量 $I_{Di}?$ 表示如下。
$I_{Di,S0,Expected} = I_{attenuatedAndGeo,S0,Di} \cdot \epsilon_{Di}$
綜合上面三個部分枪向，從S0號放射源釋放的 $\gamma?$ 光子绢涡，能夠被第i個探測器測量到的數(shù)量 $I_{Di}?$ 表示如下。
$I_{Di,S0,Expected} = I_{S0} \cdot exp(- \Sigma \cdot \left\| \boldsymbol{r_{S0}}-\boldsymbol{r_{Di}} \right\|_2^2)\cdot \frac{A_{Di}}{4\pi \left\| \boldsymbol{r_{S0}}-\boldsymbol{r_{Di}} \right\|_2^2} \cdot \epsilon_{Di}$

泊松分布和光子的觀測概率模型

探測器對單個光子的響應(yīng)過程遵循二項分布遣疯，而在觀測過程中雄可，探測器會對數(shù)量極大的光子群進行響應(yīng)，并且在本研究中認(rèn)為光子之間互相獨立缠犀，則整個觀測過程是遵循泊松分布的数苫。

第Di個探測器對第Si個放射源的響應(yīng)數(shù)量，也就是觀測到的 $\gamma?$ 光子數(shù)量為 $I_{Di,Si,Observed}?$ 辨液，而理論上應(yīng)該觀測到的數(shù)量為 $I_{Di,Si,Expected}?$ 虐急，這種事件的發(fā)生概率可以使用泊松分布進行計算。
$P(\boldsymbol r_{Si},I_{Si}|Obs_{Di}) \equiv Poisson(I_{Di,Si,Observed},I_{Di,Si,Expected})$

$P(\boldsymbol r_{Si},I_{Si}|Obs_{Di}) = \frac{{I_{Di,Si,Expected}}^{I_{Di,Si,Observed}}}{I_{Di,Si,Observed}!} \cdot e^{-I_{Di,Si,Expected}}$

單放射源多探測器觀測的概率模型

按照本研究的設(shè)定滔迈，環(huán)境中只有一個放射源止吁，探測器的數(shù)量不止一個被辑，并且分布于環(huán)境中的不同位置。假設(shè)經(jīng)過一次觀測敬惦，則每一臺探測器都有一個屬于自己的觀測數(shù)值盼理，也就是對 $\gamma?$ 光子的響應(yīng)次數(shù)。設(shè)定俄删，放射源編號S0宏怔；探測器編號從D0開始，至DN-1畴椰；編號為Di的探測器臊诊，所觀測到的響應(yīng)次數(shù)為 $I_{Di,S0,Observed}?$ ，理論上這臺探測器的響應(yīng)次數(shù)為 $I_{Di,S0,Expected}?$ 斜脂。對于這種設(shè)定下的觀測結(jié)果抓艳，可以將每一個探測器的觀測事件發(fā)生概率進行累積，并且以乘法的方式累積帚戳，用以表示本次觀測的總概率 $P(\boldsymbol r_{S0},I_{S0}|Obs_{All})?$ 玷或，其計算方法如下。
$P(\boldsymbol r_{S0},I_{S0}|Obs_{All}) = \prod_{Di=D0}^{DN-1} {P(\boldsymbol r_{S0},I_{S0}|Obs_{Di})}$

$P(\boldsymbol r_{S0},I_{S0}|Obs_{All}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Expected}}}$

似然函數(shù)的構(gòu)造

似然函數(shù)指以觀測數(shù)值為條件销斟，放射源的第Plani種分布情況為結(jié)果的概率，以 ${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})?$ 表示椒舵，可以發(fā)現(xiàn)蚂踊，似然函數(shù)就是上一節(jié)所描述的概率模型， $P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})?$ 笔宿。
${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) \equiv P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})$

${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Expected}}}$

基于極大似然估計的定位方法

假設(shè)目前存在PlanN種放射源位置的方案犁钟，則每一種方案都對應(yīng)一個似然函數(shù) ${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})?$ 。極大似然估計的思想是通過比較這些方案的似然函數(shù)泼橘，找出最大似然函數(shù)所對應(yīng)的方案涝动，就認(rèn)為該方案為放射源坐標(biāo)的解。
$\hat {\boldsymbol r}_{S0} \in \{ arg \ max \ {\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) \}$

統(tǒng)計學(xué)模型求解-方案的提出炬灭、評價和篩選

基本求解方法介紹

根據(jù)前文的分析醋粟，定位研究只剩下一個問題——方案的獲得，而獲取方案的方法有很多重归，為方便理解米愿，將方法分為數(shù)值方法和統(tǒng)計學(xué)方法，這種劃分方法不嚴(yán)謹(jǐn)鼻吮，后文可以看到育苟，一些具體算法同時采用了數(shù)值方法和統(tǒng)計學(xué)方法，如模擬退火椎木。

綜合上文介紹违柏，統(tǒng)計學(xué)模型的求解過程可以歸納為方案的獲取和方案的比較博烂。統(tǒng)計學(xué)方法首先構(gòu)造盡可能豐富的備選方案，形成一個可數(shù)的方案集合漱竖，再通過似然函數(shù)的計算和比較選取最合適的方案禽篱，也就是導(dǎo)致似然函數(shù)最大的方案。而數(shù)值方法則從似然函數(shù)出發(fā)闲孤，認(rèn)定最適合的方案存在于似然函數(shù)的極值處谆级，并且是在全局范圍內(nèi)致使似然函數(shù)最大的極值處，由此使用尋找極值的方法讼积。

數(shù)值方法將似然函數(shù)視作一個多自變量的函數(shù)肥照，在本研究背景下，自變量為放射源S0的坐標(biāo) ${\boldsymbol r}_{S0}$ 勤众。似然函數(shù)的最大值是其極值之一舆绎，只需要找到致使似然函數(shù)極值最大的自變量，就認(rèn)為是正確的解们颜。在這里吕朵，數(shù)值方法可以分成兩派，一派是首先找到所有極值窥突，再通過比較確定最大的極值努溃，也稱為全局最優(yōu)解，以模擬退火算法為代表阻问；另一派只尋找一個極值梧税，這個極值一般是最靠近自變量初始值的，以梯度下降法為代表称近。

統(tǒng)計學(xué)方法的基本思路是構(gòu)造所有可能的方案第队，在本研究背景下，就是構(gòu)造所有可能的位置刨秆，作為放射源位置的備選方案凳谦。在這種設(shè)計下，正確的放射源坐標(biāo)自然屬于構(gòu)造出的方案衡未，并且正確的方案對應(yīng)的似然函數(shù)必然最大尸执，所以只需要計算所有方案的似然函數(shù)，找到最大似然函數(shù)和其所對應(yīng)的方案缓醋，就可以認(rèn)為這個方案是正確的解剔交。

蒙特卡羅(MC)方法介紹

理論上，所有可能的位置方案有無限個改衩，如果通過列舉位置的方式岖常，是無法達(dá)到便利所有可能方案的。在研究中葫督，經(jīng)常采用隨機抽樣的方式來獲取備選方案竭鞍，雖然理論上獲得所有方案需要無限長的時間板惑，但是在真實求解中，往往認(rèn)為在有限時間內(nèi)獲得的方案即可以代表全部方案偎快，這種隨機采樣的方法稱為蒙特卡羅方法冯乘。

蒙特卡羅方法是一種基本的統(tǒng)計學(xué)方法，在其基礎(chǔ)上衍生出重要性采樣的蒙特卡羅方法晒夹、馬爾可夫鏈蒙特卡洛方法裆馒、模擬退火等方法。注意丐怯，雖然在上文介紹中喷好，將模擬退火方法作為數(shù)值方法介紹，但其算法流程中包含了蒙特卡羅的隨機采樣過程读跷，因此其本身并不能單純地歸類梗搅，而是一種綜合方法。

考慮本研究的初衷是示范放射源定位流程效览，為了確保模型構(gòu)造无切、模型求解和方案評價三個方面都能得到良好介紹，在各個方面的流程中都使用最簡潔的方法丐枉。由此哆键，最本質(zhì)的蒙特卡羅方法是適合使用的，因為最本質(zhì)的蒙特卡羅方法僅僅從隨機采樣出發(fā)瘦锹，不添加任何輔助操作籍嘹，保持了方法的簡明。

蒙特卡羅方法產(chǎn)生放射源位置方案

雖然本研究的目標(biāo)是估計放射源的位置沼本，但可以根據(jù)經(jīng)驗提出放射源位置的出現(xiàn)范圍噩峦。由于是經(jīng)驗性的锭沟，這個范圍可以很大抽兆，但這并不影響蒙特卡羅方法的使用。

使用笛卡爾坐標(biāo)系描述放射源所在平面族淮，以相互正交的X軸和Y軸組成坐標(biāo)系辫红，設(shè)定放射源的實際位置表示為 $\boldsymbol {r}_{real}(x_{real},y_{real})?$ ，經(jīng)驗性的范圍可以表示如下祝辣。
$x_{min} < x_{real} < x_{max}$

$y_{min} < y_{real} < y_{min}$

由于上文的經(jīng)驗性設(shè)置贴妻，方案的坐標(biāo) $\boldsymbol {r}_{Plani}(x_{Plani},y_{Plani})$ 必須滿足范圍限制。
$x_{Plani} \in \{ x| x_{min} < x < x_{max} \}$

$y_{Plani} \in \{ y| y_{min} < y < y_{max} \}$

$\boldsymbol {r}_{Plani}(x_{Plani},y_{Plani}) \in \{ \boldsymbol {r}(x,y) | x_{min} < x < x_{max} , y_{min} < y < y_{max}\}$

蒙特卡羅的思路是使用隨機數(shù)產(chǎn)生備選方案蝙斜，這里可以選用符合均勻分布的隨機數(shù)生成器分別為兩個分量制作隨機數(shù)名惩。
$x_{Plani} \sim Uniform(x_{min},x_{max})$

$y_{Plani} \sim Uniform(y_{min},y_{max})$

均勻分布的隨機數(shù)生成器是一個常用工具，在很多數(shù)學(xué)軟件和C/C++以及Python數(shù)學(xué)庫中都有提供孕荠，尤其是(0,1)均勻分布 $Uniform(0,1)$ 娩鹉。這里假設(shè)擁有(0,1)均勻分布的隨機數(shù)生成器 $Random::Uniform(0,1)$ 攻谁，演示如何生成指定區(qū)間的均勻分布。
$Random::Uniform(x_{min},x_{max}) = Random::Uniform(0,1) \times (x_{max}-x_{min}) + x_{min}$

$Random::Uniform(y_{min},y_{max}) = Random::Uniform(0,1) \times (y_{max}-y_{min}) + y_{min}$

至此弯予，已經(jīng)介紹了編號Plani的方案的生成方式戚宦，對于所有備選方案，通過重復(fù)該操作即可獲得锈嫩。由于隨機采樣的概率特性受楼，方案獲取的越多，越能接近真實的解呼寸。

方案評價和篩選

本研究的方案指放射源的坐標(biāo)艳汽，而且只有一個放射源，所以編號為Plani的方案是 $\boldsymbol{r}_{S0,Plaini}?$ 等舔。

方案的評價原理是依靠統(tǒng)計學(xué)模型骚灸，在本研究中就是一個泊松分布 $P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})?$ ，也就是前文所述的似然函數(shù) ${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})?$ 慌植，以方案中的放射源坐標(biāo)為期望甚牲，以探測器實際觀測得到的響應(yīng)次數(shù)為觀測結(jié)果，代入泊松分布并計算概率值 $P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})?$ 蝶柿，得到的概率值就是當(dāng)前方案的評價丈钙。

這里，再列寫一遍方案的概率計算方程交汤，實際就是似然函數(shù)的定義式雏赦，但是這里強調(diào)針對編號Plani的方案。
${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,Plani,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Plani,Expected}}}$

$I_{Di,S0,Plani,Expected} = I_{S0} \cdot exp(- \Sigma \cdot \left\| \boldsymbol{r_{S0,Plani}}-\boldsymbol{r_{Di}} \right\|_2^2)\cdot \frac{A_{Di}}{4\pi \left\| \boldsymbol{r_{S0,Plani}}-\boldsymbol{r_{Di}} \right\|_2^2} \cdot \epsilon_{Di}$

通過比較所有方案的概率值 ${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})$ 芙扎，選取最大概率值所對應(yīng)的方案星岗，假設(shè)為方案編號為PlanBest，則該方案所對應(yīng)的放射源坐標(biāo) $\boldsymbol{r}_{S0,PlainBest}?$ 就是最優(yōu)解戒洼。

最優(yōu)方案評價

方案的評價是本研究的最后一步工作俏橘，目標(biāo)是以最優(yōu)方案代表整個研究過程，以最優(yōu)方案的評價代表整個研究過程的評價圈浇，得出本次研究是否可以接受寥掐。

最優(yōu)方案的概率 $P(\boldsymbol r_{S0,PlainBest},I_{S0}|Obs_{All})$ 可以表示如下。
$P(\boldsymbol r_{S0,PlainBest},I_{S0}|Obs_{All}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,PlanBest,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,PlanBest,Expected}}}$

定位實驗

觀測數(shù)據(jù)的模擬和獲取

定位問題是以觀測數(shù)據(jù)為已知磷蜀，真實環(huán)境中召耘，觀測數(shù)據(jù)是指探測器的響應(yīng)次數(shù)，這是是通過實驗獲取的褐隆。但是我沒有探測器污它，也沒有放射源，就算有，我也懶得用衫贬。因此蜜宪，本研究使用模擬軟件Geant4對探測器的實驗過程進行模擬，并獲得觀測數(shù)據(jù)祥山。

設(shè)定放射源位于(1m, 2m, 0m)位置圃验，放射性強度為1e9，也就是在模擬中放出的 $\gamma$ 光子數(shù)缝呕，光子能量0.662MeV澳窑，發(fā)射方向在3D空間中各向同性。

設(shè)定使用HPGe探測器供常，其靈敏材料為Ge晶體摊聋，其幾何形狀為立方體，幾何尺寸為(5cm, 5cm, 8cm)栈暇。一共使用20個HPGe探測器麻裁，分別放置于20個位置，編號為Di的探測器的坐標(biāo) ${\boldsymbol r}_{Di}?$ 源祈。
${\boldsymbol r}_{Di} = [x_{Di}, y_{Di}, 0]^{T}$

$x_{Di} = -5m + (1m \times n), n = Di \ mod \ 10$

$y_{Di} = -3m, if \ \frac{Di}{10} = 0$

$y_{Di} = 3m, if \ \frac{Di}{10} = 1$

為了方便了解煎源，附上模擬程序中構(gòu)造探測器坐標(biāo)的代碼。

double tarlocationx_Detector = -5.*m;
double tarlocationy_Detector = -5.*m;
double tarlocationz_Detector = 0.;    

for(int i=0;i<10;i++)
{   
        tarlocationx_Detector += (1.*m);
        tarlocationy_Detector = -3.*m;
        std::string nameTemp1 = "Detector" + std::to_string(10*0+i);
  
        tarlocationy_Detector = 3.*m;
        std::string nameTemp2 = "Detector" + std::to_string(10*1+i);
}

蒙特卡羅模擬得到的是各個探測器內(nèi)光子發(fā)射事件的能量沉積過程香缺，再對這些能量沉積過程進行統(tǒng)計和篩選手销，最終得到每一個探測器的響應(yīng)次數(shù)，也就是觀測數(shù)據(jù)图张。

統(tǒng)計學(xué)模型構(gòu)造實例

探測器響應(yīng)模型實例

探測器響應(yīng)模型考慮三個物理過程锋拖，光子在材料中的衰減過程、探測器的有效探測面積造成的幾何探測效率祸轮、探測器的固有探測效率兽埃。這三個物理過程對應(yīng)計算探測器響應(yīng)次數(shù)期望值 $I_{Di,Si,Expected}$ 的三個部分。

首先討論第一部分适袜， $\gamma$ 射線在空氣中的衰減作用是光子與材料的相互作用的宏觀體現(xiàn)柄错。本研究所采用的光子能量為0.662MeV，這個能量的光子與空氣中的各種組成元素獨立發(fā)生反應(yīng)痪蝇，更加嚴(yán)謹(jǐn)?shù)卣f鄙陡，應(yīng)該是與空氣中的各種同位素獨立發(fā)生反應(yīng)冕房。反應(yīng)截面躏啰，也就是反應(yīng)概率，較高的反應(yīng)包括三種耙册，光電效應(yīng)给僵、康普頓散射、電子對效應(yīng)。其中帝际，電子對效應(yīng)需要光子能量在1.022MeV蔓同，也就是兩倍的電子質(zhì)量能量，以上才可能發(fā)生蹲诀，因此斑粱，本研究環(huán)境下，光子能量不足以發(fā)生電子對效應(yīng)脯爪，只能發(fā)生光電效應(yīng)和康普頓效應(yīng)则北。

由于空氣的組成成分，主要包括H痕慢、O尚揣、N，對0.662MeV光子而言掖举，發(fā)生光電效應(yīng)和康普頓散射的截面很低快骗，探測距離也沒有長到可以明顯衰減光子數(shù)量，由此塔次，在粗略估計時方篮，可以認(rèn)為不發(fā)生光電效應(yīng)和康普頓散射這些反應(yīng)，也就是無衰減励负，衰減系數(shù) $\Sigma$ 為0恭取。

實際上是我懶得再針對光子在空氣中的衰減系數(shù)進行模擬了，但如果最優(yōu)方案對觀測數(shù)據(jù)的符合程度不行熄守，那么還是要回來估計衰減系數(shù)蜈垮。

其次討論第二部分，幾何探測效率裕照。使用探測器在X和Z軸的面代表編號為Di的探測器的有效探測面攒发，則有效探測面積 $A_{Di}$ 為 $40cm^2$ 。由于設(shè)定中晋南，所有探測器都一致惠猿，因此每一個探測器的有效探測面積都是 $40cm^2$ 。
$A_{Di} = 5cm \times8cm=40cm^2$
放射源實際坐標(biāo) $\boldsymbol{r}_{S0,real}?$ 為 $(1m,2m,0)?$ 负间，放射源到編號為Di的探測器的距離 $dis(\boldsymbol{r}_{S0,real}-{\boldsymbol r}_{Di})?$ 偶妖，計算方程如下。
$dis(\boldsymbol{r}_{S0,real}-{\boldsymbol r}_{Di}) = \sqrt{(x_{S0,real}-x_{Di})^2+(y_{S0,real}-y_{Di})^2}$
使用二范數(shù)代表距離政溃，以簡化表達(dá)趾访。
$dis(\boldsymbol{r}_{S0,real}-{\boldsymbol r}_{Di}) = \left \| \boldsymbol{r}_{S0,real}-{\boldsymbol r}_{Di} \right \|_2^2$
但是，在極大似然估計中董虱，并不知道真實坐標(biāo)扼鞋，而是以每個方案下的假想坐標(biāo)代替真實坐標(biāo)申鱼，對編號為Plani的方案，其假想坐標(biāo)為 $\boldsymbol{r}_{S0,Plani}$ 云头，其距離 $dis(\boldsymbol{r}_{S0,Plani}-{\boldsymbol r}_{Di})$ 捐友。
$dis(\boldsymbol{r}_{S0,Plani}-{\boldsymbol r}_{Di}) = \left \| \boldsymbol{r}_{S0,Plani}-{\boldsymbol r}_{Di} \right \|_2^2$
最后討論第三部分，探測器固有探測效率 $\epsilon_{Di}?$ 溃槐。雖然固有探測效率與探測器接收到的光子能量有關(guān)匣砖，但是本研究簡化處理，設(shè)定固有探測效率為定值昏滴，具體設(shè)定為10%脆粥。

綜合上面三個部分，在Plani中影涉，從S0號放射源釋放的 $\gamma$ 光子变隔，能夠被第i個探測器測量到的數(shù)量 $I_{Di}$ 表示如下。
$I_{Di,S0,Expected,Plani} = I_{S0} \cdot \frac{A_{Di}}{4\pi \left\| {\boldsymbol{r}_{S0,Plani}}-{\boldsymbol{r}_{Di}} \right\|_2^2} \cdot \epsilon_{Di}$

觀測概率模型實例

以泊松分布 $Poisson()?$ 描述方案對觀測數(shù)據(jù)的符合程度蟹倾，對編號Plani的方案匣缘，其符合程度計算如下。
$P(\boldsymbol r_{S0,Plani},I_{S0}|Obs_{Di}) \equiv Poisson(I_{Di,S0,Observed},I_{Di,S0,Plani,Expected})$

$P(\boldsymbol r_{S0,Plani},I_{S0}|Obs_{Di}) = \frac{{I_{Di,S0,Plani,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Plani,Expected}}$

其中鲜棠， $I_{Di,S0,Observed}$ 是探測器響應(yīng)次數(shù)的觀測數(shù)據(jù)肌厨，在本研究中，就是通過Geant4蒙特卡羅模擬軟件獲得的豁陆； $I_{Di,S0,Plani,Expected}$ 是探測器響應(yīng)次數(shù)的期望值柑爸，根據(jù)真實的放射源坐標(biāo)計算得到。

單放射源似然函數(shù)實例

似然函數(shù)指以觀測數(shù)值為條件盒音，放射源的第Plani種分布情況為結(jié)果的概率表鳍，以 ${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})?$ 表示。
${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) \equiv P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})$

${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,Plani,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Plani,Expected}}}$

方案篩選實例

在本研究中，方案規(guī)定了放射源的坐標(biāo)，不同方案源葫，坐標(biāo)不同。但是厘熟，各個方案中，探測器的觀測數(shù)據(jù)都是相同的维哈，探測器的坐標(biāo)也是相同的绳姨。由此，不同方案的似然函數(shù)中阔挠，處于分母的階乘相同飘庄，因為這個階乘僅僅與觀測數(shù)據(jù)有關(guān)，在比較方案的符合程度時谒亦，可以省略分母部分竭宰，也就是階乘部分。

由于泊松分布是一個適合貼近冪函數(shù)形式的分布函數(shù)份招，這里的似然函數(shù)也是以累乘形式構(gòu)造切揭，整體上可以通過求對數(shù)的方法簡化函數(shù)，因此锁摔，構(gòu)造一個似然函數(shù)的自然對數(shù)廓旬， $ln {\cal L}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})?$ 。
$ln{\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \sum_{Di=D0}^{DN-1} (I_{Di,S0,Observed} \cdot ln(I_{Di,S0,Plani,Expected})-I_{Di,S0,Plani,Expected}-ln(I_{Di,S0,Observed}!))$

$ln{\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \sum_{Di=D0}^{DN-1} (I_{Di,S0,Observed} \cdot ln(I_{Di,S0,Plani,Expected})-I_{Di,S0,Plani,Expected}) -DN \cdot ln(I_{Di,S0,Observed}!)$

將階乘部分省略谐腰，得到似然函數(shù)的核心部分 $ln {\cal L}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})?$ 孕豹。
$ln{\cal{L}}_{Obs,Plani,Core}({\boldsymbol r}_{S0},I_{S0}) = \sum_{Di=D0}^{DN-1} (I_{Di,S0,Observed} \cdot ln(I_{Di,S0,Plani,Expected})-I_{Di,S0,Plani,Expected})$
這樣處理之后，可以避免似然函數(shù)的計算超過計算機存儲范圍十气。

最后編輯于：2019.04.08 07:52:28

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末励背，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子砸西，更是在濱河造成了極大的恐慌叶眉，老刑警劉巖，帶你破解...
沈念sama閱讀 219,039評論 6贊 508
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件芹枷，死亡現(xiàn)場離奇詭異衅疙，居然都是意外死亡，警方通過查閱死者的電腦和手機鸳慈，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,426評論 3贊 395
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門饱溢，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人走芋，你說我怎么就攤上這事绩郎。” “怎么了翁逞？”我有些...
開封第一講書人閱讀 165,417評論 0贊 356
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵嗽上，是天一觀的道長。經(jīng)常有香客問我熄攘，道長兽愤，這世上最難降的妖魔是什么？我笑而不...
開封第一講書人閱讀 58,868評論 1贊 295
?港島之戀（遺憾婚禮）
正文為了忘掉前任挪圾，我火速辦了婚禮浅萧，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘哲思。我一直安慰自己洼畅，他們只是感情好，可當(dāng)我...
茶點故事閱讀 67,892評論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布棚赔。她就那樣靜靜地躺著帝簇，像睡著了一般徘郭。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上丧肴，一...
開封第一講書人閱讀 51,692評論 1贊 305
城市分裂傳說
那天残揉，我揣著相機與錄音，去河邊找鬼芋浮。笑死抱环，一個胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的纸巷。我是一名探鬼主播镇草，決...
沈念sama閱讀 40,416評論 3贊 419
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼瘤旨！你這毒婦竟也來了梯啤？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 39,326評論 0贊 276
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤存哲，失蹤者是張志新（化名）和其女友劉穎条辟，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體宏胯，經(jīng)...
沈念sama閱讀 45,782評論 1贊 316
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡羽嫡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 37,957評論 3贊 337
?白月光啟示錄
正文我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了肩袍。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片杭棵。...
茶點故事閱讀 40,102評論 1贊 350
活死人
序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖氛赐，靈堂內(nèi)的尸體忽然破棺而出魂爪，到底是詐尸還是另有隱情，我是刑警寧澤艰管，帶...
沈念sama閱讀 35,790評論 5贊 346
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布滓侍，位于F島的核電站，受9級特大地震影響牲芋，放射性物質(zhì)發(fā)生泄漏撩笆。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點故事閱讀 41,442評論 3贊 331
男人毒藥：我在死后第九天來索命
文/蒙蒙一缸浦、第九天我趴在偏房一處隱蔽的房頂上張望夕冲。院中可真熱鬧，春花似錦裂逐、人聲如沸歹鱼。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,996評論 0贊 22
一樁弒父案卜高，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽弥姻。三九已至南片，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間庭敦，已是汗流浹背疼进。一陣腳步聲響...
開封第一講書人閱讀 33,113評論 1贊 272
情欲美人皮
我被黑心中介騙來泰國打工，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留螺捐，地道東北人颠悬。一個月前我還...
沈念sama閱讀 48,332評論 3贊 373
代替公主和親
正文我出身青樓矮燎，卻偏偏與公主長得像定血，于是被迫代替她去往敵國和親。傳聞我的和親對象是個殘疾皇子诞外，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 45,044評論 2贊 355

放射源定位范例-基于極大似然估計和簡單蒙特卡洛方法的單放射源定位研究

放射源定位范例-基于極大似然估計和簡單蒙特卡洛方法的單放射源定位研究

基本介紹

統(tǒng)計學(xué)模型構(gòu)造

探測器響應(yīng)模型

泊松分布和光子的觀測概率模型

單放射源多探測器觀測的概率模型

似然函數(shù)的構(gòu)造

基于極大似然估計的定位方法

統(tǒng)計學(xué)模型求解-方案的提出炬灭、評價和篩選

基本求解方法介紹

蒙特卡羅(MC)方法介紹

蒙特卡羅方法產(chǎn)生放射源位置方案

方案評價和篩選

最優(yōu)方案評價

定位實驗

觀測數(shù)據(jù)的模擬和獲取

統(tǒng)計學(xué)模型構(gòu)造實例

探測器響應(yīng)模型實例

觀測概率模型實例

單放射源似然函數(shù)實例

方案篩選實例

推薦閱讀更多精彩內(nèi)容